Inertia matrix from init_example3.dat

I1=1.000000000000000;
I2=1.012686988782515;
I3=3.306237422473038;

Initial angular momentum

L10=-0.609860759302936;
L20=0.761660947972381;
L30=0.218957654801698;

Data from quaternion code

Lq01=-0.544332842491675;
Lq02=0.729131780907662;
Lq03=-0.414811526666455;

Norm[{L10,L20,L30}]
Norm[{Lq01,Lq02,Lq03}]

Parameter d

d=L10^2/I1+L20^2/I2+L30^2/I3
dq=Lq01^2/I1+Lq02^2/I2+Lq03^2/I3
N[1/I2]

We have the case with D<1/I2

A1=Sqrt[I1*(d*I3-1)/(I3-I1)]
A2=Sqrt[I2*(d*I3-1)/(I3-I2)]
A3=Sqrt[I3*(1-d*I1)/(I3-I1)]
B=Sqrt[(1-d*I1)*(I3-I2)/(I1*I2*I3)]
m=((d*I3-1)*(I2-I1))/((1-d*I1)*(I3-I2))

This is function copied from qtraj1.nb

Q=Compile[{{t,_Real}},Module[{l1,l2,l3,Lp,ψ,T2,T1,T3,res},
l1=Re[A1*JacobiCN[B*t,m]];
l2=Re[A2*JacobiSN[B*t,m]];
l3=Re[A3*JacobiDN[B*t,m]];
Lp=Sqrt[l1^2+l2^2];
ψ=Re[(t/I3+((I3-I1)/(Sqrt[((I1(-1+d I1)I3 (I2-I3))/(I2 ))]))*EllipticPi[-(((I2-I1) I3)/(I1 (I3-I2))),JacobiAmplitude[t(Sqrt[(((1-d I1) (I3-I2))/(I1 I2 I3))]) ,m],m])];
T1={{((l1 l3 )/Lp),((l2 l3 )/Lp),-Lp },{-((l2 )/Lp),((l1 )/Lp),0},{l1,l2,l3}};
T2={{Cos[ψ],-Sin[ψ],0},{Sin[ψ],Cos[ψ],0},{0,0,1}};
T3=T2.T1;
res=T3]]

am=ArcTan[L10/A1,L20/A2]

Bt0=EllipticF[am,m]

t0=Bt0/B

Q0=Q[t0]

Inverse[Q0].Q[10+t0]

Chop[Re[A1*JacobiCN[B*t0,m]]-L10]
Chop[Re[A2*JacobiSN[B*t0,m]]-L20]
Chop[Re[A3*JacobiDN[B*t0,m]]-L30]

Transpose[Inverse[Q0].Q[10+t0]]