//  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_ex_propagate(mod, type_oe, t1, mean_oe_t1, t2, res, er, mu, j1jn)
// Orbit propagation (all analytical models)
//
// Calling Sequence
// [result1, result2] = CL_ex_propagate(mod, type_oe, t1, mean_oe_t1, t2, res [, er, mu, j1jn])
//
// Description
// <itemizedlist><listitem>
// <p>Propagates orbital elements using one analytical model. </p>
// <p></p>
// <p>The available propagation models are: </p>
// <p>"central": Central force (osculating elements = mean elements) </p>
// <p>"j2sec": Secular effects of J2 (osculating elements = mean elements by convention)</p>
// <p>"lydsec": Lyddane (mean elements include secular effects only) </p>
// <p>"lydlp": Lyddane (mean elements include secular and long period effects) </p>
// <p>"eckhech": Eckstein-Hechler  (mean elements include secular and long period effects) </p>
// <p></p></listitem>
// <listitem>
// <p>Notes:</p>
// <p>- There can be 1 or N initial times, and 1 or N final times. </p>
// <p>- Conversions take place if the type of orbital elements is not the "natural" type 
// for the model. There is no conversion for the "central" and "j2sec" models. </p>
// <p></p></listitem>
// <listitem>
// <p>See <link linkend="Propagation models">Propagation models</link> for more details.</p>
// <p></p></listitem>
// </itemizedlist>
//
// Parameters
// mod: (string) Model name: "central", "j2sec", "lydsec", "lydlp", "eckhech". (1x1)
// type_oe: (string) Type of orbital elements used for input/output: "kep", "cir", "cireq" or "equin" (1x1)
// t1: Initial time [days] (1x1 or 1xN)
// mean_oe_t1: Mean orbital elements at time t1 (6x1 or 6xN) 
// t2: Final time [days] (1xN or 1x1)
// res: (string) Type of output (mean or osculating): "m", "o", "mo", "om"  (1x1)
// er: (optional) Equatorial radius [m]. Default is %CL_eqRad
// mu: (optional) Gravitational constant [m^3/s^2]. Default is %CL_mu
// j1jn: (optional) Vector of zonal harmonics. Default is %CL_j1jn (Nz x 1)
// result1, result2: Mean or osculating orbital elements at t2 (6xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_ex_osc2mean
// CL_ex_mean2osc
//
// Examples
// mean_kep0 = [7.e6; 1.e-3; 1; %pi/2; 0.1; 0.2]; 
// t0 = 0; 
// t = 1:4; 
// osc_kep = CL_ex_propagate("eckhech","kep",t0,mean_kep0,t,"o") 
//

// Declarations:

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

Models = [ "central", "j2sec", "lydsec", "lydlp", "eckhech" ];  
// "natural" types for each model
Types_oe_nat = [ "-", "-", "kep", "kep", "cir" ]; // "-" => any
Res = ["mo", "om", "o", "m"]; 

lhs = argn(1); 
if (lhs > 2)
  CL__error("Wrong number of output arguments"); 
end

rhs = argn(2); 
if (rhs < 6)
  CL__error("Wrong number of input arguments"); 
end

I = find(res == Res); 
if (I == [])
  CL__error("Invalid value for argument ''res''"); 
end

imod = find(mod == Models); 
if (imod == [])
  CL__error("Invalid model name"); 
end

compute_osc = %f; 
if (res == "om" | res == "o"); compute_osc = %t; end
if (res == "mo" & lhs > 1); compute_osc = %t; end

if (~or(type_oe == ["kep", "cir", "cireq", "equin"]))
  CL__error("Invalid type of orbital elements"); 
end


// -------------------------------
// Main
// -------------------------------

nat_type = Types_oe_nat(imod); 
convert = (type_oe <> nat_type & nat_type <> "-"); 

mean_oe_t2 = []; 
osc_oe_t2 = []; 

if (convert)
  // converts to "natural" type 
  mean_oe_t1 = CL_oe_convert(type_oe, nat_type, mean_oe_t1, mu); 
end

if (imod == 1) 
  // central force
  mean_oe_t2 = CL_ex_kepler(t1, mean_oe_t1, t2, mu); 
  if (compute_osc); osc_oe_t2 = mean_oe_t2; end 

elseif (imod == 2) 
  // Secular J2 (osc = mean by convention) 
  // Note: use a specific function that propagates the orbital elements directly 
  // (no conversion to "kep" type) 
  j1jn = [matrix(j1jn,-1,1); 0; 0]; 
  mean_oe_t2 = CL__ex_propagJ2sec(type_oe, t1, mean_oe_t1, t2, er, mu, j1jn(2)); 
  if (compute_osc); osc_oe_t2 = mean_oe_t2; end  

elseif (imod == 3) 
  // Lyddane (mean = secular)
  [mean_oe_t2, osc_oe_t2] = CL__ex_propag_lydsec(t1, mean_oe_t1, t2, er, mu, j1jn, compute_osc); 

elseif (imod == 4) 
  // Lyddane (mean = secular + long periods)
  [mean_oe_t2, osc_oe_t2] = CL__ex_propag_lydlp(t1, mean_oe_t1, t2, er, mu, j1jn, compute_osc); 

elseif (imod == 5) 
  // Eckstein-Hechler
  [mean_oe_t2, osc_oe_t2] = CL__ex_propag_eckhech(t1, mean_oe_t1, t2, er, mu, j1jn, compute_osc); 

end

if (convert)
  // converts back to "initial" type 
  mean_oe_t2 = CL_oe_convert(nat_type, type_oe, mean_oe_t2, mu); 
  if (compute_osc)
    osc_oe_t2 = CL_oe_convert(nat_type, type_oe, osc_oe_t2, mu); 
  end
end


// generates outputs
for i = 1 : lhs
  if (part(res, i) == "m")
    varargout(i) = mean_oe_t2;
  elseif (part(res, i) == "o")
    varargout(i) = osc_oe_t2;
  else
    varargout(i) = []; 
  end
end


endfunction