//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [pos_car,vel_car,jacob] = CL_oe_cirEqua2car(cir_equa, mu)
// Circular equatorial adapted to cartesian orbital elements - DEPRECATED
//
// Calling Sequence
// [pos_car,vel_car,jacob]=CL_oe_cirEqua2car(cir_equa [,mu])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p>Replacement function: <link linkend="CL_oe_cireq2car">CL_oe_cireq2car</link></p>
// <p>"cirEqua" orbital elements are identical to "cireq" orbital elements except that 
// ix=2*sin(i/2)*cos(raan) and iy=2*sin(i/2)*sin(raan)</p>
// <p></p></listitem>
// <listitem>
// <p>Converts orbital elements adapted to near-circular near-equatorial orbits to 
// position and velocity for an elliptic orbit.</p>
// <p>The transformation jacobian is also optionally computed.</p>
// <p>See <link linkend="Orbital elements">Orbital elements</link> for more details</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// cir_equa: orbital elements adapted to near-circular near-equatorial orbits [sma;ex;ey;ix;iy;pom+raan+anm] [m,rad] (6xN)
// mu : (optional) gravitational constant [m^3/s^2] (default value is %CL_mu)
// pos_car: position [X;Y;Z] [m] (3xN)
// vel_car: velocity [Vx;Vy;Vz] [m/s] (3xN)
// jacob: (optional) transformation jacobian (See <link linkend="Orbital elements">Orbital elements</link> for more details) (6x6xN)
//
// Bibliography
// 1) CNES - MSLIB FORTRAN 90, Volume V (mv_cir_equa_car)
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Example 1
// cir_equa=[42166.712;-7.9d-6;1.1d-4;1.2d-4;-1.16d-4;5.3];
// [pos_car,vel_car,jacob1] = CL_oe_cirEqua2car(cir_equa);
// [cirequa2,jacob2] = CL_oe_car2cirEqua(pos_car,vel_car);

//!************************************************************************
//! But: Passage des parametres orbitaux dits adaptes a l'orbite CIRculaire EQUAtoriale, aux parametres CARtesiens
//! ===
//!
//! Note d'utilisation:  1) Applicable aux orbites circulaires ou elliptiques, equatoriales ou non.
//! ==================   2) Non applicable aux orbites hyperboliques ou paraboliques.
//!                      3) Pour les cas hyperboliques ou paraboliques, voir mv_kep_car.
//!                      4) L'element (i,j) de la jacobienne correspond a la derivee partielle du parametre cartesien numero i
//!                         par rapport a la derivee partielle du parametre dit adapte a l'orbite circulaire equatoriale numero j.
//!                      5) Les unites en entree doivent etre coherentes entre elles. Ex.:  mu en m**3/s**2, cir_equa en metres et
//!                         radians,ainsi pos_car sera en metres et vel_car en m/s.
//!                      6) La routine mv_car_cir_equa effectue la transformation inverse.
//!************************************************************************


// Declarations:
EPS_ORB = CL__dataGetEnv("epsOrb", internal=%t);


// Code:

CL__warnDeprecated(); // deprecated function

[lhs,rhs] = argn()
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

eps100 = 100.*%eps

N = size(cir_equa,2)

a   = cir_equa(1,:)
ex  = cir_equa(2,:)
ey  = cir_equa(3,:)
hx  = cir_equa(4,:)
hy  = cir_equa(5,:)
alm = cir_equa(6,:)

pos_car = [];
vel_car = [];
//==========================================================
//test des parametres d'entree
//==========================================================

ex2 = ex.^2
ey2 = ey.^2
e2 = ex2 + ey2

ii = find(sqrt(e2)>1-EPS_ORB.parab)
if ii~=[] then CL__error('parabolic or hyperbolic orbit'); end

ii = find(a <= eps100)
if ii~=[]
  CL__error('semi major axis negatif or null');
end

hx2 = hx.^2
hy2 = hy.^2
h2 = hx2 + hy2
ii = find(h2 > 4)
if ii~=[]
  CL__error('inclination vector norm too big');
end
ii = find(h2 > (4*(1-eps100)))
if ii~=[]                                                 //la definition du vecteur inclinaison (hx,hy) impose que l'inclinaison i soit
  CL__error('inclination bigger than (pi - epsilon_retro)'); //inferieure a (pi - epsilon_retro) -> test si h**2 <= 4*(cos(epsilon_retro/2)**2)
end                                                     //soit h**2 <= 4*(1-100*epsilon) si epsilon_retro = 20*sqrt(epsilon)

//calcul ale (=pomega+gomega+E) par equation kepler generalisee

ale = CL_kp_M2Ecir(ex,ey,alm) //call mv_kepler_gene(alm, ex, ey, ale, code_retour)

cle = cos(ale)
sle = sin(ale)

//*************************************************************
// 1- calcul vecteur position r = a (ld2*x1 + mu2*y1)
//        et vecteur vitesse rv = sqrt(mu/a) (nu2*x1 + ki2*y1)
//*************************************************************

//calcul de variables intermediaires

sq = sqrt(1-e2)   //calcul de sqrt(1-ex**2-ey**2)
                  //remarque: 1-e2 > 0 car on a teste si e etait dans [0.,1.[
exey = ex.*ey

ecceys = (ex.*cle) + (ey.*sle)   //! ->  ex*cos( petit omega + grand omega + E) + ey*sin ( petit omega + grand omega + E)
eycexs = (ey.*cle) - (ex.*sle)   //! ->  ey*cos( petit omega + grand omega + E) - ex*sin ( petit omega + grand omega + E)

hxhy = hx.*hy //calcul variable sqh = sqrt(1._pm_reel-((hx2+hy2)/4._pm_reel))

sqh = sqrt(1 - (h2/4))  //remarque: 1-h2/4 > 0 car on a teste si h2 etait > 4.

x1(1,:) = 1 - (0.5.*hy2)
x1(2,:) = 0.5.*hxhy //calcul des composantes du vecteurs x1 en fonction de hx,hy
x1(3,:) = -hy.*sqh

y1(1,:) = 0.5.*hxhy
y1(2,:) = 1 - (0.5.*hx2)  //calcul des composantes du vecteurs y1 en fonction de hx,hy
y1(3,:) = hx.*sqh

betaa = 1.0 ./ (1 + sq) //calcul variable betaa = 1 / 1+sqrt(1-ex2-ey2)

rld2 = -ex + ((1 - (betaa.*ey2)).*cle) + (betaa.*exey.*sle) //calcul de ld2
rmu2 = -ey + ((1 - (betaa.*ex2)).*sle) + (betaa.*exey.*cle) //calcul de mu2

sqmusa = sqrt(mu./a)  //mu et a > 0 a ce stade des calculs

rnu2 = (-sle + (betaa.*ey.*ecceys))./(1 - ecceys) //remarque: (1-ecceys) different de 0 car e < 1 a ce stade des calculs
rki2 = (cle - (betaa.*ex.*ecceys))./(1 - ecceys)

//calcul des composantes des vecteurs position et vitesse
pos_car = (a.*.ones(3,1)) .* ( ( (rld2.*.ones(3,1)) .* x1 ) + ( (rmu2.*.ones(3,1)) .* y1 ) )
vel_car = (sqmusa.*.ones(3,1)) .* ( ( (rnu2.*.ones(3,1)) .* x1 ) + ( (rki2.*.ones(3,1)) .* y1 ) )


//CALCUL DU JACOBIEN DE LA TRANSFORMATION
if lhs==3
  //calcul de variables intermediaires

   //calcul des derivees de sqh par rapport a hx et hy
   xdsqhx = -hx ./ (4.0.*sqh)
   xdsqhy = -hy ./ (4.0.*sqh)

   //derivees des composantes de x1 et y1 par rapport a hx et hy
   xdx1hx = zeros(3,N)
   xdx1hx(1,:) = 0
   xdx1hx(2,:) = 0.5.* hy
   xdx1hx(3,:) = -hy .* xdsqhx
   xdy1hx = zeros(3,N)
   xdy1hx(1,:) = 0.5.* hy
   xdy1hx(2,:) = -hx
   xdy1hx(3,:) = sqh + (hx.*xdsqhx)
   xdx1hy = zeros(3,N)
   xdx1hy(1,:) = -hy
   xdx1hy(2,:) = 0.5.* hx
   xdx1hy(3,:) = -sqh - (hy.*xdsqhy)
   xdy1hy = zeros(3,N)
   xdy1hy(1,:) = 0.5.* hx
   xdy1hy(2,:) = 0
   xdy1hy(3,:) = hx .* xdsqhy

   //calcul derivee de beta par rapport a ex et ey
   xdbex = (ex ./ sq) .* (1 ./ (1+sq).^2)
   xdbey = (ey ./ sq) .* (1 ./ (1+sq).^2)

   //calcul des derivees de ld2 et mu2 par rapport a ex,ey et le
   xdldex = -1 - (ey2.*cle.*xdbex) + (betaa.*ey.*sle) + (xdbex.*exey.*sle)
   xdmuex = (betaa .* (eycexs - (ex.*sle))) + (xdbex.*ex.*eycexs)
   xdldey = - (betaa .* ((ey.*cle) + eycexs)) - (xdbey.*ey.*eycexs)
   xdmuey = -1 - (ex2.*sle.*xdbey) + (betaa.*ex.*cle) + (xdbey.*exey.*cle)
   xdldle = ((1 - (betaa.*ey2)) .* (-sle)) + (betaa.*exey.*cle)
   xdmule = ((1 - (betaa.*ex2)) .* cle) - (betaa.*exey.*sle)

   unme = 1-ecceys

   //calcul des derivees de nu2 et ki2 par rapport a ex,ey et le
   xdnuex =((-cle .* (sle - (betaa.*ey.*ecceys))) ./ unme.^2)+(((xdbex.*ey.*ecceys) + (betaa.*ey.*cle)) ./ unme)
   xdkiex =((cle .* (cle - (betaa.*ex.*ecceys))) / unme.^2)-(((xdbex.*ex.*ecceys) + (betaa .* (ecceys+(ex.*cle)))) ./unme)
   xdnuey =((-sle .* (sle - (betaa.*ey.*ecceys))) / unme.^2)+(((xdbey.*ey.*ecceys) + (betaa .* (ecceys+(ey.*sle)))) ./unme)
   xdkiey = ((sle .* (cle - (betaa.*ex.*ecceys))) / unme.^2)- (((xdbey.*ex.*ecceys) + (betaa.*ex.*sle)) ./ unme)
   xdnule = ((-cle + (betaa.*ey.*eycexs)) ./ unme)+ (((-sle + (betaa.*ey.*ecceys)) .* eycexs) ./ unme.^2)
   xdkile = ((-sle - (betaa.*ex.*eycexs)) ./ unme)+ (((cle - (betaa.*ex.*ecceys)) .* eycexs) ./ unme.^2)

   //calcul derivees de le par rapport a lm , ex et ey
   //    dle/dlm = 1 / 1-ecceys, dle/dex = sinue * dle/dlm
   //    dle/dey = -cosle * dle/dlm
   xdlelm = 1.0 ./ unme
   xdleex = sle .* xdlelm
   xdleey = -cle .* xdlelm

   //===========================================================
   // calcul des derivees partielles des vecteurs position et
   // vitesse par rapport a a,ex,ey,hx,hy,lm
   //===========================================================

      //par rapport a a
      j131 = ((rld2.*.ones(3,1)) .* x1) + ((rmu2.*.ones(3,1)) .* y1)
      j461 = ((-0.5.*sqmusa./a).*.ones(3,1)) .* (((rnu2.*.ones(3,1)) .* x1) + ((rki2.*.ones(3,1)) .* y1))

      //par rapport a ex ( = derivee/ex + derivee/le * derivee le/ex)
      j132 = (a.*.ones(3,1)) .* ( ...
                                    ( ((xdldex.*.ones(3,1)).*x1) + ((xdmuex.*.ones(3,1)).*y1) ) + ...
                                    ( ((xdldle.*.ones(3,1)).*x1) + ((xdmule.*.ones(3,1)).*y1) ) .* (xdleex.*.ones(3,1)) ...
                                    );
      j462 = (sqmusa.*.ones(3,1)) .* ( ...
                                         ( ((xdnuex.*.ones(3,1)).*x1) + ((xdkiex.*.ones(3,1)).*y1) ) + ...
                                         ( ((xdnule.*.ones(3,1)).*x1) + ((xdkile.*.ones(3,1)).*y1) ) .* (xdleex.*.ones(3,1)) ...
                                         );

      //par rapport a ey ( = derivee/ey + derivee/le * derivee le/ey)
      j133 = (a.*.ones(3,1)) .* ( ...
                                    ( ((xdldey.*.ones(3,1)).*x1) + ((xdmuey.*.ones(3,1)).*y1) ) + ...
                                    ( ((xdldle.*.ones(3,1)).*x1) + ((xdmule.*.ones(3,1)).*y1) ) .* (xdleey.*.ones(3,1)) ...
                                    );
      j463 = (sqmusa.*.ones(3,1)) .* ( ...
                                         ( ((xdnuey.*.ones(3,1)).*x1) + ((xdkiey.*.ones(3,1)).*y1) ) + ...
                                         ( ((xdnule.*.ones(3,1)).*x1) + ((xdkile.*.ones(3,1)).*y1) ) .* (xdleey.*.ones(3,1)) ...
                                         );

      //par rapport a hx
      j134 = (a.*.ones(3,1)) .* ( ((rld2.*.ones(3,1)) .* xdx1hx) + ((rmu2.*.ones(3,1)) .* xdy1hx) )
      j464 = (sqmusa.*.ones(3,1)) .* (((rnu2.*.ones(3,1)) .* xdx1hx) + ((rki2.*.ones(3,1)) .* xdy1hx))

      //par rapport a hy
      j135 = (a.*.ones(3,1)) .* ( ((rld2.*.ones(3,1)) .* xdx1hy) + ((rmu2.*.ones(3,1)) .* xdy1hy) )
      j465 = (sqmusa.*.ones(3,1)) .* (((rnu2.*.ones(3,1)) .* xdx1hy) + ((rki2.*.ones(3,1)) .* xdy1hy))

      //par rapport a um ( = derivee/ue * derivee de ue/um)
      j136 = (a.*.ones(3,1)) .* (((xdldle.*.ones(3,1)) .* x1) + ((xdmule.*.ones(3,1)) .* y1)) .* (xdlelm.*.ones(3,1))
      j466 = (sqmusa.*.ones(3,1)) .* (((xdnule.*.ones(3,1)) .* x1) + ((xdkile.*.ones(3,1)) .* y1)) .* (xdlelm.*.ones(3,1))

  //jacobien
  jacob = hypermat([6 6 N],[j131;j461;j132;j462;j133;j463;j134;j464;j135;j465;j136;j466] );
  if (N==1) then jacob=jacob(:,:,1); end

end

endfunction