(*Fig 2 A*)

sigma1=10;
sigma2=10;
taus=20;
psi0=9.64062;
psi1=9.64062;
sigmapunto2=sigma2/taus;
sigmapunto1=sigma1/taus;
nu1=1/(2*Pi*taus)*Exp[-psi0^2/(2*sigma1^2)];
nu2=1/(2*Pi*taus)*Exp[-psi1^2/(2*sigma2^2)];

c[tau_]=1/Cosh[tau/taus];

rho13[tau_,r_]=r*sigma1*sigma2*c[tau];

rho14[tau_,r_]=r*sigma1*sigma2*(-((Sech[(tau/taus)] Tanh[(tau/taus)])/taus));

rho23[tau_,r_]=-rho14[tau,r];

rho24[tau_,r_]=-r*sigma1*sigma2*(-(((Sech[(tau/taus)])^3)/(taus^2))+((Sech[(tau/taus)] ((Tanh[(tau/taus)])^2))/(taus^2)));
alpha1[tau_,r_]=(rho14[tau,r]*rho13[tau,r])/(sigma1^2*sigma2^2-rho13[tau,r]^2); 
beta1[tau_,r_]=((rho14[tau,r]*sigma2^2)/(sigma2^2*sigma1^2-rho13[tau,r]^2)); 
alpha2[tau_,r_]=((-rho14[tau,r]*sigma1^2)/(sigma1^2*sigma2^2-rho13[tau,r]^2)); 
beta2[tau_,r_]=((-rho14[tau,r]*rho13[tau,r])/(sigma2^2*sigma1^2-rho13[tau,r]^2));
sigmaepsilon1[tau_,r_]=(sigmapunto1^2-(alpha1[tau,r]^2*sigma1^2+2*alpha1[tau,r]*alpha2[tau,r]*rho13[tau,r]+alpha2[tau,r]^2*sigma2^2))^(1/2);
sigmaepsilon2[tau_,r_]=(sigmapunto2^2-(beta1[tau,r]^2*sigma1^2+2*beta1[tau,r]*beta2[tau,r]*rho13[tau,r]+beta2[tau,r]^2*sigma2^2))^(1/2);
Cov[tau_,r_]=rho24[tau,r]-(beta1[tau,r]*alpha1[tau,r]*sigma1^2+(alpha1[tau,r]*beta2[tau,r]+alpha2[tau,r]*beta1[tau,r])*rho13[tau,r]+alpha2[tau,r]*beta2[tau,r]*sigma2^2);
a[tau_,r_]=(-(alpha1[tau,r]*psi0+alpha2[tau,r]*psi1)/(sigmaepsilon1[tau,r]));
b[tau_,r_]=(-(beta1[tau,r]*psi0+beta2[tau,r]*psi1)/(sigmaepsilon2[tau,r]));
hermite[x_,n_]=(2^(-n/2)*HermiteH[n,x/Sqrt[2]]);
psibara[tau_,r_]=(1/2) Erfc[((a[tau,r])/(Sqrt[2]))];
psibarb[tau_,r_]=(1/2) Erfc[((b[tau,r])/(Sqrt[2]))];
factorial[x_]=x!;
phia[tau_,r_]=((E^(-(1/2) ((a[tau,r])^2)))/(Sqrt[2 π]));
phib[tau_,r_]=((E^(-(1/2) ((b[tau,r])^2)))/(Sqrt[2 π]));
valuen0[tau_,r_]=sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r]*((b[tau,r]*psibarb[tau,r])-phib[tau,r])*((a[tau,r]*psibara[tau,r])-phia[tau,r]);
valuen1[tau_,r_]=Cov[tau,r]*psibara[tau,r]*psibarb[tau,r];

Cab[tau_,n_,r_]=Which[n==0, valuen0[tau,r],n==1,valuen0[tau,r]+valuen1[tau,r], n≥2, valuen0[tau,r]+valuen1[tau,r]+Sum[(sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r])*phia[tau,r]*phib[tau,r]/i!*hermite[a[tau,r],i-2]*hermite[b[tau,r],i-2]*(Cov[tau,r]/(sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r]))^i,{i,2,n}]];

CorMatrix[tau_,r_]={{sigma1^2,r*sigma1*sigma2*c[tau]},{r*sigma1*sigma2*c[tau],sigma2^2}};
pr[tau_,r_]=Exp[-1/2*{psi0,psi1}.Inverse[CorMatrix[tau,r]].{psi0,psi1}]/(2*Pi*Sqrt[Det[CorMatrix[tau,r]]]);

R4=0.7;
Plot[{Cab[tau*taus,0,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4],Cab[tau*taus,1,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4],Cab[tau*taus,10,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5},PlotRange->All]
N0=Table[{tau,Cab[tau*taus,0,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5,0.05}];
N1=Table[{tau,Cab[tau*taus,1,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5,0.05}];
N10=Table[{tau,Cab[tau*taus,10,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5,0.05}];

NameFile=StringJoin[{"Path/Figures_Data/",ToString[N0_R0p7],".txt"}];
Export[NameFile,N0,"Table"];
NameFile=StringJoin[{"Path/Figures_Data/",ToString[N1_R0p7],".txt"}];
Export[NameFile,N1,"Table"];
NameFile=StringJoin[{"Path/Figures_Data/",ToString[N10_R0p7],".txt"}];
Export[NameFile,N10,"Table"];

(*Fig 2 B*)

sigma1=10;
sigma2=5;
taus=20;
psi0=9.64062;
psi1=9.64062;
sigmapunto2=sigma2/taus;
sigmapunto1=sigma1/taus;
nu1=1/(2*Pi*taus)*Exp[-psi0^2/(2*sigma1^2)];
nu2=1/(2*Pi*taus)*Exp[-psi1^2/(2*sigma2^2)];

c[tau_]=1/Cosh[tau/taus];

rho13[tau_,r_]=r*sigma1*sigma2*c[tau];

rho14[tau_,r_]=r*sigma1*sigma2*(-((Sech[(tau/taus)] Tanh[(tau/taus)])/taus));

rho23[tau_,r_]=-rho14[tau,r];

rho24[tau_,r_]=-r*sigma1*sigma2*(-(((Sech[(tau/taus)])^3)/(taus^2))+((Sech[(tau/taus)] ((Tanh[(tau/taus)])^2))/(taus^2)));
alpha1[tau_,r_]=(rho14[tau,r]*rho13[tau,r])/(sigma1^2*sigma2^2-rho13[tau,r]^2); 
beta1[tau_,r_]=((rho14[tau,r]*sigma2^2)/(sigma2^2*sigma1^2-rho13[tau,r]^2)); 
alpha2[tau_,r_]=((-rho14[tau,r]*sigma1^2)/(sigma1^2*sigma2^2-rho13[tau,r]^2)); 
beta2[tau_,r_]=((-rho14[tau,r]*rho13[tau,r])/(sigma2^2*sigma1^2-rho13[tau,r]^2));
sigmaepsilon1[tau_,r_]=(sigmapunto1^2-(alpha1[tau,r]^2*sigma1^2+2*alpha1[tau,r]*alpha2[tau,r]*rho13[tau,r]+alpha2[tau,r]^2*sigma2^2))^(1/2);
sigmaepsilon2[tau_,r_]=(sigmapunto2^2-(beta1[tau,r]^2*sigma1^2+2*beta1[tau,r]*beta2[tau,r]*rho13[tau,r]+beta2[tau,r]^2*sigma2^2))^(1/2);
Cov[tau_,r_]=rho24[tau,r]-(beta1[tau,r]*alpha1[tau,r]*sigma1^2+(alpha1[tau,r]*beta2[tau,r]+alpha2[tau,r]*beta1[tau,r])*rho13[tau,r]+alpha2[tau,r]*beta2[tau,r]*sigma2^2);
a[tau_,r_]=(-(alpha1[tau,r]*psi0+alpha2[tau,r]*psi1)/(sigmaepsilon1[tau,r]));
b[tau_,r_]=(-(beta1[tau,r]*psi0+beta2[tau,r]*psi1)/(sigmaepsilon2[tau,r]));
hermite[x_,n_]=(2^(-n/2)*HermiteH[n,x/Sqrt[2]]);
psibara[tau_,r_]=(1/2) Erfc[((a[tau,r])/(Sqrt[2]))];
psibarb[tau_,r_]=(1/2) Erfc[((b[tau,r])/(Sqrt[2]))];
factorial[x_]=x!;
phia[tau_,r_]=((E^(-(1/2) ((a[tau,r])^2)))/(Sqrt[2 π]));
phib[tau_,r_]=((E^(-(1/2) ((b[tau,r])^2)))/(Sqrt[2 π]));
valuen0[tau_,r_]=sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r]*((b[tau,r]*psibarb[tau,r])-phib[tau,r])*((a[tau,r]*psibara[tau,r])-phia[tau,r]);
valuen1[tau_,r_]=Cov[tau,r]*psibara[tau,r]*psibarb[tau,r];

Cab[tau_,n_,r_]=Which[n==0, valuen0[tau,r],n==1,valuen0[tau,r]+valuen1[tau,r], n≥2, valuen0[tau,r]+valuen1[tau,r]+Sum[(sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r])*phia[tau,r]*phib[tau,r]/i!*hermite[a[tau,r],i-2]*hermite[b[tau,r],i-2]*(Cov[tau,r]/(sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r]))^i,{i,2,n}]];

CorMatrix[tau_,r_]={{sigma1^2,r*sigma1*sigma2*c[tau]},{r*sigma1*sigma2*c[tau],sigma2^2}};
pr[tau_,r_]=Exp[-1/2*{psi0,psi1}.Inverse[CorMatrix[tau,r]].{psi0,psi1}]/(2*Pi*Sqrt[Det[CorMatrix[tau,r]]]);

R4=0.7;
Plot[{Cab[tau*taus,0,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4],Cab[tau*taus,1,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4],Cab[tau*taus,10,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5},PlotRange->All]

N0=Table[{tau,Cab[tau*taus,0,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5,0.05}];
N1=Table[{tau,Cab[tau*taus,1,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5,0.05}];
N10=Table[{tau,Cab[tau*taus,10,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5,0.05}];

NameFile=StringJoin[{"Path/Figures_Data/Asym",ToString[N0_R0p7],".txt"}];
Export[NameFile,N0,"Table"];
NameFile=StringJoin[{"Path/Figures_Data/Asym",ToString[N1_R0p7],".txt"}];
Export[NameFile,N1,"Table"];
NameFile=StringJoin[{"Path/Figures_Data/Asym",ToString[N10_R0p7],".txt"}];
Export[NameFile,N10,"Table"];

(*Fig 2 C*)

sigma1=10;
sigma2=10;
taus=20;
psi0=9.64062;
psi1=9.64062;
sigmapunto2=sigma2/taus;
sigmapunto1=sigma1/taus;
nu1=1/(2*Pi*taus)*Exp[-psi0^2/(2*sigma1^2)];
nu2=1/(2*Pi*taus)*Exp[-psi1^2/(2*sigma2^2)];

c[tau_]=1/Cosh[tau/taus];

rho13[tau_,r_]=r*sigma1*sigma2*c[tau];

rho14[tau_,r_]=r*sigma1*sigma2*(-((Sech[(tau/taus)] Tanh[(tau/taus)])/taus));

rho23[tau_,r_]=-rho14[tau,r];

rho24[tau_,r_]=-r*sigma1*sigma2*(-(((Sech[(tau/taus)])^3)/(taus^2))+((Sech[(tau/taus)] ((Tanh[(tau/taus)])^2))/(taus^2)));
alpha1[tau_,r_]=(rho14[tau,r]*rho13[tau,r])/(sigma1^2*sigma2^2-rho13[tau,r]^2); 
beta1[tau_,r_]=((rho14[tau,r]*sigma2^2)/(sigma2^2*sigma1^2-rho13[tau,r]^2)); 
alpha2[tau_,r_]=((-rho14[tau,r]*sigma1^2)/(sigma1^2*sigma2^2-rho13[tau,r]^2)); 
beta2[tau_,r_]=((-rho14[tau,r]*rho13[tau,r])/(sigma2^2*sigma1^2-rho13[tau,r]^2));
sigmaepsilon1[tau_,r_]=(sigmapunto1^2-(alpha1[tau,r]^2*sigma1^2+2*alpha1[tau,r]*alpha2[tau,r]*rho13[tau,r]+alpha2[tau,r]^2*sigma2^2))^(1/2);
sigmaepsilon2[tau_,r_]=(sigmapunto2^2-(beta1[tau,r]^2*sigma1^2+2*beta1[tau,r]*beta2[tau,r]*rho13[tau,r]+beta2[tau,r]^2*sigma2^2))^(1/2);
Cov[tau_,r_]=rho24[tau,r]-(beta1[tau,r]*alpha1[tau,r]*sigma1^2+(alpha1[tau,r]*beta2[tau,r]+alpha2[tau,r]*beta1[tau,r])*rho13[tau,r]+alpha2[tau,r]*beta2[tau,r]*sigma2^2);
a[tau_,r_]=(-(alpha1[tau,r]*psi0+alpha2[tau,r]*psi1)/(sigmaepsilon1[tau,r]));
b[tau_,r_]=(-(beta1[tau,r]*psi0+beta2[tau,r]*psi1)/(sigmaepsilon2[tau,r]));
hermite[x_,n_]=(2^(-n/2)*HermiteH[n,x/Sqrt[2]]);
psibara[tau_,r_]=(1/2) Erfc[((a[tau,r])/(Sqrt[2]))];
psibarb[tau_,r_]=(1/2) Erfc[((b[tau,r])/(Sqrt[2]))];
factorial[x_]=x!;
phia[tau_,r_]=((E^(-(1/2) ((a[tau,r])^2)))/(Sqrt[2 π]));
phib[tau_,r_]=((E^(-(1/2) ((b[tau,r])^2)))/(Sqrt[2 π]));
valuen0[tau_,r_]=sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r]*((b[tau,r]*psibarb[tau,r])-phib[tau,r])*((a[tau,r]*psibara[tau,r])-phia[tau,r]);
valuen1[tau_,r_]=Cov[tau,r]*psibara[tau,r]*psibarb[tau,r];

Cab[tau_,n_,r_]=Which[n==0, valuen0[tau,r],n==1,valuen0[tau,r]+valuen1[tau,r], n≥2, valuen0[tau,r]+valuen1[tau,r]+Sum[(sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r])*phia[tau,r]*phib[tau,r]/i!*hermite[a[tau,r],i-2]*hermite[b[tau,r],i-2]*(Cov[tau,r]/(sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r]))^i,{i,2,n}]];

CorMatrix[tau_,r_]={{sigma1^2,r*sigma1*sigma2*c[tau]},{r*sigma1*sigma2*c[tau],sigma2^2}};
pr[tau_,r_]=Exp[-1/2*{psi0,psi1}.Inverse[CorMatrix[tau,r]].{psi0,psi1}]/(2*Pi*Sqrt[Det[CorMatrix[tau,r]]]);

R1=0.1;
R2=0.4;
R3=0.6;
R4=0.7;
Plot[{Cab[tau*taus,10,R1]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R1],Cab[tau*taus,10,R2]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R2],Cab[tau*taus,10,R3]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R3],Cab[tau*taus,10,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5},PlotRange->All]

r1=Table[{tau,Cab[tau*taus,10,R1]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R1]},{tau,-5,5,0.05}];
r2=Table[{tau,Cab[tau*taus,10,R2]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R2]},{tau,-5,5,0.05}];
r3=Table[{tau,Cab[tau*taus,10,R3]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R3]},{tau,-5,5,0.05}];
r4=Table[{tau,Cab[tau*taus,10,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5,0.05}];

NameFile=StringJoin[{"Path/Figures_Data/",ToString[N10_R0p1],".txt"}];
Export[NameFile,r1,"Table"];
NameFile=StringJoin[{"Path/Figures_Data/",ToString[N10_R0p4],".txt"}];
Export[NameFile,r2,"Table"];
NameFile=StringJoin[{"Path/Figures_Data/",ToString[N10_R0p6],".txt"}];
Export[NameFile,r3,"Table"];
NameFile=StringJoin[{"Path/Figures_Data/",ToString[N10_R0p7],".txt"}];
Export[NameFile,r4,"Table"];

(*Figure 2 D*)

sigma1=10;
sigma2=5;
taus=20;
psi0=9.64062;
psi1=9.64062;
sigmapunto2=sigma2/taus;
sigmapunto1=sigma1/taus;
nu1=1/(2*Pi*taus)*Exp[-psi0^2/(2*sigma1^2)];
nu2=1/(2*Pi*taus)*Exp[-psi1^2/(2*sigma2^2)];

c[tau_]=1/Cosh[tau/taus];

rho13[tau_,r_]=r*sigma1*sigma2*c[tau];

rho14[tau_,r_]=r*sigma1*sigma2*(-((Sech[(tau/taus)] Tanh[(tau/taus)])/taus));

rho23[tau_,r_]=-rho14[tau,r];

rho24[tau_,r_]=-r*sigma1*sigma2*(-(((Sech[(tau/taus)])^3)/(taus^2))+((Sech[(tau/taus)] ((Tanh[(tau/taus)])^2))/(taus^2)));
alpha1[tau_,r_]=(rho14[tau,r]*rho13[tau,r])/(sigma1^2*sigma2^2-rho13[tau,r]^2); 
beta1[tau_,r_]=((rho14[tau,r]*sigma2^2)/(sigma2^2*sigma1^2-rho13[tau,r]^2)); 
alpha2[tau_,r_]=((-rho14[tau,r]*sigma1^2)/(sigma1^2*sigma2^2-rho13[tau,r]^2)); 
beta2[tau_,r_]=((-rho14[tau,r]*rho13[tau,r])/(sigma2^2*sigma1^2-rho13[tau,r]^2));
sigmaepsilon1[tau_,r_]=(sigmapunto1^2-(alpha1[tau,r]^2*sigma1^2+2*alpha1[tau,r]*alpha2[tau,r]*rho13[tau,r]+alpha2[tau,r]^2*sigma2^2))^(1/2);
sigmaepsilon2[tau_,r_]=(sigmapunto2^2-(beta1[tau,r]^2*sigma1^2+2*beta1[tau,r]*beta2[tau,r]*rho13[tau,r]+beta2[tau,r]^2*sigma2^2))^(1/2);
Cov[tau_,r_]=rho24[tau,r]-(beta1[tau,r]*alpha1[tau,r]*sigma1^2+(alpha1[tau,r]*beta2[tau,r]+alpha2[tau,r]*beta1[tau,r])*rho13[tau,r]+alpha2[tau,r]*beta2[tau,r]*sigma2^2);
a[tau_,r_]=(-(alpha1[tau,r]*psi0+alpha2[tau,r]*psi1)/(sigmaepsilon1[tau,r]));
b[tau_,r_]=(-(beta1[tau,r]*psi0+beta2[tau,r]*psi1)/(sigmaepsilon2[tau,r]));
hermite[x_,n_]=(2^(-n/2)*HermiteH[n,x/Sqrt[2]]);
psibara[tau_,r_]=(1/2) Erfc[((a[tau,r])/(Sqrt[2]))];
psibarb[tau_,r_]=(1/2) Erfc[((b[tau,r])/(Sqrt[2]))];
factorial[x_]=x!;
phia[tau_,r_]=((E^(-(1/2) ((a[tau,r])^2)))/(Sqrt[2 π]));
phib[tau_,r_]=((E^(-(1/2) ((b[tau,r])^2)))/(Sqrt[2 π]));
valuen0[tau_,r_]=sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r]*((b[tau,r]*psibarb[tau,r])-phib[tau,r])*((a[tau,r]*psibara[tau,r])-phia[tau,r]);
valuen1[tau_,r_]=Cov[tau,r]*psibara[tau,r]*psibarb[tau,r];

Cab[tau_,n_,r_]=Which[n==0, valuen0[tau,r],n==1,valuen0[tau,r]+valuen1[tau,r], n≥2, valuen0[tau,r]+valuen1[tau,r]+Sum[(sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r])*phia[tau,r]*phib[tau,r]/i!*hermite[a[tau,r],i-2]*hermite[b[tau,r],i-2]*(Cov[tau,r]/(sigmaepsilon2[tau,r]*sigmaepsilon1[tau,r]))^i,{i,2,n}]];

CorMatrix[tau_,r_]={{sigma1^2,r*sigma1*sigma2*c[tau]},{r*sigma1*sigma2*c[tau],sigma2^2}};
pr[tau_,r_]=Exp[-1/2*{psi0,psi1}.Inverse[CorMatrix[tau,r]].{psi0,psi1}]/(2*Pi*Sqrt[Det[CorMatrix[tau,r]]]);

R1=0.1;
R2=0.4;
R3=0.6;
R4=0.7;
Plot[{Cab[tau*taus,10,R1]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R1],Cab[tau*taus,10,R2]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R2],Cab[tau*taus,10,R3]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R3],Cab[tau*taus,10,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5},PlotRange->All]

r1=Table[{tau,Cab[tau*taus,10,R1]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R1]},{tau,-5,5,0.05}];
r2=Table[{tau,Cab[tau*taus,10,R2]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R2]},{tau,-5,5,0.05}];
r3=Table[{tau,Cab[tau*taus,10,R3]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R3]},{tau,-5,5,0.05}];
r4=Table[{tau,Cab[tau*taus,10,R4]*1000/(Sqrt[nu1*nu2])*pr[tau*taus,R4]},{tau,-5,5,0.05}];

NameFile=StringJoin[{"Path/Figures_Data/Asym",ToString[N10_R0p1],".txt"}];
Export[NameFile,r1,"Table"];
NameFile=StringJoin[{"Path/Figures_Data/Asym",ToString[N10_R0p4],".txt"}];
Export[NameFile,r2,"Table"];
NameFile=StringJoin[{"Path/Figures_Data/Asym",ToString[N10_R0p6],".txt"}];
Export[NameFile,r3,"Table"];
NameFile=StringJoin[{"Path/Figures_Data/Asym",ToString[N10_R0p7],".txt"}];
Export[NameFile,r4,"Table"];