Haminiz = (px^2+py^2)/2+y*px-x*py -(1-mu)/Sqrt[(x-mu)^2+y^2]-mu/Sqrt[(x+1-mu)^2+y^2];
Sgen = px * (mu+rho*Cos[theta]) + (py - mu) * rho*Sin[theta];
Spx = D[Sgen, px];
Spy = D[Sgen,py];
Srho = D[Sgen,rho];
Stheta = D[Sgen,theta];
solsistvecchimomenti=Simplify[Expand[Solve[{prho== Srho , ptheta == Stheta}, {px,py}]]];
sostcoordinate = {{x -> Spx, y -> Spy}};
sostcompleta = Join[ solsistvecchimomenti[[1]] , sostcoordinate[[1]] ];
HamPolarElioc=Expand[Simplify[Haminiz /. sostcompleta]]

TreeForm[ HamPolarElioc ]
Head[ HamPolarElioc ];
Length[ HamPolarElioc ];
Hamsempl = HamPolarElioc;
Hamsempl -=  HamPolarElioc[[5]] + HamPolarElioc[[6]] + HamPolarElioc[[8]];
Hamsempl += - (1-mu)/rho

Ntronc = 5;
Hamsempl /. { prho -> PX, ptheta -> PY+1, rho -> X+1, theta -> Y - 2*Pi/3};
Series[%,{X,0,Ntronc}];
Series[%,{Y,0,Ntronc}];
Normal[%];
HamsemplPXY = ExpandAll[%];
MonomialList[ ExpandAll[%] , {PX,PY,X,Y}]

SvilPotPrim2 = Normal[ Series[ - mu / Sqrt[1 + u] , {u, 0, Ntronc} ] ];
uXYtmp =(2 * rho * Cos[theta] + rho^2) /. { rho -> X+1, theta -> Y - 2*Pi/3};
uXY = ExpandAll[ Normal[ Series[uXYtmp, {X,0,Ntronc}, {Y,0,Ntronc} ] ] ];
parteComplicataXY = Expand[ Normal[ SvilPotPrim2 /. {u -> uXY} ] ] ;
HamTotTmp = HamsemplPXY + parteComplicataXY;

SelPolOm[expr_,gr_]:=Module[{ris = 0, espan, nummon, j, addendo, gradd},
         espan = MonomialList[ ExpandAll[expr] , {PX,PY,X,Y} ];
	nummon = Length[ espan ];
	For[j=1, j≤nummon, j++,
	      addendo = espan[[j]];
               gradd = 0;
	      gradd += Exponent[addendo, PX ];
	      gradd += Exponent[addendo, PY ];
	      gradd += Exponent[addendo, X ];
	      gradd += Exponent[addendo, Y ];
	      If[ gradd == gr, ris += addendo ];
       ];
          ris]
HamPXY = 0;
For[ grpol = 0, grpol≤ Ntronc, grpol++,
HamPXY += SelPolOm[HamTotTmp, grpol]; ]
HamPXY
HamQuad = SelPolOm[HamTotTmp, 2];
matQ = D[ D[ HamQuad, {{PX, PY,X,Y}} ] , {{PX, PY,X,Y}} ];
matJ={{0,0,-1,0},{0,0,0,-1},{1,0,0,0},{0,1,0,0}};
matJQ = matJ . matQ;
CharacteristicPolynomial[ matJQ , lambda ]

matJQ=matJ.matQ;
matJQJQ=matJQ.matJQ;
Eigenvectors[matJQJQ];
autovet = Join[ {Eigenvectors[matJQJQ][[3]]},
                               {Eigenvectors[matJQJQ][[4]]},
                               {Eigenvectors[matJQJQ][[1]]},
                               {Eigenvectors[matJQJQ][[2]]} ];

ortonormsimpl[u1tmp_,u2tmp_]:=Module[{ris,u1,u2},u1=u1tmp/Norm[u1tmp];u2=u2tmp/(u2tmp.(matJ.u1));ris=Join[{u1},{u2}]]

autovetomega1quad =ortonormsimpl[autovet[[1]],autovet[[2]]];
autovetomega2quad =ortonormsimpl[autovet[[3]],autovet[[4]]];
matV = Transpose[
             Join[ { autovetomega1quad[[1]] } ,
                        { autovetomega2quad[[1]] } ,
                        { autovetomega1quad[[2]] } ,
                        { autovetomega2quad[[2]] } ] ];
matVtQV = Simplify[ Expand[ Transpose[matV] . matQ .matV ] ];
matLtmp = { {1, 0, beta1, 0},
               {0, 1, 0, beta2},
               {0, 0, 1, 0},
               {0, 0, 0, 1} };
matL = matLtmp /. { beta1 -> - matVtQV[[1]][[3]] / matVtQV[[1]][[1]] ,
beta2 -> - matVtQV[[2]][[4]] / matVtQV[[2]][[2]] };
matW = matV.matL;
matWtQW = Transpose[matL]. matVtQV . matL;
matD = { { (matWtQW[[3]][[3]]/matWtQW[[1]][[1]])^(1/4), 0, 0, 0 },
 { 0,(matWtQW[[4]][[4]]/matWtQW[[2]][[2]])^(1/4), 0, 0 },
 { 0, 0,(matWtQW[[1]][[1]]/matWtQW[[3]][[3]])^(1/4), 0 },
 { 0, 0, 0,(matWtQW[[2]][[2]]/matWtQW[[4]][[4]])^(1/4)} };
matU = matW . matD;
Simplify[ Expand[ Transpose[matU].matQ.matU ]] /. {mu -> 0.01}
Simplify[ Expand[ Transpose[matU].matJ.matU  - matJ ]]/. {mu -> 0.01}

rhslasttrasfcan =(matU/. {mu -> 0.01}) . {x1,x2,y1,y2};
fPX = rhslasttrasfcan[[1]];
fPY = rhslasttrasfcan[[2]];
fX = rhslasttrasfcan[[3]];
fY =  rhslasttrasfcan[[4]];
HamQuad/. {mu -> 0.01, PX -> fPX, PY -> fPY, X -> fX, Y -> fY};
Expand[%]

Hamxy = 0;
For[ grpol = 2, grpol≤ Ntronc, grpol++,
Hamtermtmp = SelPolOm[HamTotTmp, grpol];
Hamxy +=Hamtermtmp /. {mu -> 0.01, PX -> fPX, PY -> fPY, X -> fX, Y -> fY};
 ]
Hamxy;
ExpandAll[%]