(* This program construct  Poincare sections with 
a required number of point for a given energy E *)

(*Hamiltonian and Hamilton equations of motion *)


h[x1_,x2_,p1_,p2_]:=(1/2)*(om1sq*x1^2+p1^2)+g*x1^4+
(1/2)*(om2sq*x2^2+p2^2)+g*x2^4-b*x1*x2;h[x1,x2,p1,p2]

g=1; b=0.5;om1sq=1.;om2sq=1.;
v1[t_]:=D[h[x1,x2,p1,p2],p1]/.{x1->x1[t],p1->p1[t],
x2->x2[t],p2->p2[t]};
v2[t_]:=D[h[x1,x2,p1,p2],p2]/.{x1->x1[t],p1->p1[t],
x2->x2[t],p2->p2[t]};
w1[t_]:=-D[h[x1,x2,p1,p2],x1]/.{x1->x1[t],p1->p1[t],
x2->x2[t],p2->p2[t]};
w2[t_]:=-D[h[x1,x2,p1,p2],x2]/.{x1->x1[t],p1->p1[t],
x2->x2[t],p2->p2[t]};

MakePS[ax1in_,ax2in_,ap1in_,ap2in_,npoints_,step_]:=Module[{acgoal,precgoal,sectone,nn,
timint, tmin, tmax,ditotal, aa,lm,xa,pa,di,sol},

acgoal=10;
precgoal=14;
sectone={};
x1in=ax1in;
x2in=ax2in;
p1in=ap1in;
p2in=ap2in;


enin=h[x1in,x2in,p1in,p2in];


(*npoints=required number of points in the section, 
ditotal =actual number of points in the section 
timint=interval of time used iteratively until..,*)


nn=0;timint=10;ditotal=0;
While[ditotal<npoints,
tmin=nn*timint;tmax=(nn+1)*timint;
ret={};
aa=NDSolve[{x1'[t]==v1[t],p1'[t]==w1[t],
x2'[t]==v2[t],p2'[t]==w2[t],x1[tmin]==x1in,p1[tmin]==p1in,
x2[tmin]==x2in,p2[tmin]==p2in},{x1,p1,x2,p2},
	{t,tmin,tmax},MaxSteps->70000,
	AccuracyGoal->acgoal,PrecisionGoal->precgoal];
	(* Poincare sections at x2=0 p2>0*)
 lm=Floor[(tmax-tmin)/step]-1;
Do[xa=(x2[tmin+l*step]/.aa)[[1]];xb=(x2[tmin+(l+1)*step]/.aa)[[1]];
	pa=(p2[tmin+l*step]/.aa)[[1]];temp1=Sign[xa]*Sign[xb];
	(* If the variable changes sign, we interpolate linearly
	to find the time to where the condition is satisfied *)
If[temp1<0 && pa>0,to=tmin+step*(l-xa/(xb-xa));
	ret=Union[ret,{to}];xto=(x2[to]/.aa)[[1]]]
		,{l,1, lm}];
di=Dimensions[ret][[1]];
sectone=Union[sectone,Table[{(x1[ret[[l]]]/.aa)[[1]],
(p1[ret[[l]]]/.aa)[[1]]},{l,1,di}]];

ditotal=ditotal+di; x1in=(x1[tmax]/.aa)[[1]];x2in=(x2[tmax]/.aa)[[1]];
p1in=(p1[tmax]/.aa)[[1]];p2in=(p2[tmax]/.aa)[[1]];
 nn=nn+1];Print[ListPlot[sectone,Frame->True,PlotRange->All,
FrameLabel->{"x2","p2",N[enin,10] ": Energy of 2 anh.", {k,rem}}]]]

ee=1;
Do[MakePS[0,0,ll*0.1,(2*ee-(ll*0.1)^2)^(1/2),200,0.01],{ll,0,2}]