//  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 [varargout] = CL_stela_deriv(type_oe, cjd, mean_oe, params, res, frame, ut1_tref, tt_tref)
// Orbital elements time derivatives using STELA
//
// Calling Sequence
// [result1, ...] = CL_stela_deriv(type_oe, cjd, mean_oe, params, res [, frame, ut1_tref, tt_tref])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Computes the orbital time derivatives due to various perturbations using STELA. </p>
// <p></p>
// <p>The results are described by <b>res</b>: </p>
// <p>- "m": Derivatives of mean orbital elements (6xN) </p>
// <p>- "i": Information data (structure) </p>
// <p></p>
// <p>The input orbital elements are internally converted to the 
// adequate type and frame before the call to STELA. The results are then converted to 
// the same type and frame. </p>
// <p></p>
// </listitem>
// <listitem>
// <p><b>Warning</b>:</p>
// <p>The dynamic motion (that is the apparent acceleration) of the chosen reference frame with respect to 
// CIRS (which is the frame in which STELA integrates the motion) is not included in the derivatives. 
// The results must then be used with caution if the chosen frame is not ECI. This may change in future 
// versions. </p>
// <p></p>
// </listitem>
// <listitem>
// <p><b>Notes</b>:</p>
// <p>- The integration time step parameter (in the "params" structure) is not used but must exist. </p>
// <p></p>
// </listitem>
// <listitem>
// <p>See the <link linkend="STELA">STELA</link> page for more details.</p>
// <p></p>
// </listitem>
// </itemizedlist>
// 
// Parameters
// type_oe: (string) Type of orbital elements: "kep", "cir", "cireq", "equin", "pv". 
// cjd: Reference dates (CJD, time scale: TREF) (1xN)
// mean_oe: Mean orbital elements for the reference orbit (6xN) 
// params: (structure) Propagation model parameters. 
// res: (string) Wanted results: "m", "i".  
// frame: (string, optional) Input/output frame. Default is "ECI" 
// ut1_tref: (optional) UT1-TREF [seconds]. Default is %CL_UT1_TREF (1x1)
// tt_tref: (optional) TT-TREF [seconds]. Default is %CL_TT_TREF (1x1)
// result1, ...: Time derivatives or information data
//
// Authors
// CNES - DCT/SB/MS
//
// Examples
// // Generate reference orbit (Sun-synchronous, MLTAN=14h, frame: ECI)
// kep0 = [7.e6; 1.136e-3; 1.7085241; %pi/2; 0; 0];
// cjd0 = 20000;
// cjd = cjd0 + (0:365); // days, time scale: TREF
// kep = CL_ex_propagate("eckhech", "kep", cjd0, kep0, cjd, "m"); 
//
// // STELA model parameters (default values)
// params = CL_stela_params(); 
//
// // Time derivatives
// [dkepdt, info] = CL_stela_deriv("kep", cjd, kep, params, ["m", "i"]);
//
// // Plot inclination derivative (deg/year)
// scf(); 
// plot(cjd, dkepdt(3,:) * (180/%pi) * 86400 * 365.25);


// Declarations:

// Code
if (~exists("res", "local")); res=["m", "i"]; end
if (~exists("frame", "local")); frame="ECI"; end
if (~exists("ut1_tref", "local")); ut1_tref = CL__dataGetEnv("UT1_TREF"); end 
if (~exists("tt_tref", "local")); tt_tref = CL__dataGetEnv("TT_TREF"); end 

// check extension module
CL__checkExtensionModule(); 

// check res
if (typeof(res) <> "string" | size(res,1) <> 1)
  CL__error("Invalid type or size for res"); 
end

RES_REF = ["m", "i"]; 
if (setdiff(res, RES_REF) <> [])
  CL__error("Invalid values for res"); 
end

if (typeof(frame) <> "string" | size(frame, "*") <> 1)
  CL__error("Invalid argument type or size for frame"); 
end

if (size(cjd,1) > 1)
  CL__error("Invalid size for argument cjd"); 
end

// converts orbital elements to "cireq" type
// NB: uses "mu" from STELA
// converts frame to "CIRS" if necessary ("CIRS": expected by Stela)

// STELA constants (mu)
CL__stela_manageConst(); 
cst = CLx_stela_getConst(1); 

if (~(frame == "CIRS" | (frame == "ECI" & CL_configGet("ECI_FRAME") == "CIRS")))
  [M, omega] = CL_fr_convertMat(frame, "CIRS", cjd, ut1_tref=ut1_tref, tt_tref=tt_tref); 

  [pv] = CL_oe_convert(type_oe, "pv", mean_oe, cst.mu); 
  [p, v] = CL_rot_pvConvert(pv(1:3,:), pv(4:6,:), M, omega); 
  [x] = CL_oe_convert("pv", "cireq", [p; v], cst.mu); 
  
  // for the jacobian: need the inverse transformation: dx = J*dX => dX/dt = J-1 * dx/dt
  // => recompute! 
  // change direction
  [M, omega] = CL_rot_compose(M, omega, -1); 

  [pv, J1] = CL_oe_convert("cireq", "pv", x, cst.mu); 
  [p, v, J2] = CL_rot_pvConvert(pv(1:3,:), pv(4:6,:), M, omega); 
  [X, J3] = CL_oe_convert("pv", type_oe, [p;v], cst.mu); // X = mean_oe

  jact = J3*J2*J1; 

else
  // use "pv" as intermediate type (as direct transformation may not exist) 
  // for the jacobian: need the inverse transformation: dx = J*dX => dX/dt = J-1 * dx/dt
   x = CL_oe_convert(type_oe, "cireq", mean_oe, cst.mu);  
  [X, jact] = CL_oe_convert("cireq", type_oe, x, cst.mu); // X = mean_oe, unused

end

// NB: res: unused
[xdot, info] = CLx_stela_deriv(cjd, x, params, res=1, ut1_tref, tt_tref); 

// check status
if (info.status < 0)
  CL__error(info.err_msg); 
end

// converts orbital elements 
// converts frame if necessary

result = struct(); 
result.i = info; 
result.m = []; 
result.m = jact * xdot; 


for k = 1 : size(res,2)
  varargout(k) = result(res(k)); 
end


endfunction