Poincaré section for 3BP

Question

I try to plot the Poincaré section for the 3BP i use the following code that works very well:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Clear["Global`*"]
μ=0.001;
SuperStar[μ]=1-μ;
Subscript[r1, 1]=Sqrt[(Subscript[y1, 1][t]+μ)^2+Subscript[y1, 2][t]^2];
Subscript[r1, 2]=Sqrt[(Subscript[y1, 1][t]-SuperStar[μ])^2+Subscript[y1, 2][t]^2];
a=0.298;
b=0.24;
c=0.10;
d=0.40;
tin=0.0;
tfin=60.;
eqns1={{Derivative[1][Subscript[y1, 1]][t]==Subscript[y1, 3][t],Subscript[y1, 1][0]==a},
{Derivative[1][Subscript[y1, 2]][t]==Subscript[y1, 4][t],Subscript[y1, 2][0]==b},
{Derivative[1][Subscript[y1, 3]][t]==2 Subscript[y1, 4][t]+Subscript[y1, 1][t]-(SuperStar[μ] (Subscript[y1, 1][t]+μ))/
\!\(\*SubsuperscriptBox[\(r1\), \(1\), \(3\)]\)-(μ (Subscript[y1, 1][t]-SuperStar[μ]))/
\!\(\*SubsuperscriptBox[\(r1\), \(2\), \(3\)]\),Subscript[y1, 3][0]==c},
{Derivative[1][Subscript[y1, 4]][t]==-2 Subscript[y1, 3][t]+Subscript[y1, 2][t]-(SuperStar[μ] Subscript[y1, 2][t])/
\!\(\*SubsuperscriptBox[\(r1\), \(1\), \(3\)]\)-(μ Subscript[y1, 2][t])/
\!\(\*SubsuperscriptBox[\(r1\), \(2\), \(3\)]\),Subscript[y1, 4][0]==d}};
Orbit=NDSolve[
eqns1,
{Subscript[y1, 1],Subscript[y1, 2],Subscript[y1, 3],Subscript[y1, 4]},{t,tin,tfin},
Method->{StiffnessSwitching,Method->{ExplicitRungeKutta,Automatic}},
AccuracyGoal->8,PrecisionGoal->8,MaxSteps->10000];

myx[t_]:=Evaluate[Subscript[y1, 1][t]/. Orbit];
myy[t_]:=Evaluate[Subscript[y1, 2][t]/. Orbit];

myxplot=Plot[myx[t],{t,0,tfin},AxesLabel->{"time","x"}]
myyplot=Plot[myy[t],{t,0,tfin},AxesLabel->{"time","y"}]

pp1=ParametricPlot[
Evaluate[{Subscript[y1, 1][t],Subscript[y1, 2][t]}/. Orbit],{t,0,tfin},
AxesLabel->{"x","y"},PlotPoints->150]

(*pp2=ParametricPlot3D[
Evaluate[{Subscript[y1, 1][t],Subscript[y1, 2][t],Subscript[y1, 4][t]}/.\[VeryThinSpace]Orbit],{t,0,tfin},
AxesLabel\[Rule]{"x","y","vy"},PlotPoints\[Rule]5000,
DisplayFunction\[Rule]Identity]*)

pp2=ParametricPlot3D[
Evaluate[{Subscript[y1, 1][t],Subscript[y1, 2][t],Subscript[y1, 4][t]}/.Orbit],{t,0,tfin},AxesLabel->{"x","y","vy"},
PlotPoints->4000,MeshFunctions->{#3&},Mesh->{{0}},
MeshStyle->{Directive[PointSize[Large],Red]},
BoxRatios->{1,1,1},DisplayFunction->Identity]

pl1=Graphics3D[
{Green,Opacity[0.9],Polygon[{{-3,-3,0},{-3,3,0},{3,3,0},{3,-3,0}}]},
DisplayFunction->Identity];

Show[pp2,pl1,DisplayFunction->$DisplayFunction,ImageSize->{500,500}]

If you have some problem with subscript You can convert it in standard form, anyway the problem is in the following part of the code (thanks a lot to Belisarius for his very good tips! intersection between function and plane ) that show a section that is different from the one plotted in the last Show command. In fact the following command make a section on the, for istance, y-axis (the y1_2[t] variable) but i want the section on the z-axis (the y1_4[t] variable):

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[{Subscript[y1, 1][t],Subscript[y1, 2][t],Subscript[y1, 4][t]}/.Orbit],{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]

Can someone help me to modify this last part of the program to find the correct intersection.

Answer 1

Using the functions defined in my answer to your previous question here you have the intersections with all three coordinate planes:

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

SetAttributes[pp1, HoldAll];
pp1[u_] := Module[{gc, data, pp},
  pp = ParametricPlot3D[Evaluate[{x[t], y[t], p2[t]} /. sol2], {t, 0, tfin}, 
                       MeshFunctions -> {u &}, Mesh -> {{0}}, PlotStyle -> None, 
                       MeshStyle -> Red];
  gc = Cases[pp, GraphicsComplex[__], Infinity];
  data = Cases[Normal@gc, Point@x__ :> (x /. 0. :> Sequence[]), Infinity];
  Show[ListPlot@clusters@data, 
       ListLinePlot[#[[FindCurvePath[#][[1]]]] & /@ clusters[data]]]]

GraphicsRow@{pp1[#1], pp1[#2], pp1[#3]}

3D View (stealing code from Michael E2)

ParametricPlot3D[Evaluate[{x[t], y[t], p2[t]} /. sol2], {t, 0, tfin}, 
 PlotPoints -> 4000, MeshFunctions -> {#1 &, #2 &, #3 &}, 
 Mesh -> {{0}}, 
 MeshStyle -> {{Directive[PointSize[Medium], Red]}, {Directive[
     PointSize[Medium], Darker@Green]}, {Directive[PointSize[Medium], 
     Blue]}}, BoxRatios -> {1, 1, 1}, ViewPoint -> {1, 0, -2},
 PlotStyle -> {GrayLevel[.7]}]

BTW, you should try to understand the answers you receive here, not just use them.