function y=laurens3(u);

global KB


% En entrée un vecteur u à 10 composantes:
%	u(1:3) est le champ magnétique pulsé élaboré à bord à l'aide des mesures
%	u(4:10) est le vecteur d'état p q r + le quaternion Q 

% --------------------------------------------------------------------------------------
% Cette fonction calcule la rotation absolue du satellite portée par le champ magnétique
% terrestre de la position satellite au temps courant t

B=u(1:3);    

NB=norm(B);                % Champ magnétique colonne 

omega_abs=u(4:6);            % Vecteur rotation absolu par ses coposantes sur les axes satellite x y z

% Calcul de la rotation portée par B

omegaB=(omega_abs'*B)/NB;    % Valeur algébrique sur l'unitaire de B, appelé u


% Fin calcul rotation portée par B  
% ----------------------------------------------------------------------------------------

% vecteur unitaire de B

u=B/NB;    % Toujours défini

% Angle alpha entre l'axe y et la direction de B, y est l'axe 2 donc la composante de B est B(2)

cos_alpha=B(2)/NB;              % Cet angle devra être surveillé au voisinage de pi/2

alpha=acos(cos_alpha);          % angle alpha en radian

% On se débarasse tout de suite du cas où l'axe y porte le champ ( alpha=0 ou pi ), donc cos_alpha=1 ou -1
 
if abs( cos_alpha)> 0.99

  cpl2_roue= -KB*u;             % Couple C2 colonne colinéaire à B

  M2=[0 0 0]';                  % Couple C1 en colonne, nul ou très faible, donc moment magnétique M2 nul

end    % Fin du if

% -----------------------------------------

if abs( cos_alpha)<= 0.99

	y=[0 1 0]';                          % axe y

	v1=-u(2)*y(3)+u(3)*y(2);              % On calcule y vectoriel u

	v2=-u(3)*y(1)+u(1)*y(3);

	v3=-u(1)*y(2)+u(2)*y(1);

	v=[v1 v2 v3]'/sin(alpha);            % Axe v résultat du produit vectoriel y pat u

   % Calcul de M2

   if abs( cos_alpha)>0.05

	M2=KB*omegaB*tan(alpha)*v/NB;

	cpl2_roue=-KB*omegaB/cos(alpha)*[0 1 0]';

     else

         M2=[0 0 0]';                         % moment magnétique nul y ( arrêt pour éviter la singularité )

	cpl2_roue=[0 0 0 ]';                 % Couple nul sur y ( arrêt pour éviter la singularité )

   end   % fin du if secondaire

end    % Fin du if primaire


%  Sorties des résultats

y(1:3)=M2;                                    % Moment magnétique à ajouter à M1

y(4:6)=cpl2_roue;                             % Copuple roue en dimension 3


end  % Fin programme