
function [t,y]=solution(r,a,b,omegaP,omega0,theta0,tmax,te),
	c=cos(theta0),
	s=sin(theta0),
	function deriv=systeme(t,y),
		deriv(1)=(1-r)*y(2)*y(3)+a*r*(y(6)*y(7)-y(4)*y(9)), // omegaS1
		deriv(2)=(r-1)*y(1)*y(3)+a*r*(y(6)*y(8)-y(5)*y(9)), // omegaS2
		deriv(3)=-b*y(3), // 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
	endfunction
	init=[0;omegaP*s;omega0+omegaP*c;1;0;0;0;c;-s;0;s;c];
	t=[0:te:tmax];
    tolA=1d-7;
    tolR=1d-10;
	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=500;
[t,y]=solution(0.2,2*%pi*0.1,0,0,2*%pi,30*%pi/180,tmax,0.1);
[a1,a2,a3]=calcul(y);
			
plotH=scf();
plot2d(t,a1,style=5,rect=[0,0.99,tmax,1.01]);
xtitle('q11^2+q21^2+q22^2','t','');
			
plotI=scf();
plot2d(t,a2,style=5,rect=[0,0,tmax,1]);
xtitle('cos(theta)','t','');
			
plotJ=scf();
plot2d(t,a3,style=5,rect=[0,-1,tmax,1]);
xtitle('sin(phi)','t','');
			
plotK=scf();
param3d(y(6,:),y(9,:),y(12,:),35,45,"x@y@z",[3,4],[-1,1,-1,1,-1,1]);
			
