//  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 [kep,jacob] = CL_oe_cirEqua2kep(cir_equa)
// Circular equatorial adapted to keplerian orbital elements - DEPRECATED
//
// Calling Sequence
// [kep,jacob] = CL_oe_cirEqua2kep(cir_equa)
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p>Replacement function: <link linkend="CL__oe_cireq2kep">CL__oe_cireq2kep</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 
// classical keplerian orbital elements.</p>
// <p>The transformation jacobian is 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)
// kep: classical keplerian orbital elements [sma;ecc;inc;pom;raan;anm] [m,rad] (6xN)
// 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_kep)
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Example 1
// cirequa = [42166.712;-7.900e-6;1.100e-4;1.200e-4;-1.16e-4;5.300];
// [kep,jacob1]=CL_oe_cirEqua2kep(cirequa);
// [cirequa2,jacob2]=CL_oe_kep2cirEqua(kep);


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

// Code:

CL__warnDeprecated(); // deprecated function

[lhs,rhs] = argn();

if lhs==2
  compute_jacob=%t
else
  compute_jacob=%f
end

N = size(cir_equa,2);

// valeurs de test
eps_cir = EPS_ORB.cir
eps_equa = EPS_ORB.equa

// calcul de a
kep=[]
kep(1,:) = cir_equa(1,:)

// calcul de i
norme_carree = cir_equa(4,:) .* cir_equa(4,:) + cir_equa(5,:) .* cir_equa(5,:)
ii1 = find(norme_carree > 4)
if ii1~=[] then CL__error('inclination vector norm too big'); end

deux_sinus_i_sur_deux = sqrt(norme_carree)
inclinaison = 2.0 .* asin(deux_sinus_i_sur_deux./2) // resultat dans [0,+pi] car sin(i/2) > 0
kep(3,:) = inclinaison

// calcul de grand omega
jj1 = find(~(sin(inclinaison) < eps_equa))
grand_omega = zeros(kep(1,:)) // choix arbitraire: mise a 0
grand_omega(jj1) = CL__sc2angle(cir_equa(4,jj1),cir_equa(5,jj1))
kep(5,:)= grand_omega

// calcul de e
e2 = cir_equa(2,:).^2 + cir_equa(3,:).^2
eccentricite = sqrt(e2)
kep(2,:) = eccentricite

// calcul de petit omega + grand omega
jj2 = find(~(eccentricite < eps_cir)) //orbite circulaire: petit omega + grand omega indetermine
pom_plus_gom = zeros(kep(1,:))   // choix arbitraire: mise a 0
pom_plus_gom(jj2) = CL__sc2angle(cir_equa(2,jj2),cir_equa(3,jj2))

// calcul de M
kep(6,:) = cir_equa(6,:) - pom_plus_gom

// calcul de petit omega
kep(4,:) = pom_plus_gom - grand_omega

//reduce angles to 0 -> 2pi
kep(3:6,:) = CL_rMod(kep(3:6,:),2*%pi)

//jacobien calcul
if compute_jacob

   jj3 = find( ~((sin(inclinaison)<eps_equa) | (eccentricite<eps_cir)) )
   jj4 = find( (sin(inclinaison)<eps_equa) | (eccentricite<eps_cir) )

   if jj4~=[] then CL__warning("Non computable jacobian for sin(i)~0 or e~0. Zero matrix given"); end


   un_sur_e = 1 ./ eccentricite(jj3)
   un_sur_nc = 1 ./ norme_carree(jj3)
   exsure2 = cir_equa(2,jj3) ./ e2(jj3)
   eysure2 = cir_equa(3,jj3) ./ e2(jj3)
   jacob = hypermat([6 6 N]) //jacob corresponding to jj4 indexs are zero matrix
   jacob(1,1,jj3) = 1
   jacob(2,2,jj3) = cir_equa(2,jj3) .* un_sur_e
   jacob(2,3,jj3) = cir_equa(3,jj3) .* un_sur_e
   jacob(3,4,jj3) = cir_equa(4,jj3) ./ sin(inclinaison(jj3))
   jacob(3,5,jj3) = cir_equa(5,jj3) ./ sin(inclinaison(jj3))
   jacob(4,2,jj3) = - eysure2
   jacob(4,3,jj3) =   exsure2
   jacob(4,4,jj3) = cir_equa(5,jj3) .* un_sur_nc
   jacob(4,5,jj3) = - cir_equa(4,jj3) .* un_sur_nc
   jacob(5,4,jj3) = - cir_equa(5,jj3) .* un_sur_nc
   jacob(5,5,jj3) = cir_equa(4,jj3) .* un_sur_nc
   jacob(6,2,jj3) =   eysure2
   jacob(6,3,jj3) = - exsure2
   jacob(6,6,jj3) = 1

   if (N==1) then jacob=jacob(:,:,1); end

end

endfunction