intersection between function and plane

Question

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}]

Can someone help me please?

Answer 1

Use Mesh and MeshFunctions:

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]

Answer 2

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]]]