Sharp Nekhoroshev estimates for the three-body problem around a periodic torus
Author: Santiago Barbieri, PhD student, Università degli studi di Padova, Italy

Initial data and definitions (Nasa J-2000)

(* Definition Astronomical Unit (in metres) *) 
UA=149597870700.;

(* Tuning on the mass of Jupiter *)
k=8.25;

(* Masses (in kg) *)
massesolaire=1.988500*10^30; (* Sun mass *)
m0=massesolaire*10^(00);     (* Central mass *)
m1=1.89819*10^27*10^(-k); (* Jupiter mass *)
m2=0.56834*10^27*10^(-k); (* Saturn mass *)
Ggrav=6.67259*10^(-11); (* Gravitational constant (in N*m^2/(kg^2)) *)

(* Jacobi masses *)
mu1=m0*m1/(m0+m1); 
mu2=(m0+m1)*m2/(m0+m1+m2);
M1=m0+m1;
M2=m0+m1+m2;
M2p=m0+m2; 

(* Variables appearing in the Kepler problem *)
massehkep1=Ggrav^2*mu1^3*M1^2;
massehkep2=Ggrav^2*mu2^3*M2^2;
sigma0=m0/(m0+m1);
sigma1=m1/(m0+m1);

(* Choice of the resonance *) 
p=5;
q=2;

(* Semi-major axes and periods *) (* !!! WARNING !!! The semi-major axes refer to ellipses which do NOT have common foci (Jacobi coordinates) *)
demigrandaxeP1= 5.20336301*UA*10^(00); (* Semi-major axis of Jupiter at J2000 *)
periodoj= (2*Pi*(demigrandaxeP1)^(3/2))/(Sqrt[Ggrav*M1]);  (* Period of Jupiter *)
Print["Test: period of Jupiter (y): ", periodoj/(3.1536*10^7)]; 
periodot=p*periodoj; (* Period of the corresponding torus *) 
T=p/q;   (* Frequency ratio for the two planets *) 
rapportdedistance=(T)^(2/3)*(M2/M1)^(1/3);  (* Distance ratio for the two planets *) 
demigrandaxeP2=demigrandaxeP1*rapportdedistance; (* Semi-major axis of Saturn *)
periodos= (2*Pi*(demigrandaxeP2)^(3/2))/(Sqrt[Ggrav*M2p]);  (* Period of Saturn *) 
Print["Test: period of Saturn (y): ", periodos/(3.1536*10^7)]; 
Print["Test: semi-major axis of Saturn (AU): ", demigrandaxeP2/(1.4960*10^11)];
Print["Test: periodos/periodoj= ", periodos/periodoj];

(* Poincaré resonant actions *)
Lambda1=mu1*Sqrt[Ggrav*M1*demigrandaxeP1];  (* Action 1 en variable de Poincare *)
Print["Lambda1= ", Lambda1];
Lambda2=mu2*Sqrt[Ggrav*M2*demigrandaxeP2]; (* Action 2 en variable de Poincare *)
Print["Lambda2= ", Lambda2];

Analyticity widths and maximal values for the eccentricities

(* Radius of the initial conditions for the Lambda variables (R in my article) *)
ri=0*10^(-8)*Min[Lambda1,Lambda2];
Print["R= ", ri];

(* Base of the polydisk in the Lambda variables (rho in my article) *)
runo=0*10^(-7)*Min[Lambda1,Lambda2];
Print["rho= ",runo];

(* Analyticity width in the Lambda variables (4r in my article) *) 
r=1.24*10^(-5)*Min[Lambda1,Lambda2];
Print["4r= ",r];

(* Analyticity width in the definition of anisotropic norms (r in my article) *) 
erre=r/4;
Print["r= ",erre];

(* Rough estimate of the final radius for the non normalized Lambda variables *)
(* I know that the normalized Lambda variables stay inside runo+erre and I know (a priori) that the size of the normal form and its inverse will be less than erre. I also need an estimate on the final radius in the actions in order to calculate the analiticity width of the cartesian variables, which is exactly this rf *) 
rf=runo+2*erre;
Print["rho+2r= ",rf];

(* Analyticity width for the angles (eccentric longitude). See Castan (2018) for more information *)
tinit=1.68*10^(-1);

(* Eccentricities at J2000 *) 
excjup=0.04839266;  (* Jupiter *)
Print["excjupinit= ", excjup];
excsat=0.05415060;   (* Saturn *)
Print["excsatinit= ", excsat];

(* Minimal value for the angular momentum *)
minmoment=(Lambda1-ri)*Sqrt[1-excjup^2]+(Lambda2-ri)*Sqrt[1-excsat^2];  

(* Maximal values for the eccentricities *)
excjupmax=Sqrt[1-((minmoment-(Lambda2+rf))/(Lambda1+rf))^2];
Print["excjupmax= ", excjupmax];
excsatmax=Sqrt[1-((minmoment-(Lambda1+rf))/(Lambda2+rf))^2];
Print["excsatmax= ", excsatmax];

(* Tuning on the analyticity width for the cartesian variables *) 
betap= 1.044; 

(* Analyticity width for the cartesian variables (4u in my article) *)
rhop=Sqrt[r*tinit]/betap;
Print["4u= ", rhop];

(* Radius of the real domain for the cartesian variables (xi in my article) *)
alphafnxy=3*10^(-3);
rho=Sqrt[Lambda1+Lambda2+2*rf-minmoment]-2*(rhop/4-alphafnxy*rhop/4);
Print["xi= ", rho];



(* Number of iterations for the iterative lemma *)
m=18;

Norm of the perturbation in the complex field (see Castan, 2018)

(* Integrable Kepler hamiltonian *)
hamkepP[Lambda10_,Lambda20_]=-massehkep1/(2*Lambda10^2)-massehkep2/(2*Lambda20^2);  
(* Unperturbed frequencies *)
freqkep1[Lambda_]=D[hamkepP[Lambda,x],Lambda];   
freqkep2[Lambda_]=D[hamkepP[x,Lambda],Lambda]; 
(* Second derivatives *)
derhamkepPxx[x_,y_]=D[D[hamkepP[x,y],x],x];  
derhamkepPyy[x_,y_]=D[D[hamkepP[x,y],y],y];  

(* Bound on the second derivatives of the unperturbed hamiltonian *)
maxderhamkepPxx[lambdauno_,rruno_,rr_]=derhamkepPxx[lambdauno-rruno-rr,0];
maxderhamkepPyy[lambdadue_,rruno_,rr_]=derhamkepPyy[0,lambdadue-rruno-rr];

(* Various norms (?) *)
norme1rDhamkepB[Lambda10_,Lambda20_,a_]=Abs[massehkep1/(Lambda10-a)^3+massehkep2/(Lambda20-a)^3]; 
normeDDhamkepB[Lambda10_,Lambda20_,a_]=Abs[3*(Ggrav^2*M1^2*mu1^3/(Lambda10-a)^4+Ggrav^2*M2^2*mu2^3/(Lambda20-a)^4)]; 
normeDDamelioree[Lambda10_,Lambda20_,a_]=Abs[massehkep1/(Lambda10-a)^3-massehkep1/(Lambda10)^3]+Abs[massehkep2/(Lambda20-a)^3-massehkep2/(Lambda20)^3];

(* Functions appearing when considering the complex Kepler equation and the angular analyticity width *)
a1[Lambda0_,r0_,rho0_,rhop0_,t0_]=1/Sqrt[Lambda0-r0];
a2[Lambda0_,r0_,rho0_,rhop0_,t0_]=Sqrt[1+3*rhop0^2/(2*(Lambda0-r0))];
a3[Lambda0_,r0_,rho0_,rhop0_,t0_]=(rho0+2rhop0)*Cosh[t0];
b1[Lambda0_,r0_,rho0_,rhop0_,t0_]=r0/(Sqrt[Lambda0-runo](Lambda0-r0));
b2[Lambda0_,r0_,rho0_,rhop0_,t0_]=(2rhop0*(Lambda0+r0)*(rhop0+Sqrt[2]rho0)+r0*(rho0+Sqrt[2]rhop0)^2)/(4(Lambda0-r0)^2);
b3[Lambda0_,r0_,rho0_,rhop0_,t0_]=rho0*Sinh[t0]+2rhop0*Cosh[t0];

lambdapmax0[Lambda0_,r0_,rho0_,rhop0_,t0_]=t0-(a1[Lambda0,r0,rho0,rhop0,t0]*a2[Lambda0,r0,rho0,rhop0,t0]*b3[Lambda0,r0,rho0,rhop0,t0]+a2[Lambda0,r0,rho0,rhop0,t0]*a3[Lambda0,r0,rho0,rhop0,t0]*b1[Lambda0,r0,rho0,rhop0,t0]+a1[Lambda0,r0,rho0,rhop0,t0]*a3[Lambda0,r0,rho0,rhop0,t0]*b2[Lambda0,r0,rho0,rhop0,t0]+b1[Lambda0,r0,rho0,rhop0,t0]*b2[Lambda0,r0,rho0,rhop0,t0]*b3[Lambda0,r0,rho0,rhop0,t0]);

li[Lambda0_,rin_,r0_,rho0_,rhop0_,ti0_]=1/Sqrt[Lambda0-rin-r0]*Sqrt[1+3*rhop0^2/(2*(Lambda0-rin-r0))]*(rho0+2*rhop0)*Cosh[ti0];

(* Definition of the analyticity width for the angle variable (4s in my article) *)
lambdapmax=Min[lambdapmax0[Lambda1,runo+r,rho,rhop,tinit],lambdapmax0[Lambda2,runo+r,rho,rhop,tinit]];
Print["4s= ",lambdapmax]; 

(* Inequalities which must be verified in order for the complex Kepler equation to induce a diffeomorphism *)
If[Not[lambdapmax>=0 &&
4*r<3(Lambda1-runo) &&
 4*r<3(Lambda2-runo) &&
(rho+Sqrt[2]*rhop)^2/(Lambda1-r)<3 &&(rho+Sqrt[2]*rhop)^2/(Lambda2-r)<3&& (*?????????*)
rho<Sqrt[2(Min[Lambda1,Lambda2]-runo-r)]&&
rhop<Sqrt[2(Min[Lambda1,Lambda2]-runo-r)]-rho&&
tinit<ArcCosh[Sqrt[(Lambda1-runo-r)/(1+3*rhop^2/(2(Lambda1-runo-r)))]/(rho+2rhop)]&&
tinit<ArcCosh[Sqrt[(Lambda2-runo-r)/(1+3*rhop^2/(2(Lambda2-runo-r)))]/(rho+2rhop)]],Print["ERROR! Problem with the Kepler equation"]];
If[Not[lambdapmax>0],Print["ERROR! The analyticity width for the angle variable is negative!"]];
If[Not[rho<Sqrt[2(Lambda1-runo-r)]],Print["ERROR! xi is too big!"]];

(*  New analyticity widths for the eccentric longitude (pag.19 of Castan article) *)

If[(a1[Lambda1,runo+r,rho,rhop,tinit]*a2[Lambda1,runo+r,rho,rhop,tinit]*b3[Lambda1,runo+r,rho,rhop,tinit]
	+a2[Lambda1,runo+r,rho,rhop,tinit]*a3[Lambda1,runo+r,rho,rhop,tinit]*b1[Lambda1,runo+r,rho,rhop,tinit]
	+a1[Lambda1,runo+r,rho,rhop,tinit]*a3[Lambda1,runo+r,rho,rhop,tinit]*b2[Lambda1,runo+r,rho,rhop,tinit]
	+b1[Lambda1,runo+r,rho,rhop,tinit]*b2[Lambda1,runo+r,rho,rhop,tinit]*b3[Lambda1,runo+r,rho,rhop,tinit])/lambdapmax<10^(-7),
	t1=lambdapmax,
	{t1}=
		Values[
			FindRoot[t0-(a1[Lambda1,runo+r,rho,rhop,t0]*a2[Lambda1,runo+r,rho,rhop,t0]*b3[Lambda1,runo+r,rho,rhop,t0]
				+a2[Lambda1,runo+r,rho,rhop,t0]*a3[Lambda1,runo+r,rho,rhop,t0]*b1[Lambda1,runo+r,rho,rhop,t0]
				+a1[Lambda1,runo+r,rho,rhop,t0]*a3[Lambda1,runo+r,rho,rhop,t0]*b2[Lambda1,runo+r,rho,rhop,t0]
				+b1[Lambda1,runo+r,rho,rhop,t0]*b2[Lambda1,runo+r,rho,rhop,t0]*b3[Lambda1,runo+r,rho,rhop,t0])-lambdapmax,{t0,0}]]];
				
If[
	(a1[Lambda2,runo+r,rho,rhop,tinit]*a2[Lambda2,runo+r,rho,rhop,tinit]*b3[Lambda2,runo+r,rho,rhop,tinit]
		+a2[Lambda2,runo+r,rho,rhop,tinit]*a3[Lambda2,runo+r,rho,rhop,tinit]*b1[Lambda2,runo+r,rho,rhop,tinit]
		+a1[Lambda2,runo+r,rho,rhop,tinit]*a3[Lambda2,runo+r,rho,rhop,tinit]*b2[Lambda2,runo+r,rho,rhop,tinit]
		+b1[Lambda2,runo+r,rho,rhop,tinit]*b2[Lambda2,runo+r,rho,rhop,tinit]*b3[Lambda2,runo+r,rho,rhop,tinit])/lambdapmax<10^(-7),
		t2=lambdapmax,
		{t2}=
			Values[
				FindRoot[t0-(a1[Lambda2,runo+r,rho,rhop,t0]*a2[Lambda2,runo+r,rho,rhop,t0]*b3[Lambda2,runo+r,rho,rhop,t0]
					+a2[Lambda2,runo+r,rho,rhop,t0]*a3[Lambda2,runo+r,rho,rhop,t0]*b1[Lambda2,runo+r,rho,rhop,t0]
					+a1[Lambda2,runo+r,rho,rhop,t0]*a3[Lambda2,runo+r,rho,rhop,t0]*b2[Lambda2,runo+r,rho,rhop,t0]
					+b1[Lambda2,runo+r,rho,rhop,t0]*b2[Lambda2,runo+r,rho,rhop,t0]*b3[Lambda2,runo+r,rho,rhop,t0])-lambdapmax,{t0,0}]]];
					
(* Analiticity width for the eccentric longitudes *)
Print["Analyticity width in eccentric longitude for each body : t1 = ",t1,"   t2 = ",t2];

(* Calculations on the cosine of the angle between the two planets - very complicated! See the article by Castan. *)

(* Term at the denominator *)
(* Formula (25) in Castan. Difference of the eccentricities *) (*??????????*)
diffdesexc[Lambdai_,R0_,r0_,rho0_,rhop0_]=4/3*r0/(Lambdai-R0)^(3/2)*(rho0+Sqrt[2]rhop0)+Sqrt[2]rhop0/Sqrt[Lambdai-R0]+8/3*rho0*rhop0*(rhop0+Sqrt[2]rho0)*Sqrt[Lambdai+R0]/(Lambdai-R0-r0)/(4Lambdai-rho0^2);

diffdeuxiemeterme[Lambda01_,Lambda02_,r0_,rho0_,rhop0_,t01_,t02_,e01_,e02_,rin_]=((e01*(Cosh[t01]-1)+diffdesexc[Lambda01,rin,r0,rho0,rhop0])/(1-e01)/(1-li[Lambda01,rin,r0,rho0,rhop0,t01])/(1-li[Lambda01,rin,r0,rho0,rhop0,t01])+(e02*(Cosh[t02]-1)+diffdesexc[Lambda02,rin,r0,rho0,rhop0])/(1-e01)/(1-e02)/(1-li[Lambda01,rin,r0,rho0,rhop0,t01]))*(1+e01)*(1+e02);

(* ????????? Shouldn't there be runo instead of ri ?????????*)
deuxiemeterme=diffdeuxiemeterme[Lambda1,Lambda2,runo,rho,rhop,t1,t2,excjupmax,excsatmax,ri];

(* Term at the numerator *)
difflambdamoinsw[Lambdai_,rin_,r0_,rho0_,rhop0_]=1/Sqrt[Lambdai-rin-r0]*Sqrt[1+3*rhop0^2/(2*(Lambdai-r0-rin))]*(rho0+2*rhop0);

alpha1[Lambdai_,rin_,r0_,rho0_,rhop0_]=difflambdamoinsw[Lambdai,rin,r0,rho0,rhop0]*(1+difflambdamoinsw[Lambdai,rin,r0,rho0,rhop0]/(2-3rhop0^2/(Lambdai-rin-r0)));

fonctiona[rho0_,rhop0_]=rho+Sqrt[2]rhop0;
fonctionb[Lambda0_,rin_,r0_]=1/Sqrt[Lambda0-rin-r0];
fonctionc[Lambda0_,rin_,r0_,rho0_,rhop0_]=1+difflambdamoinsw[Lambda0,rin,r0,rho0,rhop0]/(2-3rhop0^2/(Lambda0-rin-r0));
fonctionf[Lambda0_,rin_]=1/Sqrt[Lambda0-rin];
fonctiong[Lambda0_,rin_,r0_,rho0_,rhop0_]=(rho0+Sqrt[2]rhop0)^2*r0/8/(Lambda0-rin)/(Lambda0-rin-r0)+rhop0*(rhop0+2rho0)/4/(Lambda0-rin);
fonctioni[Lambda0_,rin_,r0_]=1/Sqrt[Lambda0-rin-r0]*r/2/(Lambda0-rin);

fonctionbeta2[Lambda0_,rin_,r0_,rho0_,rhop0_]=r/2/(Lambda0-rin)/Sqrt[Lambda0-rin-r0]*Sqrt[1+3*rhop0^2/(2*(Lambda0-r0-rin))]*(rho0+2*rhop0)+1/Sqrt[Lambda0-rin]*fonctiong[Lambda0,rin,r0,rho0,rhop0]*(rho0+2*rhop0)+2rhop0/Sqrt[Lambda0-rin];
fonctionj[Lambda0_,rin_,r0_,rho0_,rhop0_,t0_]=2r0/(4(Lambda0-rin-r0)-2rhop0^2)*difflambdamoinsw[Lambda0,rin,r0,rho0,rhop0]+difflambdamoinsw[Lambda0,rin,r0,rho0,rhop0]*(4r0+2rhop0*(rhop0+2rho0))/(4(Lambda0-rin-r0)-2rhop0^2)+fonctionbeta2[Lambda0,rin,r0,rho0,rhop0]*Cosh[t0]+rho/Sqrt[Lambda0]*(Cosh[t0]-1);

alpha2[Lambda0_,rin_,r0_,rho0_,rhop0_,t0_]=fonctiong[Lambda0,rin,r0,rho0,rhop0]*fonctiona[rho0,rhop0]*fonctionb[Lambda0,rin,r0]*fonctionc[Lambda0,rin,r0,rho0,rhop0]+2rhop0*fonctionb[Lambda0,rin,r0]*fonctionc[Lambda0,rin,r0,rho0,rhop0]+rho*fonctionc[Lambda0,rin,r0,rho0,rhop0]*fonctioni[Lambda0,rin,r0]+rho*fonctionf[Lambda0,rin]*fonctionj[Lambda0,rin,r0,rho0,rhop0,t0];

premierterme=Cosh[t1+t2]-1+Cosh[t2]*alpha2[Lambda1,runo,r,rho,rhop,t1]+alpha1[Lambda1,runo,r,rho,rhop]*(Cosh[t2]-1)+alpha2[Lambda2,runo,r,rho,rhop,t2]*Cosh[t1]+alpha1[Lambda2,runo,r,rho,rhop]*(Cosh[t1]-1)+alpha2[Lambda1,runo,r,rho,rhop,t1]*alpha1[Lambda2,runo,r,rho,rhop]+alpha2[Lambda2,runo,r,rho,rhop,t2]*alpha1[Lambda1,runo,r,rho,rhop];

(* Difference between real and complex cosine *)
xhi=premierterme+deuxiemeterme;
Print["Xhi = ",xhi];

(* Final bound *)

f0[theta_,beta0_,beta1_,xhireel_,xhicompl_]=1/sigma0/Sqrt[1-2*(xhireel+I*xhicompl)*beta0*Exp[I*theta]+(beta0*Exp[I*theta])^2]+1/sigma1/Sqrt[1+2*(xhireel+I*xhicompl)*beta1*Exp[I*theta]+(beta1*Exp[I*theta])^2]-(1/sigma0+1/sigma1);

(* Bound on 1/(Q_2-sigma0 Q_1) in the complex domain *)
betazero=sigma0*(Lambda1+runo+r)^2/(Lambda2-runo-r)^2*mu2^2*M2/mu1^2/M1*(1+li[Lambda1,runo,r,rho,rhop,t1])/(1-li[Lambda2,runo,r,rho,rhop,t2]);
fpremierterme=1/sigma0/Sqrt[1-2*(1+xhi)*betazero+(betazero)^2]-1/sigma0;
If[betazero>1+xhi-Sqrt[2xhi+xhi^2],Print[ " ERROR! SINGULARITY CORRESPONDING TO A COMPLEX COLLISION BETWEEN THE PLANETS " ]]
(* Bound on 1/(Q_2-sigma0 Q_1) in the real domain *)
betazeroreel=sigma0*(Lambda1+runo)^2/(Lambda2-runo)^2*mu2^2*M2/mu1^2/M1*(1+excjupmax)/(1-excsatmax);
fpremiertermereel=1/sigma0/Sqrt[1-2*betazeroreel+(betazeroreel)^2]-1/sigma0;

(* Bound on 1/(Q_2+sigma1 Q_1) in the complex domain *)
betaun=sigma1*(Lambda1+runo+r)^2/(Lambda2-runo-r)^2*mu2^2*M2/mu1^2/M1*(1+li[Lambda1,runo,r,rho,rhop,t1])/(1-li[Lambda2,runo,r,rho,rhop,t2]);
fdeuxiemeterme=1/sigma1/Sqrt[1+2*(1+xhi)*betaun+(betaun)^2]-1/sigma1;
(* Bound on 1/(Q_2+sigma1 Q_1) in the real domain *)
betaunreel=sigma1*(Lambda1+runo)^2/(Lambda2-runo)^2*mu2^2*M2/mu1^2/M1*(1+excjupmax)/(1-excsatmax);
fdeuxiemetermereel=1/sigma1/Sqrt[1+2*betaunreel+(betaunreel)^2]-1/sigma1;

(* Size of the perturbation in the complex domain *)
epsilonpert=Ggrav^2*mu1*m2*mu2^2*M2/(Lambda2-runo-r)^2/(1-li[Lambda2,runo,r,rho,rhop,t2])*(fpremierterme+fdeuxiemeterme);

(* Size of the perturbation in the real domain *)
epsilonpertreel=Ggrav^2*mu1*m2*mu2^2*M2/(Lambda2-runo)^2/(1-excsatmax)*(fpremiertermereel+fdeuxiemetermereel);

(* Check of the inequalities that analyticity widths t1 and t2 must satisfy *)

If[Not[t1<ArcCosh[Sqrt[(Lambda1-runo-r)/(1+3*rhop^2/(2(Lambda1-runo-r)))]/(rho+2rhop)]&&
t2<ArcCosh[Sqrt[(Lambda2-runo-r)/(1+3*rhop^2/(2(Lambda2-runo-r)))]/(rho+2rhop)]&&
li[Lambda1,runo,r,rho,rhop,t1]<1&&
li[Lambda2,runo,r,rho,rhop,t2]<1], Print["ERROR! PROBLEM WITH THE ANALYTICITY WIDTHS IN THE ECCENTRIC LONGITUDE"]]

(* Check whether the max of the perturbation on the subset B(xhi) is actually less than the max on the whole set  *) 
maxz=(Lambda1+runo+r)^2/(Lambda2-runo-r)^2*mu2^2*M2/mu1^2/M1*(1+li[Lambda1,runo,r,rho,rhop,t1])/(1-li[Lambda2,runo,r,rho,rhop,t2]);
sxz[x_,y_,z1_,z2_]=Norm[1/sigma0/Sqrt[1-2*(x+I*y)*sigma0*(z1+I*z2)+sigma0^2*(z1+I*z2)^2]+1/sigma1/Sqrt[1+2*(x+I*y)*sigma1*(z1+I*z2)+sigma1^2*(z1+I*z2)^2]-(m0+m1)^2/m0/m1];
valeurmax=Ggrav^2*mu1*m2*mu2^2*M2/(Lambda2-runo-r)^2/(1-li[Lambda2,runo,r,rho,rhop,t2])*MaxValue[{sxz[x,y,z1,z2],-1≤x≤1,-xhi≤y≤xhi,z1^2+z2^2≤maxz^2},{x,y,z1,z2}];
Print["Norm of the perturbation on B(xhi) ", valeurmax, " < ", epsilonpert, " ? "];
Print["Central hypothesis : ", valeurmax<epsilonpert];

(* Values *)

hamkep=hamkepP[Lambda1,Lambda2];  (* Size of the Kepler hamiltonian *)
freq1=freqkep1[Lambda1];  (* First frequency *)
freq2=freqkep2[Lambda2]; (* Second frequency *)
rapportfreq=freq1/freq2; (* Ratio of the frequencies *)
Print["Size of the Kepler hamiltonian = ",hamkep];
Print["First frequency = ",freq1];
Print["Second frequency = ",freq2];
Print["Ratio of the frequencies = ",rapportfreq];
Print["Norm of the perturbation in the complex domain = ",epsilonpert];
Print["Norm of the perturbation in the real domain ", epsilonpertreel];
Print["Complex perturbation/Kepler Hamiltonian= ", epsilonpert/hamkep];

Initial estimates and iterative lemma

Print["------------------------------------------------------"];

(* Widths consumed by the iterative lemma*)
Print["r= ", erre]; (* Actions*)
esse=lambdapmax/4;
Print["s= ", esse]; (* Angles *)
u=rhop/4;
Print["u= ",u]; (* Cartesian variables *)
beta=Sqrt[erre*esse]/u;
Print["beta= ", beta]; (* Calibration on the cartesian variables *)

Print["------------------------------------------------------"];

Print["Estimates on vector fields"];

eta0=2*epsilonpert/(erre*esse); (* Double bar norm on the non-resonant field *)
Print["eta0 = ",eta0];

gamma0=eta0/2; (* Double bar norm on the resonant field *)
Print["gamma0 = ",gamma0];

xi0=(beta^2)*eta0; (* Triple bar norm on the non-resonant field *)
Print["xi0 = ",xi0];

Gamma0=xi0/2; (* Triple bar norm on the resonant field *)
Print["Gamma0 = ",Gamma0];

(* Remainder of order 2 in the developement of the hamiltonian *)
Gcorsivo[Lambda10_,Lambda20_,Lambda1ris_,Lambda2ris_]=hamkepP[Lambda10,Lambda20]-hamkepP[Lambda1ris,Lambda2ris]-freqkep1[Lambda1ris]*(Lambda10-Lambda1ris)-freqkep2[Lambda2ris]*(Lambda20-Lambda2ris);

(* Field associated to Gcorsivo *)
derGcorsivox[Lambda10_,Lambda20_,Lambda1ris_,Lambda2ris_]=D[Gcorsivo[Lambda10,Lambda20,Lambda1ris,Lambda2ris],Lambda10];
derGcorsivoy[Lambda10_,Lambda20_,Lambda1ris_,Lambda2ris_]=D[Gcorsivo[Lambda10,Lambda20,Lambda1ris,Lambda2ris],Lambda20];

(* Borne sur le reste d'ordre 2 dans le developpement du hamiltonien *)
Deltareste[Lambda1ris_,Lambda2ris_,rruno_,larg_]=Max[Abs[derGcorsivox[Lambda1ris-rruno-larg,Lambda2ris-rruno-larg,Lambda1ris,Lambda2ris]],Abs[derGcorsivoy[Lambda1ris-rruno-larg,Lambda2ris-rruno-larg,Lambda1ris,Lambda2ris]]]/esse;
delta0=Deltareste[Lambda1,Lambda2,runo,3*erre];
Print["delta0= ", Deltareste[Lambda1,Lambda2,runo,3*erre]];

Print["------------------------------------------------------"];

(* Auxiliary functions *)
Chi[eta_,gamma_,ddelta_,xi_,grangamma_]=Max[eta+gamma+ddelta,xi+grangamma];
Theta[per_,xi_,eta_]=per/2*Max[xi,eta];

(* Coefficient on the size of the functions *)
zeta[eta_,gamma_,ddelta_,xi_,grangamma_,per_,emme_]=per*emme/2*Max[gamma+ddelta,grangamma]+emme*Theta[per,xi,eta];
Print["zeta0= ", zeta[eta0,gamma0,delta0,xi0,Gamma0,periodot,m]];

(* Coefficient on the double bar norm of the fields *)
upsilon[eta_,gamma_,ddelta_,xi_,grangamma_,per_,emme_,bbeta_]=(per*emme)^2*eta*Chi[eta,gamma,ddelta,xi,grangamma]+per*emme^2*Theta[per,xi,eta]*(2*eta+2*gamma+ddelta)+per*emme/2*Chi[eta,gamma,ddelta,xi,grangamma]+Theta[per,xi,eta]*emme*(1+gamma/eta+ddelta/eta)+(per*emme*xi/bbeta)^2*Chi[eta,gamma,ddelta,xi,grangamma]/eta+(2*per*emme^2*xi*Theta[per,xi,eta]/(beta^2))*(xi/eta+grangamma/eta);
Print["upsilon0= ",upsilon[eta0,gamma0,delta0,xi0,Gamma0,periodot,m,beta]];

(* Coefficient on the triple bar norm of the fields *)
Upsilon[eta_,gamma_,ddelta_,xi_,grangamma_,per_,emme_,bbeta_]=(bbeta*per*emme*eta)^2*Chi[eta,gamma,ddelta,xi,grangamma]/xi+(bbeta*emme)^2*per*eta*Theta[per,xi,eta]*(2*eta/xi+2*gamma/xi+ddelta/xi)+per*emme/2*Chi[eta,gamma,ddelta,xi,grangamma]+Theta[per,xi,eta]*emme*(1+grangamma/xi)+(per*emme)^2*xi*Chi[eta,gamma,ddelta,xi,grangamma]+2*per*emme^2*Theta[per,xi,eta]*(xi+grangamma);
Print["Upsilon0= ",Upsilon[eta0,gamma0,delta0,xi0,Gamma0,periodot,m,beta]];

(* Coefficient on the double bar norm of the action component of the fields (differs by a factor Theta*m*ddelta/eta with respect to upsilon0) *)
upsilonI[eta_,gamma_,ddelta_,xi_,grangamma_,per_,emme_,bbeta_]=(per*emme)^2*eta*Chi[eta,gamma,ddelta,xi,grangamma]+per*emme^2*Theta[per,xi,eta]*(2*eta+2*gamma+ddelta)+per*emme/2*Chi[eta,gamma,ddelta,xi,grangamma]+Theta[per,xi,eta]*emme*(1+gamma/eta)+(per*emme*xi/bbeta)^2*Chi[eta,gamma,ddelta,xi,grangamma]/eta+(2*per*emme^2*xi*Theta[per,xi,eta]/beta^2)*(xi/eta+grangamma/eta);
Print["upsilonI0= ",upsilonI[eta0,gamma0,delta0,xi0,Gamma0,periodot,m,beta]];

(* Necessary condition for the iterative lemma to diminish the size of the perturbation *)
(* Maximal rate (0<qq\leq 2/3) *)
qq=2/3; 
If[2*Max[upsilon[eta0,gamma0,delta0,xi0,Gamma0,periodot,m,beta],Upsilon[eta0,gamma0,delta0,xi0,Gamma0,periodot,m,beta]]>qq, Print["WARNING! CHECK THAT THE ITERATIVE LEMMA ACTUALLY DIMINISHES THE PERTURBATION!"]];

Print["------------------------------------------------------"];

(* Application of the iterative lemma *)

(* Vector containing the etas after each iteration *)
Array[eta,m+1,0];
eta[0]=eta0;

(* Vector containing the action component of the non-resonant field after each iteration *)
Array[etaI,m+1,0];
etaI[0]=eta0;

(* Vector containing the gammas after each iteration *)
Array[gamma,m+1,0];
gamma[0]=gamma0;

(* Vector containing the xis after each iteration *)
Array[xi,m+1,0];
xi[0]=xi0;

(* Vector containing the Gammas after each iteration *)
Array[GGamma,m+1,0];
GGamma[0]=Gamma0;

(* Vector containing the deltas after each iteration (smaller and smaller domains) *)
Array[delta,m+1,0];
delta[0]=delta0;

(* Vector containing the size of the non-resonant part of the perturbation after each iteration *)
Array[partnonres,m+1,0];
partnonres[0]=2*epsilonpert;

(* Vector containing the size of the resonant part of the perturbation after each iteration *)
Array[partres,m+1,0];
partres[0]=epsilonpert;

(* Iterations *)

For[j=1,j<m+1,j++,
delta[j]=Deltareste[Lambda1,Lambda2,runo,(3-2*j/m)*erre];
eta[j]=2*upsilon[eta[j-1],gamma[j-1],delta[j-1],xi[j-1],GGamma[j-1],periodot,m,beta]*eta[j-1];
etaI[j]=2*upsilonI[eta[j-1],gamma[j-1],delta[j-1],xi[j-1],GGamma[j-1],periodot,m,beta]*etaI[j-1];
gamma[j]=gamma[j-1]+eta[j]/2;
xi[j]=2*Upsilon[eta[j-1],gamma[j-1],delta[j-1],xi[j-1],GGamma[j-1],periodot,m,beta]*xi[j-1];
GGamma[j]=GGamma[j-1]+xi[j]/2;
partnonres[j]=2*zeta[eta[j-1],gamma[j-1],delta[j-1],xi[j-1],GGamma[j-1],periodot,m]*partnonres[j-1];
partres[j]=partres[j-1]+partnonres[j]/2;
];

(* Displaying the values *)
For[j=0,j<m+1,j++,
Print["eta[",j,"]=",eta[j]];
];

For[j=0,j<m+1,j++,
Print["etaI[",j,"]=",etaI[j]];
];

For[j=0,j<m+1,j++,
Print["gamma[",j,"]=",gamma[j]];
];

For[j=0,j<m+1,j++,
Print["delta[",j,"]=",delta[j]];
];

For[j=0,j<m+1,j++,
Print["xi[",j,"]=",xi[j]];
];

For[j=0,j<m+1,j++,
Print["Gamma[",j,"]=",GGamma[j]];
];

For[j=0,j<m+1,j++,
Print["partnonres[",j,"]=",partnonres[j]];
];

For[j=0,j<m+1,j++,
Print["partres[",j,"]=",partres[j]];
];

(* Size of the normal form transformation *)

(* In the actions *) 
Taillefnac= periodot/2*Sum[eta[i], {i, 0, m-1}]*erre;
Print["Size of the normal form transformation in the actions= ",Taillefnac];
Print["Size w.r.t. 2r= ", Taillefnac/(2*erre)];
If[Taillefnac>2*erre,Print["ERROR! Taillefnac>2r"]];

(* In the angles *)
Taillefnan = periodot/2*Sum[eta[i], {i, 0, m-1}]*esse;
Print["Size of the normal form transformation in the angles= ",Taillefnan];
Print["Size w.r.t. 2s= ", Taillefnan/(2*esse)];
If[Taillefnan>(2*esse),Print["ERROR! Taillefnan>2s"]];

(* In the cartesian variables *)
Taillefnvd = periodot/2*Sum[xi[i], {i, 0, m-1}]*u;
Print["Size of the normal form transformation in the cartesian variables= ",Taillefnvd];
Print["Size w.r.t. 2u= ", Taillefnvd/(2*u)];
If[Taillefnvd>(2*u),Print["ERROR! Taillefnvd>2*u"]];
If[Taillefnvd>alphafnxy*u,Print["ERROR! The internal radius of analyticity in the cartesian variables is too big!"]];

Convexity and stability

(* Lower bounds on the eigenvalues of the Hessian of the unperturbed hamiltonian *)
kappagreco=Min[Abs[derhamkepPxx[Lambda1+rf-Taillefnac,0]],Abs[derhamkepPyy[0,Lambda2+rf-Taillefnac]]];
Print["kappagreco= ", kappagreco];

(* Upper bound on the Hessian matrix norm *)
K=Abs[maxderhamkepPxx[Lambda1,ri,Taillefnac]+maxderhamkepPyy[Lambda2,ri,Taillefnac]];
Print["K= ",K];
If[K<kappagreco, Print["WARNING! PROBLEM WITH THE HESSIAN OF THE KEPLER HAMILTONIAN!"]];

(* Radius of initial conditions in the normalized variables *)
Rtilde=ri+Taillefnac;
(* Hypotheses in the stability theorem *)
If[Rtilde>rf-Taillefnac, Print["ERROR! The radius of initial conditions in the normalized variables is too big!"]];

(* Maximal variation in the action variables *)
relvarI=rf/Max[Lambda1,Lambda2]; (* Maximal variation in the action variables*)

(* Constants appearing in the stability theorem *)

C1=kappagreco/2*( (rf-Taillefnac-(K/kappagreco+1)*Rtilde)^2 - (K/kappagreco*Rtilde)^2 ) - 2*(partres[m]+partnonres[m]);

C2=erre*( Abs[freqkep1[Lambda1]+freqkep2[Lambda2]] )*etaI[m];

(* Analiticity widths *)
Print["rho= ", runo];
Print["r/Max[Lambda1,Lambda2]= ", erre/Max[Lambda1,Lambda2]];
Print["m= ", m];
Print["s= ", esse];
Print["u= ", u];
Print["beta= ", beta];
Print["xi= ", rho];
Print["Rf/Max[Lambda1,Lambda2]= ", relvarI];
Print["excjupmax= ", excjupmax];
Print["excsatmax= ", excsatmax];

(* Time of stability *)

tstab=C1/C2;
Print["Time of stability (y)= ", tstab/(3.1536*10^7)];