//  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 kepculated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [kep,jacob] = CL_oe_car2kep(pos,vel, mu)
// Cartesian to classical Keplerian orbital elements
//
// Calling Sequence
// [kep, jacob] = CL_oe_car2kep(pos, vel [, mu])
//
// Description
// <itemizedlist><listitem>
// <p>Converts cartesian orbital elements 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>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// pos: Position vector [x;y;z] [m] (3xN)
// vel: Velocity vector [vx;vy;vz] [m/s] (3xN)
// mu : (optional) Gravitational constant. [m^3/s^2] (default value is %CL_mu)
// kep: Classical Keplerian orbital elements [sma;e;inc;pom;raan;M] [m,rad] (6xN)
// jacob: (optional) Transformation jacobian (See <link linkend="Orbital elements">Orbital elements</link> for more details) (6x6xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_oe_kep2car
//
// Examples
// // Example 1
// pos = [7000.e3; 1000.e3; -500.e3];
// vel = [1.e3; 2.e3; 7e3];
// kep = CL_oe_car2kep(pos, vel);
//
// // Example 2
// pos = [7000.e3; 1000.e3; -500.e3];
// vel = [1.e3; 2.e3; 7e3];
// [kep, jacob1] = CL_oe_car2kep(pos, vel);
// [pos2, vel2, jacob2] = CL_oe_kep2car(kep); 
// pos2 - pos // => zero
// vel2 - vel // => zero
// jacob2 * jacob1 // => identity

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


// Code:
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

// Handle [] cases
if (pos == [] | vel == [])
  kep = []
  jacob = [];
  return;
end

// Check validity of inputs
[isvalid,orbit_type] = CL__oe_isValid("car",[pos;vel],mu);
if (~isvalid); CL__error("Invalid position velocity (zero norm or colinear)"); end;
if (find(orbit_type == 3) <> []);
  CL__error("Invalid position velocity (parabolic orbit)"); 
end


// Flag to compute jacobian
cjac = %f;
if (argn(1) == 2)
  cjac = %t;
end

// Check input sizes
[pos,vel,N] = CL__checkInputs(pos,3, vel,3);

kep = zeros(6,N);

jacob = []; 
if (cjac)
  jacob = zeros(6,6,N); 
end

// Elliptic orbits
I = find(orbit_type == 1);
if (I <> [])
  [kep(:,I), jacob_tmp] = CL__oe_car2kep_ell(pos(:,I),vel(:,I),mu,cjac);
  if (cjac)
    jacob(:,:,I) = jacob_tmp; 
  end
end

// Hyperbolic orbits
I = find(orbit_type == 2);
if (I <> [])
  [kep(:,I), jacob_tmp] = CL__oe_car2kep_hyp(pos(:,I),vel(:,I),mu,cjac);
  if (cjac)
    jacob(:,:,I) = jacob_tmp; 
  end
end

// Parabolic orbits
// I = find(orbit_type == 3);
// if (I <> [])
//   [kep(:,I), jacob_tmp] = CL__oe_car2kep_parab(pos(:,I),vel(:,I),mu,cjac);
//   if (cjac)
//     jacob(:,:,I) = jacob_tmp; 
//   end
// end

endfunction