function y=gradgrav(u);

global Asat vitesse_rotation

% u est le vecteur d'entrée à 7 composantes
% 	u(1:3) = rotation absolue p q r en repère orbital
%	u(4;7) est la quaternion d'attitude du satellite par 
%	rapport au repère orbital

% Extraction du quaternion Q d'attitude
Q(1:4)=u(4:7);


% u est un quaternion et la fonction retourne y = Pmat la matrice de passage
% associée au changement de base


q0=Q(1);
q1=Q(2);
q2=Q(3);
q3=Q(4);

% Calcul de la matrice de passage

P11=2*(q0^2+q1^2)-1;P22=2*(q0^2+q2^2)-1;P33=2*(q0^2+q3^2)-1;
P21=2*(q1*q2+q0*q3);P12=2*(q1*q2-q0*q3);
P31=2*(q1*q3-q0*q2);P13=2*(q1*q3+q0*q2);
P32=2*(q2*q3+q0*q1);P23=2*(q2*q3-q0*q1);


% Calcul des composantes de Z sur les axes xyz du satellite

Z=[P31 P32 P33]';

% calcul du gradient

I_Z=Asat*Z;


vect(1:3)=Z(1:3);
vect(4:6)=I_Z(1:3);

grad=prodvect(vect);


grad=3*vitesse_rotation^2*grad;

y=[grad(1) grad(2) grad(3)]';



end