# intersection between function and plane

I try to solve a ODE system of equations and i want to plot the intersection of the 3D solutions with a costant plane.

To do this i use the follow code:

Clear["Global*"]
tfin = 500;
T = p1^2/2 + p2^2/2;
V = x^2/2 + y^2/2 + x^2 y - y^3/3;
H = T + V;
eq1 = x'[t] == D[H, p1] /. p1 -> p1[t];
eq2 = y'[t] == D[H, p2] /. p2 -> p2[t];
eq3 = p1'[t] == -D[H, x] /. {x -> x[t], y -> y[t]};
eq4 = p2'[t] == -D[H, y] /. {x -> x[t], y -> y[t]};
energy = (H /. {x -> -.1, y -> -.2, p2 -> -.05}) == 0.06;
sol = Solve[energy, p1];
p10 = p1 /. sol[[2]];
sol2 = NDSolve[{eq1, eq2, eq3, eq4,
x[0] == -.1,
y[0] == -.2,
p1[0] == p10,
p2[0] == -.05},
{x[t], y[t], p1[t], p2[t]}, {t, 0, tfin}, MaxSteps -> 5000];
a = ParametricPlot[Evaluate[{x[t], y[t]} /. sol2], {t, 0, 20}];
b = ParametricPlot[Evaluate[{x[t], p1[t]} /. sol2], {t, 0, 20}];
c = ParametricPlot[Evaluate[{y[t], p2[t]} /. sol2], {t, 0, 20}];

gr1 = ParametricPlot3D[Evaluate[{x[t], y[t], p2[t]} /. sol2],
{t, 0, tfin},
PlotPoints -> 4000, BoxRatios -> {1, 1, 1}, ViewPoint -> {1, 0, -2},
DisplayFunction -> Identity];


This is the 'plane' of intersection. Maybe I must use something like PLot3D[0.,{y,ymin,ymax},{z,zmin,zmax}] because i want a plane that intersect the orbit in the plane x=0.

gr2 = Graphics3D[{Green, Opacity[0.9],
Polygon[{{0, -.31, -.31}, {0, -.35, .35}, {0, .35, .35}, {0, .35, \
-.31}}]}, DisplayFunction -> Identity];


This is the 'intersection' i only plot the solutions with a very thin x-axis but this is not a plot of the points of intersection between the solutions and the plane (pl1):

ParametricPlot3D[Evaluate[{x[t], y[t], p2[t]} /. sol2], {t, 0, tfin},
PlotPoints -> 4000, PlotStyle -> Directive[Red, Thick],
PlotRange -> {{0, 0.002}, {-.4, .4}, {-.4, .4}},
ViewPoint -> {1, 0, 0}, AxesLabel -> {"", "y", "p2"},
ImageSize -> {500, 500}]

Show[gr1, gr2, DisplayFunction -> \$DisplayFunction,
ImageSize -> {500, 500}]


-

gr1 = ParametricPlot3D[
Evaluate[{x[t], y[t], p2[t]} /. sol2], {t, 0, tfin},
PlotPoints -> 4000, MeshFunctions -> {#1 &}, Mesh -> {{0}},
MeshStyle -> {Directive[PointSize[Medium], Red]},
BoxRatios -> {1, 1, 1}, ViewPoint -> {1, 0, -2},
DisplayFunction -> Identity]


-
getOneCluster[pts_List, maxDist_?NumericQ] :=(*Returns a cluster*)
Module[{f},
f = Nearest[pts];
FixedPoint[Union@Flatten[f[#, {Infinity, maxDist}] & /@ #, 1] &, {First@pts}]]
clusters[data_] :=
Module[{f, dist},(*Some Characteristic Distance, assuming no isolated points*)
f = Nearest[data];
dist = 3 Max[EuclideanDistance[Last@f[#, 2], #] & /@ data];
Flatten[ Reap[NestWhile[Complement[#,
Sow@getOneCluster[#, dist]] &,  data, # != {} &]][[2]], 1]]

p1 = ParametricPlot3D[ Evaluate[{x[t], y[t], p2[t]} /. sol2], {t, 0, tfin},
MeshFunctions -> {#1 &}, Mesh -> {{0}}, PlotStyle -> None,
MeshStyle -> Red];
gc = Cases[p1, GraphicsComplex[__], Infinity];
data = First[(Normal@gc)[[1, 2, 1, 2 ;;]] /. Point -> Sequence][[All, 2 ;;]];
clusters[data];
Show[ListPlot@clusters@data, ListLinePlot[#[[FindCurvePath[#][[1]]]] & /@ clusters[data]]]


-
Thanks a lot, but if i want another section respect the one i found with mesh` how can i modify Your code? – Panichi Pattumeros PapaCastoro Mar 4 '14 at 20:13