ϵ=0;
G[{z1_,z2_,z3_,z4_}]:=z1 z3+z2 z4
F1[{z1_,z2_,z3_,z4_}]:=z1 z4-z2 z3
R1[{z1_,z2_,z3_,z4_}]:=z1 ^2+z2 ^2
R2[{z1_,z2_,z3_,z4_}]:=z3^2+z4 ^2
ψ1[{z1_,z2_,z3_,z4_}]:=If[z1 z2 == 0,1,ArcTan[z1,z2]]
ψ2[{z1_,z2_,z3_,z4_}]:=If[z3 z4 == 0,1,ArcTan[z3,z4]]
F2[{z1_,z2_,z3_,z4_}]:=(R1[{z1,z2,z3,z4}]-R2[{z1,z2,z3,z4}])/2+R2[{z1,z2,z3,z4}]^2/4
f[{z1_,z2_,z3_,z4_}]:={-z2+v z3 (z3^2+z4^2-1)+5 ϵ z3^4,z1+v z4 (z3^2+z4^2-1),-z4-v z1, z3-v z2}
(*runge kutta*)
ff[{z1_,z2_,z3_,z4_},h_]:=Module[{k1,k2,k3,k4},
k1=f[{z1,z2,z3,z4}];
k2=f[{z1,z2,z3,z4}+k1* h/2];
k3=f[{z1,z2,z3,z4}+k2* h/2];
k4=f[{z1,z2,z3,z4}+k3*h];
(k1+2k2+2k3+k4)/6
]

PS={Red,Green};PL={"Position","Velocity"};

(*stable evolution from z0*)

iter=700;n1=80;v=0.2;n1*iter
z0={0,0,1,0.9};
Z={z0}~Join~Table[{a,a,a,a},{i,2,n1 iter}];
For[i=2,i≤n1 iter,i++,Z[[i]]=Z[[i-1]]+ff[Z[[i-1]],1/iter]/iter]

Beep[]

Export["D:\\UNIVERSITAT\\TFG\\numdata\\evol_unpert_from_0_0_1_0dot9_unst.dat",Z]

ListPointPlot3D[Table[{Z[[i]][[1]],Z[[i]][[2]],Z[[i]][[3]]},{i,1,n1 iter,25}], PlotRange->All,AspectRatio->1,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])","(Subscript[x, 4])"}]

ExtendQuadrant[L_List]:=Join[L,L.{{-1,0},{0,1}},L.{{-1,0},{0,-1}},L.{{1,0},{0,-1}}]

ListPlot[ExtendQuadrant[Table[{(Sqrt[R2[Z[[i]]]])(-1),(Sqrt[R1[Z[[i]]]]) (-1)},{i,1,30 iter}]],PlotRange->All,AspectRatio->1,PlotStyle->PointSize[0.0001],AxesLabel->{"(Subscript[R, 2])","(Subscript[R, 1])"},AxesOrigin->{0,0}]

ListPlot[Table[{Z[[i]][[1]],Z[[i]][[2]]},{i,1,n1 iter}], PlotRange->All,AspectRatio->1,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotStyle->PS[[1]]]

f1=OpenRead["D:\\UNIVERSITAT\\TFG\\numdata\\evol_unpert_from_0_0_1_1_unst.dat"];
f0=OpenRead["D:\\UNIVERSITAT\\TFG\\numdata\\evol_unpert_from_0_0_1_0dot9_unst.dat"];
f2=OpenRead["D:\\UNIVERSITAT\\TFG\\numdata\\evol_unpert_from_0_0_1_1dot1_unst.dat"];
Z1=ReadList[f1,{Number,Number,Number,Number}];
Z0=ReadList[f0,{Number,Number,Number,Number}];
Z2=ReadList[f2,{Number,Number,Number,Number}];
Close[f0];Close[f1];Close[f2];
n1=Length[Z1[[All,1]]];n0=Length[Z0[[All,1]]];n2=Length[Z2[[All,1]]];

ListPlot[{ExtendQuadrant[Table[{(Sqrt[R2[Z0[[i]]]])(-1),(Sqrt[R1[Z0[[i]]]]) (-1)},{i,1,n0/2}]],ExtendQuadrant[Table[{(Sqrt[R2[Z1[[i]]]])(-1),(Sqrt[R1[Z1[[i]]]]) (-1)},{i,1,n1/2}]],ExtendQuadrant[Table[{(Sqrt[R2[Z2[[i]]]])(-1),(Sqrt[R1[Z2[[i]]]]) (-1)},{i,1,n2/2}]]},PlotRange->All,AspectRatio->1,AxesLabel->{Style["(Subscript[R, 2])",Large,Bold],Style["(Subscript[R, 1])",Large,Bold]},AxesOrigin->{0,0},ImageSize->Large,Ticks->{{-1,-(Sqrt[2]),1,(Sqrt[2])},Automatic},TicksStyle->Directive[Thick]]

ListPointPlot3D[Table[{Z2[[i]][[1]],Z2[[i]][[2]],Z2[[i]][[3]]},{i,1,n2,25}], PlotRange->All,AspectRatio->1,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])","(Subscript[x, 3])"}]

ListPlot[Table[{Z1[[i]][[1]],Z1[[i]][[2]]},{i,1,n1  0.35}], PlotRange->All,AspectRatio->1,PlotStyle->Red,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"}]
ListPlot[Table[{Z1[[i]][[3]],Z1[[i]][[4]]},{i,1,n1  0.35}], PlotRange->All,AspectRatio->1,AxesLabel->{"(Subscript[y, 1])","(Subscript[y, 2])"}]

f1=OpenRead["D:\\UNIVERSITAT\\TFG\\numdata\\unst_intersec_pert_bis.dat"];
f=OpenRead["D:\\UNIVERSITAT\\TFG\\numdata\\stab_intersec_pert_bis.dat"];
data=ReadList[f,{Number,Number,Number,Number}];
data1=ReadList[f1,{Number,Number,Number,Number}];
Close[f];Close[f1];
n=Length[data[[All,1]]];
Z1=data1;
Z=data;
Curve=Table[{ArcTan[Z[[i]][[4]],Z[[i]][[3]]],(Sqrt[Z[[i]][[3]]^2+Z[[i]][[4]]^2])-(Sqrt[2])(*Z[[i]][[3]]^2+Z[[i]][[4]]^2-2*)},{i,1,n}];
Curve1=Table[{ArcTan[Z1[[i]][[4]],Z1[[i]][[3]]],(Sqrt[Z1[[i]][[3]]^2+Z1[[i]][[4]]^2])-(Sqrt[2])(*Z1[[i]][[3]]^2+Z1[[i]][[4]]^2-2*)},{i,1,n}];
Curve=SortBy[Curve,First];
Curve1=SortBy[Curve1,First];

ListLinePlot[{Table[{Z[[i]][[1]],Z[[i]][[2]]},{i,1,n}],Table[{Z1[[i]][[1]],Z1[[i]][[2]]},{i,1,n}]}, PlotRange->All,AspectRatio->1,PlotLegends->Placed[{"stable","unstable"},{1,0.8}],PlotStyle->{Green,Red},AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotLabel->"Poincaré intersections"]

ListPlot[{Curve,Curve1},AspectRatio->1,PlotLegends->{"stable","unstable"},PlotStyle->{Red,Green}]
ListPlot[{Table[{Curve[[i]][[1]],Curve[[i]][[2]]-Curve1[[i]][[2]]},{i,1,n}]},AspectRatio->1,PlotLegends->{"differentce"},PlotStyle->{Red,Green}]

ListPointPlot3D[{Table[{Z1[[i]][[1]],Z1[[i]][[2]],Z1[[i]][[3]]},{i,1,n}]}, PlotRange->All,AspectRatio->1,PlotStyle->{Green}]
ListPointPlot3D[{Table[{Z[[i]][[1]],Z[[i]][[2]],Z[[i]][[4]]},{i,1,n}],Table[{Z1[[i]][[1]],Z1[[i]][[2]],Z1[[i]][[4]]},{i,1,n}]}, PlotRange->All,AspectRatio->1,PlotLegends->Placed[{"stable","unstable"},{1,0.8}],PlotStyle->{Green,Red}]

Curve1[[n]]

Curve1[[1]]

R=Interpolation[{{-2Pi,0}+Curve[[n]]}~Join~Curve~Join~{{2Pi,0}+Curve[[1]]}];
R1=Interpolation[{{-2Pi,0}+Curve1[[n]]}~Join~Curve1~Join~{{2Pi,0}+Curve1[[1]]}];

Plot[R1[t]-R[t],{t,-Pi,Pi},AxesLabel->{"θ","((Subscript[R, 1])^U)(θ=(Subscript[θ, 1]))-((Subscript[R, 1])^S)(θ=(Subscript[θ, 2]))"},Ticks->{Range[-Pi,Pi,Pi/2],Automatic}]

ListPlot[Abs[Fourier[Table[R[Pi t/512]-R1[Pi t/512],{t,-512,512}]]],Filling->Axis]

(*calcul monodromia: cas integrable*)