(*integrate evolution dA/dt=DF(z(t))*A in a closed circular path around singularity t=t0:i*Pi/(2\nu) of the unperturbed solution z(t) *)
(Subscript[R, 1])[t_]:=(Sqrt[2])((Sinh[ν t])/((Cosh[ν t])^2))
(Subscript[R, 2])[t_]:=((Sqrt[2])/(Cosh[ν t]))
(Subscript[ψ, 2])[t_]:=t+β
(Subscript[ψ, 1])[t_]:=(Subscript[ψ, 2])[t]+π
(Subscript[x, 1])[t_]:=(Subscript[R, 1])[t] Cos[(Subscript[ψ, 1])[t]]
(Subscript[x, 2])[t_]:=(Subscript[R, 1])[t]Sin[(Subscript[ψ, 1])[t]]
(Subscript[y, 1])[t_]:=(Subscript[R, 2])[t] Cos[(Subscript[ψ, 2])[t]]
(Subscript[y, 2])[t_]:=(Subscript[R, 2])[t]Sin[(Subscript[ψ, 2])[t]]
z[t_]:={(Subscript[x, 1])[t],(Subscript[x, 2])[t],(Subscript[y, 1])[t],(Subscript[y, 2])[t]}
t0=((I  π)/(2 ν ));
t[τ_]:=-t0(1- (E^(I τ)))(*switch sign only when taking gamma 2, also to the derivative*)

M0={{0,-1,-ν,0},{1,0,0,-ν},{-ν,0,0,-1},{0,-ν,1,0}};
M1[t_]:={{0,0,2+Cos[2 (Subscript[ψ, 2])[t]] ,Sin[2 (Subscript[ψ, 2])[t]]},{0,0,Sin[2 (Subscript[ψ, 2])[t]],2-Cos[2 (Subscript[ψ, 2])[t]]},{0,0,0,0},{0,0,0,0}}
M2[t_]:={{0,0,20ϵ (Subscript[y, 1])[t]^3,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
DF[t_]:=M0+ν (((Subscript[R, 2])[t])^2)M1[t]+M2[t]
D[t[τ],τ]
MatrixForm[M0]//TraditionalForm
MatrixForm[M1[t]]//TraditionalForm
MatrixForm[M2[t]]//TraditionalForm
MatrixForm[DF[t]]//TraditionalForm
D[t[τ],τ]

ν=1/2;

Clear[A]

DA[τ_,A_List]:=-(E^(I τ))  π DF[t[τ]].Table[A[[i]][[j]],{i,4},{j,4}]
DA[τ,Table[(Subscript[a, i]^j),{i,4},{j,4}]]//MatrixForm//TraditionalForm

FullSimplify[DA[τ,Table[(Subscript[a, i]^j),{i,4},{j,4}]]]//MatrixForm//TraditionalForm

DARK[τ_,A_List,h_]:=Module[{k1,k2,k3,k4},
k1=DA[τ,A];
k2=DA[τ+h/2,A+k1 h/2];
k3=DA[τ+h/2,A+k2 h/2];
k4=DA[τ+h/2,A+k3 h];
h (k1+2k2+2k3+k4)/6
]

n=(10^5)(*;ρ=π;ν=0.5*);ϵ=(10^(-10));β=0;h=2.π/n;plotn=1000;
A=IdentityMatrix[4];
MatrixForm[DA[π,A]]
N[%]//MatrixForm
plotA=Table[A,{i,1,plotn}];

For[i=0,i<  n,i++,A+=DARK[h i,A,h];If[Mod[i plotn,n]== 0,plotA[[1+i  plotn /n]]=A]]
MatrixForm[A]
MatrixForm[IdentityMatrix[4]-A]

Beep[]

TraditionalForm[MatrixForm[IdentityMatrix[4]-A]]

Norm[IdentityMatrix[4]-A]

Export["D:\\UNIVERSITAT\\TFG\\numdata\\plotA_gamma_2_eps_1.dat",plotA]

1-Det[A]

Export["D:\\UNIVERSITAT\\TFG\\numdata\\A_gamma_2_eps_1.dat",A]

plotA[[400]]//MatrixForm

f=OpenRead["D:\\UNIVERSITAT\\TFG\\numdata\\plotA_gamma_2_eps_1.dat"];
plotA=ReadList[f,Table[Number,{i,1,32}]];
Close[f];
plotn=Length[plotA[[All,1]]];

Table[Table[ListLinePlot[{Table[plotA[[k]][[i]],{k,plotn}],Table[plotA[[k]][[i+1]],{k,plotn}]},ImageSize-> 200,PlotRange->All,PlotStyle-> Thick,PlotLegends->Placed[{Re,Im},{0.2,0.8}],Ticks->{None,Automatic}],{i,1+8 (j-1),8j,2}],{j,1,4}]//Grid
Table[Table[ListLinePlot[Table[{plotA[[k]][[i]],plotA[[k]][[i+1]]},{k,plotn}],PlotRange->All,PlotStyle->{Thick},AxesLabel->{Re,Im},AspectRatio->1],{i,1+8 (j-1),8j,2}],{j,1,4}]//Grid

Table[Table[ListLinePlot[{Table[plotA[[k]][[i]],{k,plotn}],Table[plotA[[k]][[i+1]],{k,plotn}]},PlotRange->All,PlotStyle->{Thick},PlotLegends->Placed[{Re,Im},{0.2,0.8}],Ticks->{None,Automatic}],{i,1+8 (j-1),8j,2}],{j,1,4}]//Grid

Table[Table[ListLinePlot[Table[{plotA[[k]][[i]],plotA[[k]][[i+1]]},{k,plotn}],PlotRange->All,PlotStyle->{Thick},AxesLabel->{Re,Im},AspectRatio->1],{i,1+8 (j-1),8j,2}],{j,1,4}]//Grid

Table[Table[ListPlot[Table[{plotA[[k]][[i]],plotA[[k]][[i+1]]},{k,plotn}],PlotStyle->{Thick},Ticks-> {None,None}],{i,1+8 (j-1),8j,2}],{j,1,4}]//Grid