
Q={{q11[t],q12[t],q13[t]},{q21[t],q22[t],q23[t]},{q31[t],q32[t],q33[t]}};
n={-q23[t],q13[t],0};
ns=Transpose[Q].n;
fs={0,0,omegaS3[t]};
euler={omegaS1'[t]==(1-r)*omegaS2[t]*omegaS3[t]+a*r*ns[[1]]-b*r*fs[[1]],
       omegaS2'[t]==(r-1)*omegaS3[t]*omegaS1[t]+a*r*ns[[2]]-b*r*fs[[2]],
	   omegaS3'[t]==a*ns[[3]]-b*fs[[3]]}
			 

OmegaS={{0,-omegaS3[t],omegaS2[t]},{omegaS3[t],0,-omegaS1[t]},{-omegaS2[t],omegaS1[t],0}};
d = Q.OmegaS;
dQ = {q11'[t]==d[[1,1]],q12'[t]==d[[1,2]],q13'[t]==d[[1,3]],
      q21'[t]==d[[2,1]],q22'[t]==d[[2,2]],q23'[t]==d[[2,3]],
	  q31'[t]==d[[3,1]],q32'[t]==d[[3,2]],q33'[t]==d[[3,3]]}
			 

Q0=RotationMatrix[phi,{0,0,1}].RotationMatrix[theta,{1,0,0}].RotationMatrix[psi,{0,0,1}];
dQ0 = D[Q0/.{phi->phi[t],theta->theta[t],psi->psi[t]},t];
Qinit = Q0/.{phi->0,theta->theta0,psi->0};
dQinit = dQ0/.{phi'[t]->omegaP,psi'[t]->omega0,theta'[t]->0,phi[t]->0,
               theta[t]->theta0,psi[t]->0};
OmegaS = Transpose[Qinit].dQinit;
init = {omegaS1[0]==OmegaS[[3,2]],omegaS2[0]==OmegaS[[1,3]],omegaS3[0]==OmegaS[[2,1]],
        q11[0]==Qinit[[1,1]],q12[0]==Qinit[[1,2]],q13[0]==Qinit[[1,3]],
		q21[0]==Qinit[[2,1]],q22[0]==Qinit[[2,2]],q23[0]==Qinit[[2,3]],
		q31[0]==Qinit[[3,1]],q32[0]==Qinit[[3,2]],q33[0]==Qinit[[3,3]]}
			 

equations=Join[euler,dQ]
			 

systeme=Join[equations,init]/.{r->0.2,a->N[2*Pi*0.1],b->0,omegaP->0,omega0->N[2*Pi],
                               theta0->N[30*Pi/180]}
				

tmax=500;
variables={omegaS1,omegaS2,omegaS3,q11,q12,q13,q21,q22,q23,q31,q32,q33}; 
solution = NDSolve[systeme,variables,{t,0,tmax},AccuracyGoal->5,
                   Method->"ImplicitRungeKutta",MaxSteps->Infinity]
				
Plot[q11[t]*q11[t]+q21[t]*q21[t]+q31[t]*q31[t]/.solution,{t,0,tmax},
		PlotRange->{{0,tmax},{0.99,1.01}},AxesLabel->{"t","q11^2+q21^2+q31^2"}]
Plot[q33[t]/.solution,{t,0,tmax},
		PlotRange->{{0,tmax},{0,1}},AxesLabel->{"t","cos(theta)"},PlotPoints->10^4]
Plot[q13[t]/Sqrt[1-q33[t]*q33[t]]/.solution,{t,0,tmax},
		PlotRange->{{0,tmax},{-1,1}},AxesLabel->{"t","sin(phi)"},PlotPoints->10^4]
ParametricPlot3D[{q13[t],q23[t],q33[t]}/.solution,{t,0,100},PlotRange->{{-1,1},{-1,1},{-1,1}}]

systeme=Join[equations,init]/.{r->0.2,a->N[2*Pi*0.1],b->0.001,omegaP->0,
                               omega0->N[2*Pi],theta0->N[30*Pi/180]}
tmax=1000;
solution = NDSolve[systeme,variables,{t,0,tmax},AccuracyGoal->5,
                   Method->"ImplicitRungeKutta",MaxSteps->Infinity]
				
Plot[q11[t]*q11[t]+q21[t]*q21[t]+q31[t]*q31[t]/.solution,{t,0,tmax},
		PlotRange->{{0,tmax},{0.99,1.01}},AxesLabel->{"t","q11^2+q21^2+q31^2"}]
Plot[q33[t]/.solution,{t,0,tmax},
		PlotRange->{{0,tmax},{0,1}},AxesLabel->{"t","cos(theta)"},PlotPoints->10^4]
Plot[q13[t]/Sqrt[1-q33[t]*q33[t]]/.solution,{t,0,tmax},
		PlotRange->{{0,tmax},{-1,1}},AxesLabel->{"t","sin(phi)"},PlotPoints->10^4]
