Twist coefficients for the periodic orbit of period 2 of the billiard map

1. Derivatives of the billiard map (s,u)->((Subscript[s, 1]),(Subscript[u, 1]))->((Subscript[s, 2]),(Subscript[u, 2]))

Functions for the map from point A to point B

Clear[rad0,rad1,tau,s1,u1];
(*set up the radius function*)
rad0[0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=rad0[0,t,ra,ra2,ra4,rb,rb2,rb4]=ra;
Derivative[1,0,0,0,0,0,0,0][rad0][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=0;
Derivative[2,0,0,0,0,0,0,0][rad0][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=ra2;
Derivative[3,0,0,0,0,0,0,0][rad0][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=0;
Derivative[4,0,0,0,0,0,0,0][rad0][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=ra4;
rad1[0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=rad1[0,t,ra,ra2,ra4,rb,rb2,rb4]=rb;
Derivative[1,0,0,0,0,0,0,0][rad1][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=0;
Derivative[2,0,0,0,0,0,0,0][rad1][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=rb2;
Derivative[3,0,0,0,0,0,0,0][rad1][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=0;
Derivative[4,0,0,0,0,0,0,0][rad1][0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=rb4;
(*set up the orbit length function and the billiard map*)tau[0,0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau[0,0,t,ra,ra2,ra4,rb,rb2,rb4]=t;
u1[0,0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=u1[0,0,t,ra,ra2,ra4,rb,rb2,rb4]=0;
s1[0,0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=s1[0,0,t,ra,ra2,ra4,rb,rb2,rb4]=0;
Derivative[1,0,0,0,0,0,0,0,0][s1][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad0[s,t,ra,ra2,ra4,rb,rb2,rb4]*Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])-((Sqrt[1-u^2])/(Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]));
Derivative[0,1,0,0,0,0,0,0,0][s1][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=((tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4])/(Sqrt[1-u^2]*Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]));
Derivative[1,0,0,0,0,0,0,0,0][u1][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad0[s,t,ra,ra2,ra4,rb,rb2,rb4]*rad1[s1[s,u,t,ra,ra2,ra4,rb,rb2,rb4],t,ra,ra2,ra4,rb,rb2,rb4])-Sqrt[1-u^2]/rad1[s1[s,u,t,ra,ra2,ra4,rb,rb2,rb4],t,ra,ra2,ra4,rb,rb2,rb4] -((Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])/(rad0[s,t,ra,ra2,ra4,rb,rb2,rb4]));
Derivative[0,1,0,0,0,0,0,0,0][u1][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad1[s1[s,u,t,ra,ra2,ra4,rb,rb2,rb4],t,ra,ra2,ra4,rb,rb2,rb4]*Sqrt[1-u^2])-((Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])/(Sqrt[1-u^2]));
Derivative[1,0,0,0,0,0,0,0,0][tau][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=u-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]*(tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad0[s,t,ra,ra2,ra4,rb,rb2,rb4]*Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])-((Sqrt[1-u^2])/(Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])));
Derivative[0,1,0,0,0,0,0,0,0][tau][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]*((tau[s,u,t,ra,ra2,ra4,rb,rb2,rb4])/(Sqrt[1-u^2]*Sqrt[1-u1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]));

Functions for the map from point B to point A

Clear[tau1,s2,u2];
tau1[0,0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau1[0,0,t,ra,ra2,ra4,rb,rb2,rb4]=t;
u2[0,0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=u2[0,0,t,ra,ra2,ra4,rb,rb2,rb4]=0;
s2[0,0,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=s2[0,0,t,ra,ra2,ra4,rb,rb2,rb4]=0;
Derivative[1,0,0,0,0,0,0,0,0][s2][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad1[s,t,ra,ra2,ra4,rb,rb2,rb4]*Sqrt[1-u2[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])-((Sqrt[1-u^2])/(Sqrt[1-u2[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]));
Derivative[0,1,0,0,0,0,0,0,0][s2][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=((tau1[s,u,t,ra,ra2,ra4,rb,rb2,rb4])/(Sqrt[1-u^2]*Sqrt[1-u2[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]));
Derivative[1,0,0,0,0,0,0,0,0][u2][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad1[s,t,ra,ra2,ra4,rb,rb2,rb4]*rad0[s2[s,u,t,ra,ra2,ra4,rb,rb2,rb4],t,ra,ra2,ra4,rb,rb2,rb4])-Sqrt[1-u^2]/rad0[s2[s,u,t,ra,ra2,ra4,rb,rb2,rb4],t,ra,ra2,ra4,rb,rb2,rb4] -((Sqrt[1-u2[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])/(rad1[s,t,ra,ra2,ra4,rb,rb2,rb4]));
Derivative[0,1,0,0,0,0,0,0,0][u2][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=tau1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad0[s2[s,u,t,ra,ra2,ra4,rb,rb2,rb4],t,ra,ra2,ra4,rb,rb2,rb4]*Sqrt[1-u^2])-((Sqrt[1-u2[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])/(Sqrt[1-u^2]));
Derivative[1,0,0,0,0,0,0,0,0][tau1][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=u-u2[s,u,t,ra,ra2,ra4,rb,rb2,rb4]*(tau1[s,u,t,ra,ra2,ra4,rb,rb2,rb4]/(rad1[s,t,ra,ra2,ra4,rb,rb2,rb4]*Sqrt[1-u2[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])-((Sqrt[1-u^2])/(Sqrt[1-u2[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2])));
Derivative[0,1,0,0,0,0,0,0,0][tau1][s_,u_,t_,ra_,ra2_,ra4_,rb_,rb2_,rb4_]:=-u2[s,u,t,ra,ra2,ra4,rb,rb2,rb4]*((tau1[s,u,t,ra,ra2,ra4,rb,rb2,rb4])/(Sqrt[1-u^2]*Sqrt[1-u2[s,u,t,ra,ra2,ra4,rb,rb2,rb4]^2]));

Evaluate along the orbit from A to B and then back to A

(*v={t,ra,ra2,ra4,rb,rb2,rb4}*)
Clear[t,ra,ra2,ra4,rb,rb2,rb4];
t[v__]:=t[v]=v[[1]];
ra[v__]:=ra[v]=v[[2]];
ra2[v__]:=ra2[v]=v[[3]];
ra4[v__]:=ra4[v]=v[[4]];
rb[v__]:=rb[v]=v[[5]];
rb2[v__]:=rb2[v]=v[[6]];
rb4[v__]:=rb4[v]=v[[7]];

Clear[scomp,ucomp,eta,lambda];
(*j-th derivative of s and k-th derivative of u at the periodic point O*)
scomp[j_,k_,v__]:=scomp[j,k,v]=(D[s2[s1[svariable,uvariable,t[v],ra[v],ra2[v],ra4[v],rb[v],rb2[v],rb4[v]],u1[svariable,uvariable,t[v],ra[v],ra2[v],ra4[v],rb[v],rb2[v],rb4[v]],t[v],ra[v],ra2[v],ra4[v],rb[v],rb2[v],rb4[v]],{svariable,j},{uvariable,k}]/.svariable->0)/.uvariable->0;
ucomp[j_,k_,v__]:=ucomp[j,k,v]=(D[u2[s1[svariable,uvariable,t[v],ra[v],ra2[v],ra4[v],rb[v],rb2[v],rb4[v]],u1[svariable,uvariable,t[v],ra[v],ra2[v],ra4[v],rb[v],rb2[v],rb4[v]],t[v],ra[v],ra2[v],ra4[v],rb[v],rb2[v],rb4[v]],{svariable,j},{uvariable,k}]/.svariable->0)/.uvariable->0;
eta[v__]:=eta[v]=(-scomp[0,1,v]/ucomp[1,0,v])^(1/4);
(*the preferred eigenvalue for the preferred eigenvector*)
lambda[v__]:=lambda[v]=scomp[1,0,v]+I*Sqrt[-scomp[0,1,v]*ucomp[1,0,v]];
(*only consider the case t<ra and t<rb*)
(*the general case: If[scomp[0,1,v]<0,scomp[1,0,v]+I*Sqrt[-scomp[0,1,v]*ucomp[1,0,v]],scomp[1,0,v]-I*Sqrt[-scomp[0,1,v]*ucomp[1,0,v]]]*)

2. Taylor expansion of the billiard map in the coordinate system (x,y)=(s/η, ηu)

Coefficients of (Subscript[x, 1])=∑(Subscript[a, jk])(x^j)(y^k), (Subscript[y, 1])=∑(Subscript[b, jk])(x^j)(y^k)

Clear[a10,a01,b10,b01,a30,a21,a12,a03,b30,b21,b12,b03,a50,a41,a32,a23,a14,a05,b50,b41,b32,b23,b14,b05];
a10[v__]:=a10[v]=scomp[1,0,v];a01[v__]:=a01[v]=eta[v]^(-2)*scomp[0,1,v];a30[v__]:=a30[v]=eta[v]^2*scomp[3,0,v]/6;a21[v__]:=a21[v]=scomp[2,1,v]/2;a12[v__]:=a12[v]=eta[v]^(-2)*scomp[1,2,v]/2;a03[v__]:=a03[v]=eta[v]^(-4)*scomp[0,3,v]/6;a50[v__]:=a50[v]=eta[v]^4*scomp[5,0,v]/5!;a41[v__]:=a41[v]=eta[v]^2*scomp[4,1,v]/24;a32[v__]:=a32[v]=scomp[3,2,v]/12;a23[v__]:=a23[v]=eta[v]^(-2)*scomp[2,3,v]/12;a14[v__]:=a14[v]=eta[v]^(-4)*scomp[1,4,v]/24;a05[v__]:=a05[v]=eta[v]^(-6)*scomp[0,5,v]/5!;
b10[v__]:=b10[v]=eta[v]^(2)*ucomp[1,0,v];b01[v__]:=b01[v]=ucomp[0,1,v];b30[v__]:=b30[v]=eta[v]^4*ucomp[3,0,v]/6;b21[v__]:=b21[v]=eta[v]^2*ucomp[2,1,v]/2;b12[v__]:=b12[v]=ucomp[1,2,v]/2;b03[v__]:=b03[v]=eta[v]^(-2)*ucomp[0,3,v]/6;b50[v__]:=b50[v]=eta[v]^6*ucomp[5,0,v]/5!;b41[v__]:=b41[v]=eta[v]^4*ucomp[4,1,v]/24;b32[v__]:=b32[v]=eta[v]^2*ucomp[3,2,v]/12;b23[v__]:=b23[v]=ucomp[2,3,v]/12;b14[v__]:=b14[v]=eta[v]^(-2)*ucomp[1,4,v]/24;b05[v__]:=b05[v]=eta[v]^(-4)*ucomp[0,5,v]/5!;

3. Coefficients of the coordinate transform (x,y)->(X,Y)=(∑(Subscript[p, jk])(x^j)(y^k),∑(Subscript[q, jk])(x^j)(y^k))

Clear[p30,p21,p12,p03,q30,q21,q12,q03,p50,p41,p32,p23,p14,p05,q50,q41,q32,q23,q14,q05];
p30[v__]:=p30[v]=-(2scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](a30[v]+b21[v]-a12[v]-b03[v])-(scomp[1,0,v]^2+scomp[0,1,v]*ucomp[1,0,v])(b30[v]-a21[v]-b12[v]+a03[v]))/(8Sqrt[-scomp[0,1,v]*ucomp[1,0,v]])+(Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](a30[v]-b21[v]-a12[v]+b03[v])+scomp[1,0,v](b30[v]+a21[v]-b12[v]-a03[v]))/(32scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]] );
p21[v__]:=p21[v]=3(Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](b30[v]+a21[v]-b12[v]-a03[v])-scomp[1,0,v](a30[v]-b21[v]-a12[v]+b03[v]))/(32scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]]);
p12[v__]:=p12[v]=-3(2scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](a30[v]+b21[v]-a12[v]-b03[v])-(scomp[1,0,v]^2+scomp[0,1,v]*ucomp[1,0,v])(b30[v]-a21[v]-b12[v]+a03[v]))/(8Sqrt[-scomp[0,1,v]*ucomp[1,0,v]])-3(Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](a30[v]-b21[v]-a12[v]+b03[v])+scomp[1,0,v](b30[v]+a21[v]-b12[v]-a03[v]))/(32scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]] );
p03[v__]:=p03[v]=(2scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](b30[v]-a21[v]-b12[v]+a03[v])+(scomp[1,0,v]^2+scomp[0,1,v]*ucomp[1,0,v])(a30[v]+b21[v]-a12[v]-b03[v]))/(4Sqrt[-scomp[0,1,v]*ucomp[1,0,v]])-(Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](b30[v]+a21[v]-b12[v]-a03[v])-scomp[1,0,v](a30[v]-b21[v]-a12[v]+b03[v]))/(32scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]] );
q30[v__]:=q30[v]=(2scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](b30[v]-a21[v]-b12[v]+a03[v])+(scomp[1,0,v]^2+scomp[0,1,v]*ucomp[1,0,v])(a30[v]+b21[v]-a12[v]-b03[v]))/(4Sqrt[-scomp[0,1,v]*ucomp[1,0,v]])+(Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](b30[v]+a21[v]-b12[v]-a03[v])-scomp[1,0,v](a30[v]-b21[v]-a12[v]+b03[v]))/(32scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]] );
q21[v__]:=q21[v]=3(2scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](a30[v]+b21[v]-a12[v]-b03[v])-(scomp[1,0,v]^2+scomp[0,1,v]*ucomp[1,0,v])(b30[v]-a21[v]-b12[v]+a03[v]))/(8Sqrt[-scomp[0,1,v]*ucomp[1,0,v]])-3(Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](a30[v]-b21[v]-a12[v]+b03[v])+scomp[1,0,v](b30[v]+a21[v]-b12[v]-a03[v]))/(32scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]] );
q12[v__]:=q12[v]=-3(Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](b30[v]+a21[v]-b12[v]-a03[v])-scomp[1,0,v](a30[v]-b21[v]-a12[v]+b03[v]))/(32scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]] );
q03[v__]:=q03[v]=(2scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](a30[v]+b21[v]-a12[v]-b03[v])-(scomp[1,0,v]^2+scomp[0,1,v]*ucomp[1,0,v])(b30[v]-a21[v]-b12[v]+a03[v]))/(8Sqrt[-scomp[0,1,v]*ucomp[1,0,v]])+(Sqrt[-scomp[0,1,v]*ucomp[1,0,v]](a30[v]-b21[v]-a12[v]+b03[v])+scomp[1,0,v](b30[v]+a21[v]-b12[v]-a03[v]))/(32scomp[1,0,v]Sqrt[-scomp[0,1,v]*ucomp[1,0,v]] );
p50[v__]:=p50[v]=p21[v]*q30[v];
p41[v__]:=p41[v]=p21[v]*q21[v]+2p12[v]*q30[v];
p32[v__]:=p32[v]=p21[v]*q12[v]+2p12[v]*q21[v]+3p03[v]*q30[v];
p23[v__]:=p23[v]=p21[v]*q03[v]+2p12[v]*q12[v]+3p03[v]*q21[v];
p14[v__]:=p14[v]=2p12[v]*q03[v]+3p03[v]*q12[v];
p05[v__]:=p05[v]=3p03[v]*q03[v];q50[v__]:=q50[v]=q21[v]*q30[v];
q41[v__]:=q41[v]=q21[v]*q21[v]+2q12[v]*q30[v];
q32[v__]:=q32[v]=q21[v]*q12[v]+2q12[v]*q21[v]+3q03[v]*q30[v];
q23[v__]:=q23[v]=q21[v]*q03[v]+2q12[v]*q12[v]+3q03[v]*q21[v];
q14[v__]:=q14[v]=2q12[v]*q03[v]+3q03[v]*q12[v];
q05[v__]:=q05[v]=3q03[v]*q03[v];

4. Taylor expansion of the billiard map (F^2) in terms of new coordinate system (X,Y): (X,Y)->((Subscript[X, 1]),(Subscript[Y, 1]))=(∑(Subscript[A, jk])(X^j)(Y^k),∑(Subscript[B, jk])(X^j)(Y^k))

Clear[aa30,aa21,aa12,aa03,bb30,bb21,bb12,bb03];
aa30[v__]:=aa30[v]=a30[v]+((b10[v])^3) p03[v]+a10[v] ((b10[v])^2) p12[v]+((a10[v])^2) b10[v]* p21[v]+((a10[v])^3) p30[v]-(a10[v]* p30[v]+a01[v]* q30[v]);
aa21[v__]:=aa21[v]=a21[v]+b01[v] (3 ((b10[v])^2) p03[v]+2 a10[v]* b10[v]* p12[v]+((a10[v])^2) p21[v])+a01[v] (((b10[v])^2) p12[v]+2 a10[v]* b10[v]* p21[v]+3 ((a10[v])^2) p30[v])-(a10[v]* p21[v]+a01[v]* q21[v]);aa12[v__]:=aa12[v]=a12[v]+((b01[v])^2) (3 b10[v]* p03[v]+a10[v]* p12[v])+2 a01[v]* b01[v] (b10[v]* p12[v]+a10[v]* p21[v])+((a01[v])^2) (b10[v]* p21[v]+3 a10[v]* p30[v])-(a10[v]* p12[v]+a01[v]* q12[v]);aa03[v__]:=aa03[v]=a03[v]+((b01[v])^3) p03[v]+a01[v] ((b01[v])^2) p12[v]+((a01[v])^2) b01[v]* p21[v]+((a01[v])^3) p30[v]-(a10[v]* p03[v]+a01[v]* q03[v]);
bb30[v__]:=bb30[v]=b30[v]+((b10[v])^3) q03[v]+a10[v] ((b10[v])^2) q12[v]+((a10[v])^2) b10[v]* q21[v]+((a10[v])^3) q30[v]-(b10[v]* p30[v]+b01[v]* q30[v]);
bb21[v__]:=bb21[v]=b21[v]+b01[v] (3 ((b10[v])^2) q03[v]+2 a10[v]* b10[v]* q12[v]+((a10[v])^2) q21[v])+a01[v] (((b10[v])^2) q12[v]+2 a10[v]* b10[v]* q21[v]+3 ((a10[v])^2) q30[v])-(b10[v]* p21[v]+b01[v]* q21[v]);bb12[v__]:=bb12[v]=b12[v]+((b01[v])^2) (3 b10[v]* q03[v]+a10[v]* q12[v])+2 a01[v]* b01[v] (b10[v]* q12[v]+a10[v]* q21[v])+((a01[v])^2) (b10[v]* q21[v]+3 a10[v]* q30[v])-(b10[v]* p12[v]+b01[v]* q12[v]);bb03[v__]:=bb03[v]=b03[v]+((b01[v])^3) q03[v]+a01[v] ((b01[v])^2) q12[v]+((a01[v])^2) b01[v]* q21[v]+((a01[v])^3) q30[v]-(b10[v]* p03[v]+b01[v]* q03[v]);

Clear[aa50,aa41,aa32,aa23,aa14,aa05];
aa50[v__]:=aa50[v]=a50[v]+3 ((b10[v])^2) b30[v] p03[v]+((b10[v])^5) p05[v]+(a30[v] ((b10[v])^2)+2 a10[v] b10[v] b30[v]) p12[v]+a10[v] ((b10[v])^4) p14[v]+(2 a10[v] a30[v] b10[v]+((a10[v])^2) b30[v]) p21[v]+((a10[v])^2) ((b10[v])^3) p23[v]+3 ((a10[v])^2) a30[v] p30[v]+((a10[v])^3) ((b10[v])^2) p32[v]+((a10[v])^4) b10[v] p41[v]+((a10[v])^5) p50[v]-(3 aa30[v] p30[v]+a10[v] p50[v]+aa21[v] q30[v]+a01[v] q50[v]);
aa41[v__]:=aa41[v]=a41[v]+((b10[v])^4) (5 b01[v] p05[v]+a01[v] p14[v])+2 a10[v] ((b10[v])^3) (2 b01[v] p14[v]+a01[v] p23[v])+((b10[v])^2) (3 b21[v] p03[v]+a21[v] p12[v]+3 ((a10[v])^2) (b01[v] p23[v]+a01[v] p32[v]))+2 b10[v] (a10[v] b21[v] p12[v]+a01[v] b30[v] p12[v]+a10[v] a21[v] p21[v]+a01[v] a30[v] p21[v]+b01[v] (3 b30[v] p03[v]+a30[v] p12[v]+((a10[v])^3) p32[v])+2 a01[v] ((a10[v])^3) p41[v])+a10[v] (a10[v] b21[v] p21[v]+2 a01[v] b30[v] p21[v]+3 a10[v] a21[v] p30[v]+6 a01[v] a30[v] p30[v]+b01[v] (2 b30[v] p12[v]+2 a30[v] p21[v]+((a10[v])^3) p41[v])+5 a01[v] ((a10[v])^3) p50[v])-(3 aa30[v] p21[v]+2 aa21[v] p30[v]+a10[v] p41[v]+aa21[v] q21[v]+2 aa12[v] q30[v]+a01[v] q41[v]);
aa32[v__]:=aa32[v]=a32[v]+3 ((b01[v])^2) b30[v] p03[v]+a30[v] ((b01[v])^2) p12[v]+2 a10[v] b01[v] b21[v] p12[v]+2 a01[v] b01[v] b30[v] p12[v]+2 a10[v] a21[v] b01[v] p21[v]+2 a01[v] a30[v] b01[v] p21[v]+((a10[v])^2) b12[v] p21[v]+2 a01[v] a10[v] b21[v] p21[v]+((a01[v])^2) b30[v] p21[v]+((b10[v])^3) (10 ((b01[v])^2) p05[v]+4 a01[v] b01[v] p14[v]+((a01[v])^2) p23[v])+3 ((a10[v])^2) a12[v] p30[v]+6 a01[v] a10[v] a21[v] p30[v]+3 ((a01[v])^2) a30[v] p30[v]+((a10[v])^3) ((b01[v])^2) p32[v]+((b10[v])^2) (3 b12[v] p03[v]+a12[v] p12[v]+6 a10[v] ((b01[v])^2) p14[v]+6 a01[v] a10[v] b01[v] p23[v]+3 ((a01[v])^2) a10[v] p32[v])+4 a01[v] ((a10[v])^3) b01[v] p41[v]+b10[v] (3 ((a10[v])^2) ((b01[v])^2) p23[v]+2 b01[v] (3 b21[v] p03[v]+a21[v] p12[v]+3 a01[v] ((a10[v])^2) p32[v])+2 (a10[v] b12[v] p12[v]+a01[v] b21[v] p12[v]+a10[v] a12[v] p21[v]+a01[v] a21[v] p21[v]+3 ((a01[v])^2) ((a10[v])^2) p41[v]))+10 ((a01[v])^2) ((a10[v])^3) p50[v]-(3 aa30[v] p12[v]+2 aa21[v] p21[v]+aa12[v] p30[v]+a10[v] p32[v]+aa21[v] q12[v]+2 aa12[v] q21[v]+3 aa03[v] q30[v]+a01[v] q32[v]);
aa23[v__]:=aa23[v]=a23[v]+6 b01[v] b10[v] b12[v] p03[v]+3 ((b01[v])^2) b21[v] p03[v]+10 ((b01[v])^3) ((b10[v])^2) p05[v]+a21[v] ((b01[v])^2) p12[v]+2 a12[v] b01[v] b10[v] p12[v]+a03[v] ((b10[v])^2) p12[v]+2 a10[v] b01[v] b12[v] p12[v]+2 a01[v] b10[v] b12[v] p12[v]+2 a01[v] b01[v] b21[v] p12[v]+4 a10[v] ((b01[v])^3) b10[v] p14[v]+6 a01[v] ((b01[v])^2) ((b10[v])^2) p14[v]+2 a10[v] a12[v] b01[v] p21[v]+2 a01[v] a21[v] b01[v] p21[v]+2 a03[v] a10[v] b10[v] p21[v]+2 a01[v] a12[v] b10[v] p21[v]+2 a01[v] a10[v] b12[v] p21[v]+((a01[v])^2) b21[v] p21[v]+b03[v] (3 ((b10[v])^2) p03[v]+2 a10[v] b10[v] p12[v]+((a10[v])^2) p21[v])+((a10[v])^2) ((b01[v])^3) p23[v]+6 a01[v] a10[v] ((b01[v])^2) b10[v] p23[v]+3 ((a01[v])^2) b01[v] ((b10[v])^2) p23[v]+3 a03[v] ((a10[v])^2) p30[v]+6 a01[v] a10[v] a12[v] p30[v]+3 ((a01[v])^2) a21[v] p30[v]+3 a01[v] ((a10[v])^2) ((b01[v])^2) p32[v]+6 ((a01[v])^2) a10[v] b01[v] b10[v] p32[v]+((a01[v])^3) ((b10[v])^2) p32[v]+6 ((a01[v])^2) ((a10[v])^2) b01[v] p41[v]+4 ((a01[v])^3) a10[v] b10[v] p41[v]+10 ((a01[v])^3) ((a10[v])^2) p50[v]-(3 aa30[v] p03[v]+2 aa21[v] p12[v]+aa12[v] p21[v]+a10[v] p23[v]+aa21[v] q03[v]+2 aa12[v] q12[v]+3 aa03[v] q21[v]+a01[v] q23[v]);
aa14[v__]:=aa14[v]=a14[v]+((b01[v])^4) (5 b10[v] p05[v]+a10[v] p14[v])+2 a01[v] ((b01[v])^3) (2 b10[v] p14[v]+a10[v] p23[v])+((b01[v])^2) (3 b12[v] p03[v]+a12[v] p12[v]+3 ((a01[v])^2) (b10[v] p23[v]+a10[v] p32[v]))+2 b01[v] (3 b03[v] b10[v] p03[v]+a10[v] b03[v] p12[v]+a03[v] b10[v] p12[v]+a01[v] b12[v] p12[v]+a03[v] a10[v] p21[v]+a01[v] a12[v] p21[v]+((a01[v])^3) b10[v] p32[v]+2 ((a01[v])^3) a10[v] p41[v])+a01[v] (2 a03[v] b10[v] p21[v]+a01[v] b12[v] p21[v]+2 b03[v] (b10[v] p12[v]+a10[v] p21[v])+6 a03[v] a10[v] p30[v]+3 a01[v] a12[v] p30[v]+((a01[v])^3) b10[v] p41[v]+5 ((a01[v])^3) a10[v] p50[v])-(2 aa21[v] p03[v]+aa12[v] p12[v]+a10[v] p14[v]+2 aa12[v] q03[v]+3 aa03[v] q12[v]+a01[v] q14[v]);
aa05[v__]:=aa05[v]=a05[v]+((b01[v])^5) p05[v]+a01[v] ((b01[v])^4) p14[v]+((a01[v])^2) ((b01[v])^3) p23[v]+((b01[v])^2) (3 b03[v] p03[v]+a03[v] p12[v]+((a01[v])^3) p32[v])+a01[v] b01[v] (2 b03[v] p12[v]+2 a03[v] p21[v]+((a01[v])^3) p41[v])+((a01[v])^2) (b03[v] p21[v]+3 a03[v] p30[v]+((a01[v])^3) p50[v])-(aa12[v] p03[v]+a10[v] p05[v]+3 aa03[v] q03[v]+a01[v] q05[v]);

Clear[bb50,bb41,bb32,bb23,bb14,bb05];
bb50[v__]:=bb50[v]=b50[v]+3 ((b10[v])^2) b30[v] q03[v]+((b10[v])^5) q05[v]+(a30[v] ((b10[v])^2)+2 a10[v] b10[v] b30[v]) q12[v]+a10[v] ((b10[v])^4) q14[v]+(2 a10[v] a30[v] b10[v]+((a10[v])^2) b30[v]) q21[v]+((a10[v])^2) ((b10[v])^3) q23[v]+3 ((a10[v])^2) a30[v] q30[v]+((a10[v])^3) ((b10[v])^2) q32[v]+((a10[v])^4) b10[v] q41[v]+((a10[v])^5) q50[v]-(3 bb30[v] p30[v]+b10[v] p50[v]+bb21[v] q30[v]+b01[v] q50[v]);
bb41[v__]:=bb41[v]=b41[v]+((b10[v])^4) (5 b01[v] q05[v]+a01[v] q14[v])+2 a10[v] ((b10[v])^3) (2 b01[v] q14[v]+a01[v] q23[v])+((b10[v])^2) (3 b21[v] q03[v]+a21[v] q12[v]+3 ((a10[v])^2) (b01[v] q23[v]+a01[v] q32[v]))+2 b10[v] (a10[v] b21[v] q12[v]+a01[v] b30[v] q12[v]+a10[v] a21[v] q21[v]+a01[v] a30[v] q21[v]+b01[v] (3 b30[v] q03[v]+a30[v] q12[v]+((a10[v])^3) q32[v])+2 a01[v] ((a10[v])^3) q41[v])+a10[v] (a10[v] b21[v] q21[v]+2 a01[v] b30[v] q21[v]+3 a10[v] a21[v] q30[v]+6 a01[v] a30[v] q30[v]+b01[v] (2 b30[v] q12[v]+2 a30[v] q21[v]+((a10[v])^3) q41[v])+5 a01[v] ((a10[v])^3) q50[v])-(3 bb30[v] p21[v]+2 bb21[v] p30[v]+b10[v] p41[v]+bb21[v] q21[v]+2 bb12[v] q30[v]+b01[v] q41[v]);
bb32[v__]:=bb32[v]=b32[v]+3 ((b01[v])^2) b30[v] q03[v]+a30[v] ((b01[v])^2) q12[v]+2 a10[v] b01[v] b21[v] q12[v]+2 a01[v] b01[v] b30[v] q12[v]+2 a10[v] a21[v] b01[v] q21[v]+2 a01[v] a30[v] b01[v] q21[v]+((a10[v])^2) b12[v] q21[v]+2 a01[v] a10[v] b21[v] q21[v]+((a01[v])^2) b30[v] q21[v]+((b10[v])^3) (10 ((b01[v])^2) q05[v]+4 a01[v] b01[v] q14[v]+((a01[v])^2) q23[v])+3 ((a10[v])^2) a12[v] q30[v]+6 a01[v] a10[v] a21[v] q30[v]+3 ((a01[v])^2) a30[v] q30[v]+((a10[v])^3) ((b01[v])^2) q32[v]+((b10[v])^2) (3 b12[v] q03[v]+a12[v] q12[v]+6 a10[v] ((b01[v])^2) q14[v]+6 a01[v] a10[v] b01[v] q23[v]+3 ((a01[v])^2) a10[v] q32[v])+4 a01[v] ((a10[v])^3) b01[v] q41[v]+b10[v] (3 ((a10[v])^2) ((b01[v])^2) q23[v]+2 b01[v] (3 b21[v] q03[v]+a21[v] q12[v]+3 a01[v] ((a10[v])^2) q32[v])+2 (a10[v] b12[v] q12[v]+a01[v] b21[v] q12[v]+a10[v] a12[v] q21[v]+a01[v] a21[v] q21[v]+3 ((a01[v])^2) ((a10[v])^2) q41[v]))+10 ((a01[v])^2) ((a10[v])^3) q50[v]-(3 bb30[v] p12[v]+2 bb21[v] p21[v]+bb12[v] p30[v]+b10[v] p32[v]+bb21[v] q12[v]+2 bb12[v] q21[v]+3 bb03[v] q30[v]+b01[v] q32[v]);
bb23[v__]:=bb23[v]=b23[v]+6 b01[v] b10[v] b12[v] q03[v]+3 ((b01[v])^2) b21[v] q03[v]+10 ((b01[v])^3) ((b10[v])^2) q05[v]+a21[v] ((b01[v])^2) q12[v]+2 a12[v] b01[v] b10[v] q12[v]+a03[v] ((b10[v])^2) q12[v]+2 a10[v] b01[v] b12[v] q12[v]+2 a01[v] b10[v] b12[v] q12[v]+2 a01[v] b01[v] b21[v] q12[v]+4 a10[v] ((b01[v])^3) b10[v] q14[v]+6 a01[v] ((b01[v])^2) ((b10[v])^2) q14[v]+2 a10[v] a12[v] b01[v] q21[v]+2 a01[v] a21[v] b01[v] q21[v]+2 a03[v] a10[v] b10[v] q21[v]+2 a01[v] a12[v] b10[v] q21[v]+2 a01[v] a10[v] b12[v] q21[v]+((a01[v])^2) b21[v] q21[v]+b03[v] (3 ((b10[v])^2) q03[v]+2 a10[v] b10[v] q12[v]+((a10[v])^2) q21[v])+((a10[v])^2) ((b01[v])^3) q23[v]+6 a01[v] a10[v] ((b01[v])^2) b10[v] q23[v]+3 ((a01[v])^2) b01[v] ((b10[v])^2) q23[v]+3 a03[v] ((a10[v])^2) q30[v]+6 a01[v] a10[v] a12[v] q30[v]+3 ((a01[v])^2) a21[v] q30[v]+3 a01[v] ((a10[v])^2) ((b01[v])^2) q32[v]+6 ((a01[v])^2) a10[v] b01[v] b10[v] q32[v]+((a01[v])^3) ((b10[v])^2) q32[v]+6 ((a01[v])^2) ((a10[v])^2) b01[v] q41[v]+4 ((a01[v])^3) a10[v] b10[v] q41[v]+10 ((a01[v])^3) ((a10[v])^2) q50[v]-(3 bb30[v] p03[v]+2 bb21[v] p12[v]+bb12[v] p21[v]+b10[v] p23[v]+bb21[v] q03[v]+2 bb12[v] q12[v]+3 bb03[v] q21[v]+b01[v] q23[v]);
bb14[v__]:=bb14[v]=b14[v]+((b01[v])^4) (5 b10[v] q05[v]+a10[v] q14[v])+2 a01[v] ((b01[v])^3) (2 b10[v] q14[v]+a10[v] q23[v])+((b01[v])^2) (3 b12[v] q03[v]+a12[v] q12[v]+3 ((a01[v])^2) (b10[v] q23[v]+a10[v] q32[v]))+2 b01[v] (3 b03[v] b10[v] q03[v]+a10[v] b03[v] q12[v]+a03[v] b10[v] q12[v]+a01[v] b12[v] q12[v]+a03[v] a10[v] q21[v]+a01[v] a12[v] q21[v]+((a01[v])^3) b10[v] q32[v]+2 ((a01[v])^3) a10[v] q41[v])+a01[v] (2 a03[v] b10[v] q21[v]+a01[v] b12[v] q21[v]+2 b03[v] (b10[v] q12[v]+a10[v] q21[v])+6 a03[v] a10[v] q30[v]+3 a01[v] a12[v] q30[v]+((a01[v])^3) b10[v] q41[v]+5 ((a01[v])^3) a10[v] q50[v])-(2 bb21[v] p03[v]+bb12[v] p12[v]+b10[v] p14[v]+2 bb12[v] q03[v]+3 bb03[v] q12[v]+b01[v] q14[v]);
bb05[v__]:=bb05[v]=b05[v]+3 ((b01[v])^2) b03[v] q03[v]+((b01[v])^5) q05[v]+(a03[v] ((b01[v])^2)+2 a01[v] b01[v] b03[v]) q12[v]+a01[v] ((b01[v])^4) q14[v]+(2 a01[v] a03[v] b01[v]+((a01[v])^2) b03[v]) q21[v]+((a01[v])^2) ((b01[v])^3) q23[v]+3 ((a01[v])^2) a03[v] q30[v]+((a01[v])^3) ((b01[v])^2) q32[v]+((a01[v])^4) b01[v] q41[v]+((a01[v])^5) q50[v]-(bb12[v] p03[v]+b10[v] p05[v]+3 bb03[v] q03[v]+b01[v] q05[v]);

5. Twist Coefficients

Clear[twist1,twist2];
twist1[v__]:=twist1[v]=(a10[v](3b30[v]-a21[v]+b12[v]-3a03[v])-Sqrt[-a01[v]*b10[v]](3a30[v]+b21[v]+a12[v]+3b03[v]))/8;
twist2[v__]:=twist2[v]=(a10[v](-2aa41[v]-2aa23[v]-10aa05[v]+10bb50[v]+2bb32[v]+2bb14[v])-Sqrt[-a01[v]*b10[v]](10aa50[v]+2aa32[v]+2aa14[v]+2bb41[v]+2bb23[v]+10bb05[v]))/32;