
function [t,y]=solution(i,a,b,r,alpha,f,omegaP,omega0,theta0,tmax,te),
	c=cos(theta0),
	s=sin(theta0),
	Nz0=cos(alpha)
	sa=sin(alpha)
	function deriv=systeme(t,y),
		Nz=Nz0,
		e=y(10)*y(1)+y(11)*y(2)+y(12)*y(3),
		vpx=y(15)-(y(7)*y(1)+y(8)*y(2)+y(9)*y(3))*(y(12)+r)+e*y(9),
		vpy=y(16)+(y(4)*y(1)+y(5)*y(2)+y(6)*y(3))*(y(12)+r)-e*y(6),
		vp=sqrt(vpx*vpx+vpy*vpy),
		Tx=-f*Nz*vpx/vp,
		Ty=-f*Nz*vpy/vp,
		n1=-y(9)*Nz+(y(12)+r)*Ty,
		n2=y(6)*Nz-(y(12)+r)*Tx,
		n3=-y(6)*Ty+y(9)*Tx,
		deriv(1)=(1-i)*y(2)*y(3)+a*i*(y(4)*n1+y(7)*n2+y(10)*n3), // omegaS1
		deriv(2)=(i-1)*y(1)*y(3)+a*i*(y(5)*n1+y(8)*n2+y(11)*n3), // omegaS2
		deriv(3)=a*(y(6)*n1+y(9)*n2+y(12)*n3), // omegaS3
		deriv(4)=y(3)*y(5)-y(2)*y(6), // q11
		deriv(5)=y(1)*y(6)-y(3)*y(4), // q12
		deriv(6)=y(2)*y(4)-y(1)*y(5), // q13
		deriv(7)=y(3)*y(8)-y(2)*y(9), // q21
		deriv(8)=y(1)*y(9)-y(3)*y(7), // q22
		deriv(9)=y(2)*y(7)-y(1)*y(8), // q23
		deriv(10)=y(3)*y(11)-y(2)*y(12), // q31
		deriv(11)=y(1)*y(12)-y(3)*y(10), // q32
		deriv(12)=y(2)*y(10)-y(1)*y(11), // q33
		deriv(13)=y(15), // x
		deriv(14)=y(16), // y
		deriv(15)=b*(Tx+sa), // vx
		deriv(16)=b*Ty, // vy
	endfunction
	init=[0;omegaP*s;omega0+omegaP*c;1;0;0;0;c;-s;0;s;c;0;0;0;0];
	t=[0:te:tmax];
    tolA=1d-6;
    tolR=1d-9;
	%ODEOPTIONS=[1,0,0,%inf,0,2,10000,12,5,0,-1,-1];
	y=ode(init,0,t,tolR,tolA,systeme);
endfunction
			

function [a1,a2,a3]=calcul(y)
	a1 = y(1,:);
	a2 = a1
	a3 = a1;
	s=size(a1);
	for i=1:s(2),
		a1(i)=y(4,i)*y(4,i)+y(7,i)*y(7,i)+y(10,i)*y(10,i),
		a2(i)=y(12,i),
		a3(i)=y(6,i)/sqrt(1-y(12,i)*y(12,i)),
	end;
endfunction
			

tmax=10000;
i=0.2;
f=0.01;
a=2*%pi*0.01;
b=a*0.1;
r=0.01;
alpha=0;
[t,y]=solution(i,a,b,r,alpha,f,0,2*%pi,20*%pi/180,tmax,1);
[a1,a2,a3]=calcul(y);
			
plotA=scf();
plot2d(t,a1,style=5,rect=[0,0.9,tmax,1.1]);
xtitle('q11^2+q21^2+q22^2','t','');
			
plotB=scf();
plot2d(t,a2,style=5,rect=[0,0,tmax,1]);
xtitle('cos(theta)','t','');
			
plotC=scf();
plot2d(t,a3,style=5,rect=[0,-1,tmax,1]);
xtitle('sin(phi)','t','');
			
plotD=scf();
param3d(y(6,:),y(9,:),y(12,:),35,45,"x@y@z",[3,4],[-1,1,-1,1,-1,1]);
			
plotE=scf();
plot2d(y(13,:),y(14,:),style=5,rect=[-2,-4,2,0])
plotE.axes_size=[450,500];
xtitle('','xG','yG');
			
