Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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`*"]
\[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.

share|improve this question
3  
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

1 Answer 1

up vote 6 down vote accepted

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

Mathematica graphics

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

Mathematica graphics

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.