# Poincaré section for 3BP

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

Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
Clear["Global*"]
\[Mu]=0.001;
SuperStar[\[Mu]]=1-\[Mu];
Subscript[r1, 1]=Sqrt[(Subscript[y1, 1][t]+\[Mu])^2+Subscript[y1, 2][t]^2];
Subscript[r1, 2]=Sqrt[(Subscript[y1, 1][t]-SuperStar[\[Mu]])^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[\[Mu]] (Subscript[y1, 1][t]+\[Mu]))/
\!$$\*SubsuperscriptBox[\(r1$$, $$1$$, $$3$$]\)-(\[Mu] (Subscript[y1, 1][t]-SuperStar[\[Mu]]))/
\!$$\*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[\[Mu]] Subscript[y1, 2][t])/
\!$$\*SubsuperscriptBox[\(r1$$, $$1$$, $$3$$]\)-(\[Mu] 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.

-
Somebody already told you not to use Subscripted symbols in this site. Please edit your code accordingly, that makes the code too difficult to read –  belisarius Mar 4 at 20:25

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.

-