w1 = 1;
w2 = N[(Sqrt[5] - 1) / 2,16];
Znf0 = -I * w1 * z1 * izb1 - I * w2 * z2* izb2;
R1 = 4; (* R1 è l'indice del resto *)
fs = Table[ 0 , R1 -1 ];
cubterm = x1^2*x2 - x2^3/3;
fs[[1]] = Expand[cubterm /. {x1 -> (z1-I*izb1)/N[Sqrt[2],16],  x2 -> (z2-I*izb2)/N[Sqrt[2],16]} ];
chi = Table[ 0 , R1 -1 ];

poisson[f_,g_]:=Module[{ris},ris=0;
ris=ris+D[f,z1]*D[g,izb1]-D[f,izb1]*D[g,z1];
ris=ris+D[f,z2]*D[g,izb2]-D[f,izb2]*D[g,z2];
ris];

omosolv[f_] :=
Module[{ris=0,termpert,listamonomi,j,jfin,addendo,l1,l2,l1b,l2b,denom},
	     termpert = ExpandAll[f];
	     listamonomi = MonomialList[termpert, {z1,z2,izb1,izb2}];
	     jfin = Length[listamonomi];
	     For[j=1, j<=jfin, j++,
		  addendo = listamonomi[[j]];
		  l1 = Exponent[addendo, z1];        l2  = Exponent[addendo, z2];
		  l1b = Exponent[addendo, izb1];   l2b =  Exponent[addendo, izb2];
		   denom = I * w1 * (l1-l1b) + I * w2 * (l2 - l2b);
		   If[ Abs[l1-l1b] + Abs[l2-l2b] > 0,
			ris -= addendo/denom; ]
	  ];
	    ris ]

passonorm[fs_, r_] :=
  Module[{ris, j, jfin, lieterm, summand,s, fix},
          ris = ExpandAll[fs];
	jfin = Floor[(R1 -1) / r];
         lieterm = Znf0;
	For[j=1, j≤jfin, j++,
	   summand = poisson[ lieterm, chi[[r]] ];
	    lieterm = summand / j;
	    ris[[j*r]] += lieterm;
      ];
	For[s=1, s<=R1-1, s++,
	    jfin = Floor[(R1 -1 - s) / r];
             lieterm = fs[[s]];
	    For[j=1, j≤jfin, j++,
	      summand = poisson[ lieterm, chi[[r]] ];
	      lieterm = summand / j;
	      ris[[s+j*r]] += lieterm;
          ];
     ];
        fix = ExpandAll[ ris ];
        fix[[r]] = Chop[ fix[[r]] , 10^-7];
        ris = fix
]

For[r=1, r≤R1-1, r++,
      chi[[r]] = omosolv[ fs[[r]] ];
      fs = passonorm[fs, r];
]

fs;
fs[[1]]
fs[[2]]
fs[[3]]

Znftot = Znf0 /. {izb1 -> I*az1, izb2 -> I*az2, z1 -> 1, z2 -> 1};
For[s=1, s<R1, s++,
Znftot += fs[[s]] /. {izb1 -> I*az1, izb2 -> I*az2, z1 -> 1, z2 -> 1};
]
Znftot
w1tot = D[ Znftot, az1 ]
w2tot = D[ Znftot, az2 ]

Znf0 = (z1 - I * izb1) / N[Sqrt[2],16];
fs=Table[0,{l,1,R1-1}];
For[r=1, r≤R1-1, r++,
      fs = passonorm[fs, r];
]
x1 = Znf0;
For[s=1, s<R1, s++,
x1 += fs[[s]];
]

Znf0 = (z2 - I * izb2) / N[Sqrt[2],16];
fs=Table[0,{l,1,R1-1}];
For[r=1, r≤R1-1, r++,
      fs = passonorm[fs, r];
]
x2 = Znf0;
For[s=1, s<R1, s++,
x2 += fs[[s]];
]

Znf0 = - I *(z1 + I * izb1) / N[Sqrt[2],16];
fs=Table[0,{l,1,R1-1}];
For[r=1, r≤R1-1, r++,
      fs = passonorm[fs, r];
]
y1 = Znf0;
For[s=1, s<R1, s++,
y1 += fs[[s]];
]

Znf0 = - I *(z2 + I * izb2) / N[Sqrt[2],16];
fs=Table[0,{l,1,R1-1}];
For[r=1, r≤R1-1, r++,
      fs = passonorm[fs, r];
]
y2 = Znf0;
For[s=1, s<R1, s++,
y2 += fs[[s]];
]

energia[x1_,x2_,y1_,y2_]:=w1*(x1^2+y1^2)/2+w2*(x2^2+y2^2)/2+x1^2*x2-x2^3/3;

az1star = 0.00016; (* Sono i valori di I1(0) e I2(0), con cui fissiamo *)
az2star = 0.00016; (* le condizioni iniziali,perché poniamo theta1(0)=theta2(0)=0 *)
w1star = N[ w1tot /. {az1 -> az1star, az2 -> az2star} , 16 ]
w2star = N[ w2tot /. {az1 -> az1star, az2 -> az2star} , 16 ]
Sqaz1star=Sqrt[az1star];
Sqaz2star=Sqrt[az2star];

(* Facciamo il test di conservazione dell'energia, per tempi t = 0, 0.1, 0.2, 0.3, ... , 100; Ci calcoliamo l'errore massimo *)
errmax = 0;
listpoints={};
For[t=0, t≤1000, t++,
z1t=Sqaz1star*Exp[-I*(w1star*t/10.)];z2t=Sqaz2star*Exp[-I*(w2star*t/10.)];
izb1t = I * Conjugate[z1t];
izb2t = I * Conjugate[z2t];
  x1t=x1/.{z1->z1t,z2->z2t,izb1->izb1t,izb2->izb2t};
x2t=x2/.{z1->z1t,z2->z2t,izb1->izb1t,izb2->izb2t};
y1t=y1/.{z1->z1t,z2->z2t,izb1->izb1t,izb2->izb2t};
y2t=y2/.{z1->z1t,z2->z2t,izb1->izb1t,izb2->izb2t};
Et = Re[ energia[x1t,x2t,y1t,y2t] ];
If[t == 0, 
E0 = Et; x10 = Re[x1t]; x20 = Re[x2t]; y10 = Re[y1t]; y20 = Re[y2t]; ];
If[ Abs[Et-E0] > errmax, errmax = Abs[Et-E0] ];
listpoints=Join[ listpoints,{{t/10.,Et-E0}}];
]
E0
errmax
x10
x20
y10
y20

ListLinePlot[ listpoints ]

N[2^(16)]

(9.99925139041130451/10- w1star) / (8.549940871072081 * 10^(-6))
(6.17988674052193154/10 - w2star) / (6.31665 * 10^(-7))
w1star
w2star