//  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_extrap(type_oe, cjd0, mean_oe0, cjd, params, res, frame, ut1_tref, tt_tref)
// Orbit propagation using STELA
//
// Calling Sequence
// [result1, result2, ...] = CL_stela_extrap(type_oe, cjd0, mean_oe0, cjd, params, res [, frame, ut1_tref, tt_tref])
//
// Description
// <itemizedlist>
// <listitem>
// <p>Propagates a mean orbital state using STELA. </p>
// <p></p>
// <p>The results are described by <b>res</b>: </p>
// <p>- "m": Mean orbital elements (6xN) </p>
// <p>- "tm": Transition matrix for the mean orbital elements (6x6xN) </p>
// <p>- "i": Information data about the propagation, including an XML description of the inputs for 
// the STELA tool  </p>
// <p></p>
// </listitem>
// <listitem>
// <p><b>Notes</b>:</p>
// <p>- The dates offsets relative to the initial date (that is: cjd-cjd0) must be proportional to the 
// integration time step. </p>
// <p>- The integration time step is typically chosen equal to 1 day. </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 back to the specified type and frame. </p>
// <p>- The time scale offsets (ut1_tref and tt_tref) should be constant values (size = 1x1). </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". 
// cjd0: Initial date (CJD, time scale: TREF) (1x1)
// mean_oe0: Initial mean orbital elements (6x1) 
// cjd: Final dates (CJD, time scale: TREF) (1xN)
// params: (structure) Propagation model parameters. 
// res: (string) Wanted results: "m", "tm", "i" (1xP).  
// 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, result2, ...: Mean orbital elements, transition matrix, or information data (see above)
//
// Authors
// CNES - DCT/SB/MS
//
// Examples
//
// // Initial date (cjd, time scale: TREF)
// cjd0 = 20000; 
// 
// // Keplerian mean orbital elements (frame: ECI)
// mean_kep0 = [7.e6; 1.e-3; 98*(%pi/180); %pi/2; 0; 0]; 
//
// // Final dates
// cjd = cjd0 + (0:60); 
//
// // STELA model parameters (default values)
// params = CL_stela_params(); 
//
// // Propagation
// [mean_kep, info] = CL_stela_extrap("kep", cjd0, mean_kep0, cjd, params, ["m", "i"]);
//
// // Plot inclination (deg)
// scf(); 
// plot(cjd, mean_kep(3,:) * (180/%pi));


// Declarations:

// Code
if (~exists("res", "local")); res=["m", "tm", "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", "tm", "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)
// NB: To simplify the code, the jacobians are always computed. 


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

convfr = %f; 
if (~(frame == "CIRS" | (frame == "ECI" & CL_configGet("ECI_FRAME") == "CIRS")))
  // if frame <> CIRS => compute frame conversion for all dates
  [M, omega] = CL_fr_convertMat(frame, "CIRS", [cjd0, cjd], ut1_tref=ut1_tref, tt_tref=tt_tref); 
  convfr = %t; 
end

if (convfr)
  [pv0, J1] = CL_oe_convert(type_oe, "pv", mean_oe0, cst.mu); 
  [p0, v0, J2] = CL_rot_pvConvert(pv0(1:3), pv0(4:6), M(:,:,1), omega(:,1)); 
  [x0, J3] = CL_oe_convert("pv", "cireq", [p0; v0], cst.mu); 
  jac1 = J3*J2*J1; 

else
  [x0, jac1] = CL_oe_convert(type_oe, "cireq", mean_oe0, cst.mu); 
end

// if [] => remains []
delta_t = ((cjd - cjd0) .* ones(cjd)) * 86400; 

// compute transition matrix only if necessary
// NB: xosc is not used (and not computed) ! 
sres = 0; 
if (find(res == "tm") <> []); sres = 1; end
[x, xosc, dxd0, info] = CLx_stela_extrap(cjd0, x0, delta_t, params, res=sres, 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.tm = []; 
result.m = []; 


if (convfr)
  // 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(:,:,2:$), omega(:,2:$)); 
  [mean_oe, J3] = CL_oe_convert("pv", type_oe, [p;v], cst.mu); 
  jac2 = J3*J2*J1; 

else
  [mean_oe, jac2] = CL_oe_convert("cireq", type_oe, x, cst.mu); 
      
end


result.m = mean_oe; 

N = size(cjd, 2); 
jac1 = hypermat([6,6,N], repmat(jac1,1, N)); 
if (dxd0 <> []) 
  result.tm = jac2 * dxd0(1:6,1:6,1:N) * jac1;     
end

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

      
endfunction