//  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 [oe2,jacob] = CL_oe_convert(type_oe1,type_oe2,oe1, mu, cjac)
// Conversion of orbital elements
//
// Calling Sequence
// [oe2, jacob] = CL_oe_convert(type_oe1, type_oe2, oe1 [, mu, cjac])
//
// Description
// <itemizedlist><listitem>
// <p>Converts orbital elements from a given type to another.</p>
// <p>Available types of orbital elements are: "kep", "cir", "cireq", "equin" and "pv".</p>
// <p>Orbital elements of type "pv" are position and velocity in a 6xN vector ([pos;vel]). </p>
// <p>See <link linkend="Orbital elements">Orbital elements</link> for more details on orbital elements.</p>
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>
// <p>The gravitational parameter "mu" is only used when converting to or from "pv" type. </p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// type_oe1: (string) Type of input orbital elements ("kep", "cir", "cireq", "equin" or "pv") (1x1)
// type_oe2: (string) Type of output orbital elements ("kep", "cir", "cireq", "equin" or "pv") (1x1)
// oe1: Input orbital elements (6xN)
// oe2: Output orbital elements (6xN)
// mu: (optional) Gravitational constant [m^3/s^2]. Default value is %CL_mu
// cjac: (boolean, optional) Jacobian computation indicator: %t if jacobian should be computed. Default is %t
// jacob: (optional) Transformation jacobian (See <link linkend="Orbital elements">Orbital elements</link> for more details) (6x6xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Example 1
// pos = [7000.e3; 1000.e3; -500.e3];
// vel = [1.e3; 2.e3; 7e3];
// kep = CL_oe_convert("pv", "kep", [pos; vel])
//
// // Example 2
// cir = [7000.e3; 0.1; 0.2; 1; 2; 3];
// [kep, jacob1] = CL_oe_convert("cir", "kep", cir);
// [cir2, jacob2] = CL_oe_convert("kep", "cir", kep);
// cir2 - cir // => zero
// jacob2 * jacob1 // => identity


// Declarations:


// Existing orbital elements types
TYPES_OE = ["kep", "cir", "cireq", "equin", "pv"];

// Internal conversion function for simplification 
// - call existing functions (uses evstr to generate function from strings)
// - handles computation of jacobian 
// - handles arguments for "pv" type (2 args: pos and vel)
// NB: direct conversion function must exist  
function [oe2, jac] = oeconv_convert(type_oe1, type_oe2, oe1, mu, cjac)
  jac = []; 
  
  // generates function to be called
  // NB: "pv" replaced by "car"
  if (type_oe1 == "pv"); type_oe1 = "car"; end
  if (type_oe2 == "pv"); type_oe2 = "car"; end
  
  // function from/to cireq or equin are hidden functions (CL__oe...) 
  prefix  = "CL_oe_"; 
  if (or(type_oe1 == ["cireq", "equin"]) |  or(type_oe2 == ["cireq", "equin"]))
     prefix = "CL__oe_"
  end
  
  fun = evstr(prefix + type_oe1 + "2" + type_oe2); 
  
  if (type_oe1 == "car")
    // special case: 2 input args: pos and vel
    if (cjac)
      [oe2, jac] = fun(oe1(1:3,:), oe1(4:6,:), mu); 
    else
      [oe2] = fun(oe1(1:3,:), oe1(4:6,:), mu); 
    end

  elseif (type_oe2 == "car")
    // special case: 2 output args: pos and vel
    if (cjac)
      [pos, vel, jac] = fun(oe1, mu); 
    else
      [pos, vel] = fun(oe1, mu); 
    end
    oe2 = [pos; vel]; 

  else
    // Note: mu not passed to function
    if (cjac)
      [oe2, jac] = fun(oe1);  
    else
      [oe2] = fun(oe1); 
    end
  end
  
endfunction


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

if (argn(2) < 3)
  CL__error("Invalid number of input arguments");
end

// check orbital elements types   
if ~(or(type_oe1 == TYPES_OE) & or(type_oe2 == TYPES_OE))
  CL__error("Invalid type of orbital elements");
end

if (size(oe1,1) <> 6 & size(oe1,1) <> 0)
  CL__error("Invalid orbital elements (number of rows)");
end

// cjac = %f if no corresponding output argument exists (optimization)
if (argn(1) == 1); cjac = %f; end

// Special cases
if (oe1 == [])
  oe2 = []; 
  jacob = []; 

elseif (type_oe1 == type_oe2)
  oe2 = oe1;
  jacob = []; 
  if (cjac)
    N = size(oe1, 2);
    jacob = matrix(repmat(eye(6,6),1,N),6,6,N); // identity hypermatrix
  end

else
  // conversion 
  [oe2, jacob] = oeconv_convert(type_oe1, type_oe2, oe1, mu, cjac); 

end  

endfunction