function y=champmgn(u);

% U est le temps de la simulation

% Cette fonction donne les composantes approchées du champ magnétique terrestre  pour:
%	Une orbite circulaire d'altitude Z
%	Inclinaison orbitale i
% 	Passage à l'équateur au temps t=0 s
%
% Sur le repère orbital, à savoir :
%	Axe X tangent à l'orbite dans le même sens que la vitesse
%	Axe Y portant le moment orbital ou la vitesse de rotation
%	Axe Z verticale ascendante

% Caractéristiques de l'orbite du satellite

 altitude_orbite=817;		% altitude satellite en km
 rayon_orbital=(6378+altitude_orbite)*1000;
 vitesse_rotation=sqrt(39.86e13/rayon_orbital^3);
				% pulsation orbitale en rd/s
 inclinaison=98.7*pi/180;



% Caractéristiques du champ magnétique terrestre

 K=6.413e21;
 mu0=4e-7*pi;
 constante_magnetic=K*mu0/rayon_orbital^3;

% Calcul du champ magnétique terrestre

 t=u;

 periode_orbitale=2*pi/vitesse_rotation;

 angle_polaire=vitesse_rotation*t;

 latitude=asin(sin(inclinaison)*sin(angle_polaire));

 beta=real(asin(cos(inclinaison)./cos(latitude)));


 BX=constante_magnetic*sign(cos(angle_polaire)).*cos(latitude).*cos(beta);
 

 BY=constante_magnetic*cos(inclinaison)*sin(pi/2+t*1e-25);;

 
 BZ=-2*constante_magnetic*sin(inclinaison)*sin(vitesse_rotation*t);



  % Sortie des composantes

 
 y=[BX;BY;BZ];

 

  end