Sharp Nekhoroshev estimates for the three-body problem around periodic orbits
Author: Santiago Barbieri, PhD student

Astronomical data (Nasa database)

(* Definition Astronomical Unit (m) *) 
AU=149597870700;

(* Newton's constant (N*m^2/kg^2) *)
Ggrav=6.67259*10^(-11);

(* Jupiter's orbit radius (the orbit is supposed to be circular) *)
aj=5.204*AU;

(* Jupiter's real eccentricity *)
ej=4.89*10^(-2);

(* Tuning on the mass of Jupiter *)
k=2.5;
Print["We are in the m_jupiter=10^(",-k-3,") m_sun case"];

(* Masses (in kg) *)
sunmass=1.9885*10^30; (* Mass of the Sun *)
m0=sunmass*10^(00);     (* Central Mass *)
m1=1.89819*10^27*10^(-k); (* Mass of Jupiter *)
M=m0+m1;
mu=1/(Ggrav*m0)^(2/3);

(* Perturbative parameter *)
epsilon=m1/((m0)^(2/3))*1/(m0+m1)^(1/3);
Print["epsilon= ", epsilon];

(* Angular velocity of Jupiter *)
omegac=Sqrt[Ggrav*M/aj^3];

(*Period of Jupiter *)
Tj=2*Pi/omegac;
Print["Test. Period of Jupiter (y): ", Tj/(3.1536*10^7)];

Hamiltonian structure for the restricted three-body problem

(****************************** Delaunay variables **********************************************)
L0[aa_]=mu*Sqrt[Ggrav*m0*aa];
G0[aa_,ee_]=L0[aa]*Sqrt[1-ee^2];
a0[LL_]=((LL/mu)^2)/(Ggrav*m0); (* semi-major axis *)
e0[LL_,GG_]=Sqrt[1-(GG/LL)^2]; (* eccentricity *)
Lj=L0[aj]; (* First Delaunay variable for Jupiter *)
Gj=G0[aj,ej];
Print["Test. Jupiter semi-major axis: ", a0[Lj]/(AU), " AU"];
Print["Test. Jupiter eccentricity: ", e0[Lj,Gj]];

(******************************** Restricted three-body problem hamiltonian and its associated fields **************************************)

(* Keplerian part *)
H0[L_,G_]=-1/(2*L^2)-omegac*G;
(* Keplerian frequencies *)
freq1[L_]=D[H0[L,G],L];
freq2=-omegac;

(* Semi-major axes ratio *)
ratioa[LL1_,LL2_]=a0[LL1]/a0[LL2]; 

(* Perturbation (from the article "KAM stability and celestial mechanics" by Celletti and Chierchia, AMS, 2003) *)
Pert[L_,G_,l_,g_]=1+1/4*(ratioa[L,Lj])^2+9/64*(ratioa[L,Lj])^4+3/8*(ratioa[L,Lj])^2*(e0[L,G])^2-(1/2+9/16*(ratioa[L,Lj])^2)*(ratioa[L,Lj])^2*e0[L,G]*Cos[l]+(3/8*(ratioa[L,Lj])^3+15/64*(ratioa[L,Lj])^5)*Cos[l+g]-(9/4+5/4*(ratioa[L,Lj])^2)*(ratioa[L,Lj])^2*e0[L,G]*Cos[l+2g]+(3/4*(ratioa[L,Lj])^2+5/16*(ratioa[L,Lj])^4)*Cos[2l+2g]+3/4*(ratioa[L,Lj])^2*e0[L,G]*Cos[3l+2g]+(5/8*(ratioa[L,Lj])^3+35/128*(ratioa[L,Lj])^5)*Cos[3l+3g]+35/64*(ratioa[L,Lj])^4*Cos[4l+4g]+63/128*(ratioa[L,Lj])^5*Cos[5l+5g];
H1[L_,G_,l_,g_]=-(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)*Pert[L,G,l,g];

(* Total Hamiltonian *)
H[L_,G_,l_,g_]=H0[L,G]+H1[L,G,l,g];

(* Resonant part of the perturbation (the 1/3 resonance that we use to make our Nekhoroshev calculations is in the discarded harmonics, so that the only resonant part of the perturbation is the polynomial one) *)
RH[L_,G_]=Simplify[-(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)*(1+1/4*(ratioa[L,Lj])^2+9/64*(ratioa[L,Lj])^4+3/8*(ratioa[L,Lj])^2*(e0[L,G])^2)];

(* Non-resonant part of the perturbation (composed by the trigonometric part of the Celletti-Chierchia perturbation) *)
NRH[L_,G_,l_,g_]=Simplify[-(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)*(-(1/2+9/16*(ratioa[L,Lj])^2)*(ratioa[L,Lj])^2*e0[L,G]*Cos[l]+(3/8*(ratioa[L,Lj])^3+15/64*(ratioa[L,Lj])^5)*Cos[l+g]-(9/4+5/4*(ratioa[L,Lj])^2)*(ratioa[L,Lj])^2*e0[L,G]*Cos[l+2g]+(3/4*(ratioa[L,Lj])^2+5/16*(ratioa[L,Lj])^4)*Cos[2l+2g]+3/4*(ratioa[L,Lj])^2*e0[L,G]*Cos[3l+2g]+(5/8*(ratioa[L,Lj])^3+35/128*(ratioa[L,Lj])^5)*Cos[3l+3g]+35/64*(ratioa[L,Lj])^4*Cos[4l+4g]+63/128*(ratioa[L,Lj])^5*Cos[5l+5g])];

(* Remainder of order two in the developement of the integrable hamiltonian (by Taylor's formula) *)
Rem2[LL_]=H0[LL,GG]-H0[Linit,GGris]-freq1[Linit]*(LL-Linit)-freq2*(GG-GGris);

(* Resonant hamiltonian field *)
XRHL=0;
XRHG=0;
XRHl[L_,G_]=Simplify[D[RH[L,G],L]];
XRHg[L_,G_]=Simplify[D[RH[L,G],G]];

(* Non-resonant hamiltonian field *)
XNRHL[L_,G_,l_,g_]=Simplify[-D[NRH[L,G,l,g],l]];
XNRHG[L_,G_,l_,g_]=Simplify[-D[NRH[L,G,l,g],g]];
XNRHl[L_,G_,l_,g_]=Simplify[D[NRH[L,G,l,g],L]];
XNRHg[L_,G_,l_,g_]=Simplify[D[NRH[L,G,l,g],G]];

(* Remainder of order 2 hamiltonian field *)
XRem2L=0;
XRem2G=0;
XRem2l[LL_]=Simplify[D[Rem2[LL],LL]];
XRem2g=0;

Resonance and initial eccentricity

(* Mean motion resonance with Jupiter *)
pp=1;
qq=3;
freqratio=pp/qq;

(* Choice of the initial eccentricity for the asteroid (arbitrary) *)

einit=0;

(* Associated Delaunay variables *)
Linit=(1/omegac*freqratio)^(1/3);
Ginit=Linit*Sqrt[1-einit^2];

(* Initial, maximum and minimum values of the semi-major axes on the considered domain *)
ainit=a0[Linit];

Print["One has ","a_init=",ainit/AU," AU for the ", pp/qq," resonance"];

(* Resonant frequencies *)
omegaLinit=Abs[freq1[Linit]];
omegaGinit=Abs[freq2];
Print["Test. omegaG_init/omegaL_init= ", Rationalize[omegaGinit/omegaLinit]];

(* Period of the resonant torus *)

T=qq*2*Pi/omegaLinit;
Print["Test. T/Tj= ", T/Tj];

Analyticity widths

(* Radius of the disk in the L action *)
rL=2*10^(-3)*Linit;
(* Width consumed normal form lemma in the L action *)
erreL=rL/3;
(* Fraction of the analiticity width consumed by the initial numerical normal form *)
alphaNFL=2.4*10^(-1);
erreNFL=alphaNFL*erreL;

(* Analiticity width in the angle l conjugated to action L *)
sl=1.05;
(* Width consumed by the normal form lemma in the l angle *)
essel=sl/3;
(* Fraction of essel consumed by the initial numerical normal form *)
alphaNFl=3*10^(-5);
esseNFl=alphaNFl*essel;

(* It is necessary to consider a polydisk in the G variables, not just a disk! *)
(* Basis for the polydisk in the G action *)
rGbasis=8.5*10^(-3)*Ginit;
(* Radius of the polydisk in the G action *)
rG=1.8*10^(-3)*Ginit;
(* Width consumed by the normal form lemma in the G action *)
erreG=rG/3;
(* Fraction of erreG consumed by the initial numerical normal form *)
alphaNFG=3*10^(-2);
erreNFG=alphaNFG*erreG;

(* Angle g conjugated to action G *)
sg=8.1*10^(-4);
(* Width consumed by the normal form lemma in the g angle *)
esseg=sg/3;
(* Fraction of essel consumed by the initial numerical normal form *)
alphaNFg=7.6*10^(-1);
esseNFg=alphaNFg*esseg;

(* Number of steps in the iterative lemma *)

z=3;

Validity of the Celletti-Chierchia model and extremum for the eccentricity

(************************************ Extrema for the variables in the considered domain **************************************)

Lmax=Linit+2*rL;
Lmin=Linit-2*rL;
Gmax=Ginit+rGbasis+2*rG;
Gmin=Ginit-rGbasis-2*rG;

amax=a0[Lmax];
amin=a0[Lmin];

Print["For the resonance ", pp/qq, " one gets ","a_max=",amax/AU," AU"];
Print["For the resonance ", pp/qq, " one gets ","a_min=",amin/AU," AU"];

(* Initial conditions *)
(* Magnitude of the biggest discarded value in Celletti and Chierchia (AMS, 2007), p.79 *)
minharm=1.99*10^(-2);

(* a^6 is discarded by Celletti and Chierchia *)
If[(ainit/aj)^6>minharm,Print["!!!! WARNING INITIAL CONDITIONS !!!!!!a^6 not negligible!!!!"]];

(* e^3 is discarded by Celletti and Chierchia *)
If[(einit)^3>minharm,Print["!!!! WARNING INITIAL CONDITIONS !!!!!!e^3 not negligible!!!!"]];

(* a^2*e^2 is discarded by Celletti and Chierchia *)
If[(ainit/aj)^2*einit^2>minharm,Print["!!!! WARNING INITIAL CONDITIONS !!!!!!a^2*e^2 not negligible!!!!"]];

(* a^3*e is discarded by Celletti and Chierchia *)
If[(ainit/aj)^3*einit>minharm,Print["!!!! WARNING INITIAL CONDITIONS !!!!!!a^3*e not negligible!!!!"]];

(* Final conditions *)
(* The final configuration of the system must be in a region where the Celletti-Chierchia model still holds *)

(* a^6 is discarded by Celletti and Chierchia *)
If[(amax/aj)^6>minharm,Print["!!!! WARNING FINAL CONDITIONS !!!!!!a^6 not negligible!!!!"]];

(* The other three discarded terms impose an upper bound on the final eccentricity *)
ealpha=(minharm)^(1/3);
ebeta=minharm*(amin/aj)^(-3);
egamma=Sqrt[minharm]*(amin/aj)^(-1);

(* This is the bound on the final eccentricity imposed by the model *)
emodel=Min[ealpha,ebeta,egamma];
Print["The bound on the final eccentricity imposed by the model is: ", emodel];

(* The displacement between initial and final conditions also implies an upper bound on the eccentricities by the mean of the conservation of energy *)

(* Choice of the random angles to calculate the difference in value between H1 calculated in ainit and H1 calculated in amax/amin *)
H1points=Table[ RandomReal[{-Pi,Pi},{4}] , 10000];
(* Difference between H1 calculated in ainit and in amax *)
DeltaH1=Max[Max@Abs@Map[ Simplify[H1[Linit+(#[[3]]/Pi)*(Lmax-Linit),Ginit+(#[[4]]/Pi)*(Gmax-Ginit),#[[1]],#[[2]]]- H1[Linit,Ginit,#[[1]],#[[2]]]]&, H1points],Max@Abs@Map[ H1[Linit,Ginit,#[[1]],#[[2]]]- H1[Lmax,Lmax*Sqrt[1-emodel^2],#[[1]],#[[2]]]&, H1points],Max@Abs@Map[ H1[Linit,Ginit,#[[1]],#[[2]]]- H1[Lmin,Lmin*Sqrt[1-emodel^2],#[[1]],#[[2]]]&, H1points]];

(* Variation of the L actions *)
K1=Max[Abs[1/(2*Lmin^2)-1/(2*Linit^2)],Abs[1/(2*Lmax^2)-1/(2*Linit^2)]];
(* Maximal variation in the G actions *)
maxvarG=1/omegac*(K1+DeltaH1);
Print["Ginit+maxvarG= ", Ginit+maxvarG];
Print["Gmax= ", Gmax];

(* Check that the maximal variation in the G actions in well contained in the domain of analyticity *)
If[Ginit+maxvarG>Gmax, Print["!!!!WARNING!!!! The G variables vary too much and go outside the considered domain."]];
yminus=1/amax*(ainit*(1-einit^2)+(K1+DeltaH1)^2/((omegac*mu)^2*Ggrav*m0))-2*(K1+DeltaH1)/(amax*omegac*mu)*Sqrt[ainit*(1-einit^2)/(Ggrav*m0)];

(* Constraint on the eccentricity *)
emax=Sqrt[1-yminus];
Print["Bound on the eccentricity imposed by the conservation of energy: ",emax];
If[emax>emodel,Print["!!!!WARNING!!!! The maximal eccentricity is greater than the one allowed by the model"]];

(************************* Check that the G/g variables stay in the right domain ***************************)

(* Values of the different terms in the Celletti and Chierchia model (applied to Victoria and Iris) *)
Print["Celletti-Chierchia"];
Print["Victoria"];
av=0.449;
ev=0.22;
Print["Retained terms"];
Print["av^5= ", av^5];
Print["av^2*ev= ", av^2*ev];
Print["Discarded terms"];
Print["av^6= ", av^6];
Print["av^2*e^2= ", av^2*ev^2];
Print["ev^3= ", ev^3];
Print["av^3*ev= ", av^3*ev];
Print["-------------------------"];
Print["Me"];
Print["Retained terms"];
Print["amax^5= ", (amax/aj)^5];
Print["amax^2*emax= ", (amax/aj)^2*emax];
Print["Discarded terms"];
Print["amax^6= ", (amax/aj)^6];
Print["amax^2*emax^2= ", (amax/aj)^2*emax^2];
Print["emax^3= ", emax^3];
Print["amax^3*emax= ", (amax/aj)^3*emax];
Print["-------------------------"];

Numerical normal form 1st step (following Gallavotti, 1983)

(* Before applying Nekhoroshev theorem we reduce the size of the perturbation in the style of Celletti and Ferrara (1996) *)

(* Generator of the normal form trasformation in a mixed set of coordinates *)
omega1=freq1[Linit];
omega2=-omegac;
Phione[Lp_,Gp_,l_,g_]=-(1/2+9/16*(ratioa[Lp,Lj])^2)*(ratioa[Lp,Lj])^2*e0[Lp,Gp]*Sin[l]/omega1+(3/8*(ratioa[Lp,Lj])^3+15/64*(ratioa[Lp,Lj])^5)*Sin[l+g]/(omega1+omega2)-(9/4+5/4*(ratioa[Lp,Lj])^2)*(ratioa[Lp,Lj])^2*e0[Lp,Gp]*Sin[l+2g]/(omega1+2*omega2)+(3/4*(ratioa[Lp,Lj])^2+5/16*(ratioa[Lp,Lj])^4)*Sin[2l+2g]/(2*omega1+2*omega2)+3/4*(ratioa[Lp,Lj])^2*e0[Lp,Gp]*Sin[3l+2g]/(3*omega1+2*omega2)+(5/8*(ratioa[Lp,Lj])^3+35/128*(ratioa[Lp,Lj])^5)*Sin[3l+3g]/(3*omega1+3*omega2)+35/64*(ratioa[Lp,Lj])^4*Sin[4l+4g]/(4*omega1+4*omega2)+63/128*(ratioa[Lp,Lj])^5*Sin[5l+5g]/(5*omega1+5*omega2);

(* New perturbation in a mixed set of coordinates *)
H2[Lp_,Gp_,l_,g_]=Simplify[1/2*D[D[H0[Lp,Gp],Lp],Lp]*(D[Phione[Lp,Gp,l,g],l])^2+D[Pert[Lp,Gp,l,g],Lp]*D[Phione[Lp,Gp,l,g],l]+D[Pert[Lp,Gp,l,g],Gp]*D[Phione[Lp,Gp,l,g],g]];

(* Rescaling with respect to the parameters *)
H2rescaled[Lp_,Gp_,l_,g_]=(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^2*H2[L,G,l,g];

Numerical normal form 2nd step (following Celletti - Ferrara, 1996)

(* Extracting the Fourier coefficients of H2 *)

H2exp[L_,G_,l_,g_]=Simplify[TrigToExp[H2[L,G,l,g]]];

Print["I am calculating H2exp00"];
H2exp00[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{0,0}];
Print["I am calculating H2exp10"];
H2exp10[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{1,0}];
Print["I am calculating H2exp20"];
H2exp20[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{2,0}];
Print["I am calculating H2exp11"];
H2exp11[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{1,1}];
Print["I am calculating H2exp12"];
H2exp12[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{1,2}];
Print["I am calculating H2exp21"];
H2exp21[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{2,1}];
Print["I am calculating H2exp22"];
H2exp22[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{2,2}];
Print["I am calculating H2exp23"];
H2exp23[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{2,3}];
Print["I am calculating H2exp24"];
H2exp24[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{2,4}];
Print["I am calculating H2exp32"];
H2exp32[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{3,2}];
Print["I am calculating H2exp33"];
H2exp33[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{3,3}];
Print["I am calculating H2exp34"];
H2exp34[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{3,4}];
Print["I am calculating H2exp42"];
H2exp42[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{4,2}];
Print["I am calculating H2exp43"];
H2exp43[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{4,3}];
Print["I am calculating H2exp44"];
H2exp44[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{4,4}];
Print["I am calculating H2exp45"];
H2exp45[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{4,5}];
Print["I am calculating H2exp54"];
H2exp54[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{5,4}];
Print["I am calculating H2exp55"];
H2exp55[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{5,5}];
Print["I am calculating H2exp64"];
H2exp64[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{6,4}];
Print["I am calculating H2exp01"];
H2exp01[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{0,1}];
Print["I am calculating H2exp02"];
H2exp02[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{0,2}];
Print["I am calculating H2exp56"];
H2exp56[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{5,6}];
Print["I am calculating H2exp65"];
H2exp65[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{6,5}];
Print["I am calculating H2exp66"];
H2exp66[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{6,6}];
Print["I am calculating H2exp67"];
H2exp67[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{6,7}];
Print["I am calculating H2exp77"];
H2exp77[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{7,7}];
Print["I am calculating H2exp76"];
H2exp76[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{7,6}];
Print["I am calculating H2exp87"];
H2exp87[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{8,7}];
Print["I am calculating H2exp88"];
H2exp88[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{8,8}];
Print["I am calculating H2exp99"];
H2exp99[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{9,9}];
Print["I am calculating H2exp1010"];
H2exp1010[L_,G_]=FourierCoefficient[H2exp[L,G,l,g],{l,g},{10,10}];

(* Since the resonant part of H1 only depends on the actions, we can rescale the frequencies up to a factor of order epsilon and the non-resonant remainder decreases as epsilon^4 after the second step of the normal form *)

(* Frequencies rescaled *)

freq3[L_,G_]=freq1[L]+D[RH[L,G],L];
freq4[L_,G_]=freq2+D[RH[L,G],G];

omega3=freq3[Linit,Ginit];
omega4=freq4[Linit,Ginit];


(* Generator of the normal form trasformation in a mixed set of coordinates *)
Phitwo[Lp_,Gp_,l_,g_]=-2*(
H2exp10[Lp,Gp]/omega3*Sin[l]+
H2exp20[Lp,Gp]/omega3*Sin[2*l]+
H2exp11[Lp,Gp]/(omega3+omega4)*Sin[l+g]+
H2exp12[Lp,Gp]/(omega3+2*omega4)*Sin[l+2*g]+
H2exp22[Lp,Gp]/(2*omega3+2*omega4)*Sin[2*l+2*g]+
H2exp32[Lp,Gp]/(3*omega3+2*omega4)*Sin[3*l+2*g]+
H2exp33[Lp,Gp]/(3*omega3+3*omega4)*Sin[3*l+3*g]+
H2exp44[Lp,Gp]/(4*omega3+4*omega4)*Sin[4*l+4*g]+
H2exp55[Lp,Gp]/(5*omega3+5*omega4)*Sin[5*l+5*g]+
H2exp24[Lp,Gp]/(2*omega3+4*omega4)*Sin[2*l+4*g]+
H2exp64[Lp,Gp]/(6*omega3+4*omega4)*Sin[6*l+4*g]+
H2exp66[Lp,Gp]/(6*omega3+6*omega4)*Sin[6*l+6*g]+
H2exp88[Lp,Gp]/(8*omega3+8*omega4)*Sin[8*l+8*g]+
H2exp1010[Lp,Gp]/(10*omega3+10*omega4)*Sin[10*l+10*g]+
H2exp21[Lp,Gp]/(2*omega3+omega4)*Sin[2*l+g]+
H2exp42[Lp,Gp]/(4*omega3+2*omega4)*Sin[4*l+2*g]+
H2exp43[Lp,Gp]/(4*omega3+3*omega4)*Sin[4*l+3*g]+
H2exp54[Lp,Gp]/(5*omega3+4*omega4)*Sin[5*l+4*g]+
H2exp65[Lp,Gp]/(6*omega3+5*omega4)*Sin[6*l+5*g]+
H2exp23[Lp,Gp]/(2*omega3+3*omega4)*Sin[2*l+3*g]+
H2exp34[Lp,Gp]/(3*omega3+4*omega4)*Sin[3*l+4*g]+
H2exp45[Lp,Gp]/(4*omega3+5*omega4)*Sin[4*l+5*g]+
H2exp56[Lp,Gp]/(5*omega3+6*omega4)*Sin[5*l+6*g]+
H2exp67[Lp,Gp]/(6*omega3+7*omega4)*Sin[6*l+7*g]+
H2exp77[Lp,Gp]/(7*omega3+7*omega4)*Sin[7*l+7*g]+
H2exp76[Lp,Gp]/(7*omega3+6*omega4)*Sin[7*l+6*g]+
H2exp87[Lp,Gp]/(8*omega3+7*omega4)*Sin[8*l+7*g]+
H2exp99[Lp,Gp]/(9*omega3+9*omega4)*Sin[9*l+9*g]+
H2exp01[Lp,Gp]/(omega4)*Sin[g]+
H2exp02[Lp,Gp]/(2*omega4)*Sin[2*g]
);

(* New epsilon^4 perturbation in a mixed set of coordinates *)
H3aux[Lp_,Gp_,l_,g_]=1/2*(D[D[H0[Lp,Gp],Lp],Lp]+D[D[RH[Lp,Gp],Lp],Lp])*(D[Phitwo[Lp,Gp,l,g],l])^2+1/2*(D[D[RH[Lp,Gp],Gp],Gp])*(D[Phitwo[Lp,Gp,l,g],g])^2+1/2*(D[D[RH[Lp,Gp],Lp],Gp])*(D[Phitwo[Lp,Gp,l,g],g])*(D[Phitwo[Lp,Gp,l,g],l])+D[H2[Lp,Gp,l,g],Lp]*D[Phitwo[Lp,Gp,l,g],l]+D[H2[Lp,Gp,l,g],Gp]*D[Phitwo[Lp,Gp,l,g],g];

(* Rescaling the new perturbation *)
H3[Lp_,Gp_,l_,g_]=(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^4*H3aux[Lp,Gp,l,g];

Initial estimates

(* Vector field of H3 *)
XH3L[Lp_,Gp_,l_,g_]=D[H3[Lp,Gp,l,g],l];
XH3G[Lp_,Gp_,l_,g_]=D[H3[Lp,Gp,l,g],g];
XH3l[Lp_,Gp_,l_,g_]=D[H3[Lp,Gp,l,g],Lp];
XH3g[Lp_,Gp_,l_,g_]=D[H3[Lp,Gp,l,g],Gp];

(* Resonant vector field of the previous remainder *)
XRH2l[L_,G_]=D[H2exp00[L,G],L];
XRH2g[L_,G_]=D[H2exp00[L,G],G];

(* Size of the total normal form *)
funcsizenfL[Lp_,Gp_,l_,g_]=(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^2*Abs[D[Phione[Lp,Gp,l,g],l]]+(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^4*Abs[D[Phitwo[Lp,Gp,l,g],l]];
funcsizenfG[Lp_,Gp_,l_,g_]=(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^2*Abs[D[Phione[Lp,Gp,l,g],g]]+(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^4*Abs[D[Phitwo[Lp,Gp,l,g],g]];
funcsizenfl[Lp_,Gp_,l_,g_]=(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^2*Abs[D[Phione[Lp,Gp,l,g],Lp]]+(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^4*Abs[D[Phitwo[Lp,Gp,l,g],Lp]];
funcsizenfg[Lp_,Gp_,l_,g_]=(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^2*Abs[D[Phione[Lp,Gp,l,g],Gp]]+(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^4*Abs[D[Phitwo[Lp,Gp,l,g],Gp]];

(* Choice of the random boundary points *)
boundarypoints=Table[ RandomReal[{-Pi,Pi},{7}] , 100];

(* Size of the normal form in the actions L*)
sizenfL=Max[
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfL[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints]];

(* Size of the normal form in the actions G*)
sizenfG=Max[
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfG[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints]];


(* Size of the normal form in the angles l*)
sizenfl=Max[
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints]];

(* Size of the normal form in the angles g*)
sizenfg=Max[
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  funcsizenfg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints]];

(* Check that the size of the normal form is not too big *)
If[sizenfL>erreNFL,Print["ERROR! THE SIZE OF THE NORMAL FORM IS TOO BIG IN L!"]];
If[sizenfG>erreNFG,Print["ERROR! THE SIZE OF THE NORMAL FORM IS TOO BIG IN G!"]];
If[sizenfl>esseNFl,Print["ERROR! THE SIZE OF THE NORMAL FORM IS TOO BIG IN l!"]];
If[sizenfg>esseNFg,Print["ERROR! THE SIZE OF THE NORMAL FORM IS TOO BIG IN g!"]];

(* Size on the non-resonant part of H3 *)

NormNRH3=2*Max[
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  H3[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints]];

(* Size of the resonant part of H3 - it is given by the resonant part RH plus the resonant part of the first remainder plus the size of the 
second one *)

NormRH3=Max[Max@Abs@Map[  RH[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  RH[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  RH[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2)))]&, boundarypoints],
Max@Abs@Map[  RH[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2)))]&, boundarypoints]]+(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^2*Max[Max@Abs@Map[  H2exp00[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  H2exp00[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  H2exp00[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2)))]&, boundarypoints],
Max@Abs@Map[  H2exp00[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2)))]&, boundarypoints]]+NormNRH3/2;

(* Size of the non-resonant vector fields *)

NormXNRH3L=2*Max[
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3L[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints]];


NormXNRH3G=2*Max[
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3G[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints]];

NormXNRH3l=2*Max[
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints]];

NormXNRH3g=2*Max[
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]+sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]+sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I,#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints],
Max@Abs@Map[  XH3g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2))),#[[3]]-sl*I,#[[4]]-sg*I]&, boundarypoints]];


(* Size of the resonant vector fields *)

NormXRH3L=0;

NormXRH3G=0;

NormXRH3l=Max[
Max@Abs@Map[  XRHl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  XRHl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  XRHl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2)))]&, boundarypoints],
Max@Abs@Map[  XRHl[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2)))]&, boundarypoints]]+(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^2*Max[
Max@Abs@Map[  XRH2l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  XRH2l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  XRH2l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2)))]&, boundarypoints],
Max@Abs@Map[  XRH2l[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2)))]&, boundarypoints]]+NormXNRH3l/2;

NormXRH3g=Max[
Max@Abs@Map[  XRHg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  XRHg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  XRHg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2)))]&, boundarypoints],
Max@Abs@Map[  XRHg[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2)))]&, boundarypoints]]+(epsilon*(Ggrav*(m0+m1))^(1/3)/aj)^2*Max[
Max@Abs@Map[  XRH2g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+#[[5]]/Pi*rGbasis+(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  XRH2g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-#[[5]]/Pi*rGbasis-(Gmax-Ginit-rGbasis-erreNFG)I]&, boundarypoints],
Max@Abs@Map[  XRH2g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit+rGbasis+(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[6]]/2)))]&, boundarypoints],
Max@Abs@Map[  XRH2g[Linit+(Lmax-Linit-erreNFL) (E^(I #[[1]])),Ginit-rGbasis-(Gmax-Ginit-rGbasis-erreNFG)*(E^(I (#[[7]]/2)))]&, boundarypoints]]+NormXNRH3g/2;




Print["---------------------------------------------------------"];
Print["Initial bounds on functions and vector fields"];
Print["---------------------------------------------------------"];
(* Bound on the field associated to the remainder of order 2 in the developement of the unperturbed hamiltonian *)
delta0=Max[Abs[XRem2l[Linit-rL]]/essel,Abs[XRem2l[Linit+rL]]/essel];
Print["delta0= ",delta0];
eta0L=NormXNRH3L/erreL;
eta0G=NormXNRH3G/erreG;
eta0l=NormXNRH3l/essel;
eta0g=NormXNRH3g/esseg;
Print["eta0L=",eta0L];
Print["eta0G=",eta0G];
Print["eta0l=",eta0l];
Print["eta0g=",eta0g];
gamma0L=NormXRH3L/erreL;
gamma0G=NormXRH3G/erreG;
gamma0l=NormXRH3l/essel;
gamma0g=NormXRH3g/esseg;
Print["gamma0L=",gamma0L];
Print["gamma0G=",gamma0G];
Print["gamma0l=",gamma0l];
Print["gamma0g=",gamma0g];

Decreasing terms

Print["---------------------------------------------------------"];
Print["Decreasing terms"];

chifunc[etaElle_,etaelle_, etaGi_, etagi_ ,gammaElle_,gammaelle_, gammaGi_, gammagi_,ddelta_]=Max[etaElle+gammaElle,etaGi+gammaGi,etaelle+gammaelle+ddelta,etagi+gammagi];
chi0=chifunc[eta0L,eta0l, eta0G, eta0g ,gamma0L,gamma0l, gamma0G, gamma0g,delta0];
Print["chi0= ", chi0];

thetafunc[etaElle_,etaelle_, etaGi_, etagi_ , per_]=per/2*Max[etaElle,etaelle,etaGi,etagi];
theta0=thetafunc[eta0L,eta0l, eta0G, eta0g , T];
Print["theta0= ", theta0];

upsilonL[etaElle_,etaelle_, etaGi_, etagi_ ,gammaElle_,gammaelle_, gammaGi_, gammagi_,ddelta_, per_, emme_,cchi_,ttheta_]=((per*emme)^2*etaelle*cchi)/2+per*emme*cchi/2+(1+gammaElle/etaElle)*emme*ttheta+sg/sl*rG/rL*per*emme^2/(2*etaElle)*(per*etagi*etaGi*cchi+(etagi*(etaGi+gammaGi)+etaGi*(etagi+gammagi))*ttheta)+per*emme^2/2*(etaelle*(1+gammaElle/etaElle)+etaelle+gammaelle+ddelta)*ttheta;
upsilon0L=upsilonL[eta0L,eta0l,eta0G,eta0g,gamma0L,gamma0l,gamma0G,gamma0g,delta0,T,z,chi0,theta0];
Print["upsilon0L= ",upsilon0L];

upsilonG[etaElle_,etaelle_, etaGi_, etagi_ ,gammaElle_,gammaelle_, gammaGi_, gammagi_,ddelta_, per_, emme_,cchi_,ttheta_]=((per*emme)^2*etagi*cchi)/2+per*emme*cchi/2+(1+gammaGi/etaGi)*emme*ttheta+sl/sg*rL/rG*per*emme^2/(2*etaGi)*(per*etaelle*etaElle*cchi+(etaelle*(etaElle+gammaElle)+etaElle*(etaelle+gammaelle+ddelta))*ttheta)+per*emme^2/2*(etagi*(1+gammaGi/etaGi)+etagi+gammagi)*ttheta;
upsilon0G=upsilonG[eta0L,eta0l,eta0G,eta0g,gamma0L,gamma0l,gamma0G,gamma0g,delta0,T,z,chi0,theta0];
Print["upsilon0G= ",upsilon0G];

upsilonl[etaElle_,etaelle_, etaGi_, etagi_ ,gammaElle_,gammaelle_, gammaGi_, gammagi_,ddelta_, per_, emme_,cchi_,ttheta_]=((per*emme)^2*etaElle*cchi)/2+per*emme*cchi/2+(1+gammaelle/etaelle+ddelta/etaelle)*emme*ttheta+sg/sl*rG/rL*per*emme^2/(2*etaelle)*(per*etagi*etaGi*cchi+(etagi*(etaGi+gammaGi)+etaGi*(etagi+gammagi))*ttheta)+per*emme^2/2*(etaElle*(1+gammaelle/etaelle+ddelta/etaelle)+etaElle+gammaElle)*ttheta;
upsilon0l=upsilonl[eta0L,eta0l,eta0G,eta0g,gamma0L,gamma0l,gamma0G,gamma0g,delta0,T,z,chi0,theta0];
Print["upsilon0l= ",upsilon0l];

upsilong[etaElle_,etaelle_, etaGi_, etagi_ ,gammaElle_,gammaelle_, gammaGi_, gammagi_,ddelta_, per_, emme_,cchi_,ttheta_]=((per*emme)^2*etaGi*cchi)/2+per*emme*cchi/2+(1+gammagi/etagi)*emme*ttheta+sl/sg*rL/rG*per*emme^2/(2*etagi)*(per*etaelle*etaElle*cchi+(etaelle*(etaElle+gammaElle)+etaElle*(etaelle+gammaelle+ddelta))*ttheta)+per*emme^2/2*(etaGi*(1+gammagi/etagi)+etaGi+gammaGi)*ttheta;
upsilon0g=upsilong[eta0L,eta0l,eta0G,eta0g,gamma0L,gamma0l,gamma0G,gamma0g,delta0,T,z,chi0,theta0];
Print["upsilon0g= ",upsilon0g];

zeta[per_,emme_,cchi_]=(per*emme/2)*cchi;
zeta0=zeta[T,z,chi0];
Print["zeta0= ",zeta0];

(***************************** Iterative lemma ****************************************)

(* Arrays containing the sizes of the vector fields and functions after successive applications of the iterative lemma*)

Array[etaL,z,0];
etaL[0]=eta0L;

Array[etaG,z,0];
etaG[0]=eta0G;

Array[etal,z,0];
etal[0]=eta0l;

Array[etag,z,0];
etag[0]=eta0g;

Array[gammaL,z,0];
gammaL[0]=gamma0L;

Array[gammaG,z,0];
gammaG[0]=gamma0G;

Array[gammal,z,0];
gammal[0]=gamma0l;

Array[gammag,z,0];
gammag[0]=gamma0g;

Array[delta,z,0];
delta[0]=delta0;

Array[zeta,z,0];
zeta[0]=zeta0;

Array[partnonres,z,0];
partnonres[0]=NormNRH3;

Array[partres,z,0];
partres[0]=NormRH3;

Array[chi,z,0];
chi[0]=chi0;

Array[theta,z,0];
theta[0]=theta0;

(* Size of the iterative normal form (components) *)

itsizenfL=0;
itsizenfG=0;
itsizenfl=0;
itsizenfg=0;

(* Gives a warning if the hypotheses of the iterative lemma are not fulfilled *)

If[T*Max[eta0L,eta0l,eta0G,eta0G]*z/2>1,Print["!!!!Warning!!!! The hypotheses of the iterative lemma are not fulfilled"]];

(* Iterations *)


For[j=1,j<z+1,j++,
delta[j]=Max[Abs[XRem2l[Linit-(3-2*j/z)*erreL]]/essel,Abs[XRem2l[Linit+(3-2*j/z)*erreL]]/essel];
etaL[j]=2*upsilonL[etaL[j-1],etal[j-1], etaG[j-1], etag[j-1] ,gammaL[j-1],gammal[j-1], gammaG[j-1], gammag[j-1],delta[j-1], T, z,chi[j-1],theta[j-1]]*etaL[j-1];
etaG[j]=2*upsilonG[etaL[j-1],etal[j-1], etaG[j-1], etag[j-1] ,gammaL[j-1],gammal[j-1], gammaG[j-1], gammag[j-1],delta[j-1], T, z,chi[j-1],theta[j-1]]*etaG[j-1];
etal[j]=2*upsilonl[etaL[j-1],etal[j-1], etaG[j-1], etag[j-1] ,gammaL[j-1],gammal[j-1], gammaG[j-1], gammag[j-1],delta[j-1], T, z,chi[j-1],theta[j-1]]*etal[j-1];
etag[j]=2*upsilong[etaL[j-1],etal[j-1], etaG[j-1], etag[j-1] ,gammaL[j-1],gammal[j-1], gammaG[j-1], gammag[j-1],delta[j-1], T, z,chi[j-1],theta[j-1]]*etag[j-1];
gammaL[j]=gammaL[j-1]+etaL[j]/2;
gammaG[j]=gammaG[j-1]+etaG[j]/2;
gammal[j]=gammal[j-1]+etal[j]/2;
gammag[j]=gammag[j-1]+etag[j]/2;
partnonres[j]=2*zeta[j-1]*partnonres[j-1];
partres[j]=partres[j-1]+partnonres[j]/2;
chi[j]=chifunc[etaL[j],etal[j], etaG[j], etag[j] ,gammaL[j],gammal[j], gammaG[j], gammag[j],delta[j]];
theta[j]=thetafunc[etaL[j],etal[j], etaG[j], etag[j] ,T];
zeta[j]=zeta[T,z,chi[j]];
itsizenfL=itsizenfL+T*etaL[j-1]/2*erreL;
itsizenfG=itsizenfG+T*etaG[j-1]/2*erreG;
itsizenfl=itsizenfl+T*etal[j-1]/2*essel;
itsizenfg=itsizenfg+T*etag[j-1]/2*esseg;
];


(* Gives a warning if the iterative lemma doesn't work *)

If[etaL[z]>etaL[0] && z≠0, Print["!!!!Warning!!!! The iterative lemma may not work well enough on the L variables!"]];
If[etaG[z]>etaG[0] && z≠0, Print["!!!!Warning!!!! The iterative lemma may not work well enough on the G variables!"]];
If[etal[z]>etal[0] && z≠0, Print["!!!!Warning!!!! The iterative lemma may not work well enough on the l variables!"]];
If[etag[z]>etag[0] && z≠0, Print["!!!!Warning!!!! The iterative lemma may not work well enough on the g variables!"]];


(* Gives a warning if the iterative lemma goes out of the right domains *)
If[Max[sizenfL/erreL,sizenfG/erreG,sizenfl/essel,sizenfg/esseg]>1,Print["ERROR! THE ITERATIVE NORMAL FORM GOES OUTSIDE OF THE ANALYTICITY DOMAIN!"]];

(* Print values after application of the iterative lemma *)

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

For[j=0,j<z+1,j++,
Print["etaL[",j,"]=",etaL[j]];
];

Print["---------------------------------------------------------"];
For[j=0,j<z+1,j++,
Print["etaG[",j,"]=",etaG[j]];
];

Print["---------------------------------------------------------"];
For[j=0,j<z+1,j++,
Print["etal[",j,"]=",etal[j]];
];

Print["---------------------------------------------------------"];
For[j=0,j<z+1,j++,
Print["etag[",j,"]=",etag[j]];
];

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

For[j=0,j<z+1,j++,
Print["gammaL[",j,"]=",gammaL[j]];
];

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

For[j=0,j<z+1,j++,
Print["gammaG[",j,"]=",gammaG[j]];
];

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

For[j=0,j<z+1,j++,
Print["gammal[",j,"]=",gammal[j]];
];

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

For[j=0,j<z+1,j++,
Print["gammag[",j,"]=",gammag[j]];
];

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

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

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

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

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

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

(* Initial radius in the normalized variables (since the initial condition is the point Linit, the initial radius is just the size of the normal form) *)
riLtilde=sizenfL+itsizenfL;
rfLtilde=2*erreL-sizenfL-itsizenfL;
rfL=2*erreL;

(* Inferior bound on the L eigenvalue of the hessian matrix of H0 *)
Dfreq1[LL_]=D[D[H0[LL,GG],LL],LL];
kgreek=Abs[Dfreq1[Linit+rfLtilde]];
Print["kgreek= ", kgreek];

(* Upper bound on the L eigenvalue of the hessian matrix of H0 *)
Dfreq1[LL_]=D[D[H0[LL,GG],LL],LL];
K=Abs[Dfreq1[Linit-riLtilde]];
Print["K= ", K];

(* Constants (see my article) *)
C3=kgreek/2*((rfLtilde-(K/kgreek+1)*riLtilde)^2-(K/kgreek*riLtilde)^2)-partnonres[z]-partres[z];
C4=Abs[freq1[Linit]]*etaL[z]*erreL+Abs[freq2]*etaG[z]*erreG;

(* Time of stability *)
tstab=C3/C4;

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

(* Variation in the action variables *)
relvarL=rfL/Linit; (* Maximal variation in the action variables*)
Print["Maximal relative variation in the action variable L= ", relvarL];