//  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 [mean_t2,osc_t2] = CL_ex_lyddaneMan(t1,mean_t1,t2, tman,dvman,er,mu,j1jn)
// Lyddane orbit propagation analytical model with maneuvers - DEPRECATED
//
// Calling Sequence
// [mean_t2,osc_t2] = CL_ex_lyddaneMan(t1,mean_t1,t2 [,tman,dvman,er,mu,j1jn])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p>Replacement function: <link linkend="CL_ex_propagateMan">CL_ex_propagateMan</link></p>
// <p></p>
// </listitem>
// <listitem>
// <p>Computes the mean and osculating orbital elements <b>mean_t2</b> 
// and <b>osc_t2</b> at time <b>t2</b> given the mean elements 
// <b>mean_t1</b> at time <b>t1</b>. </p>
// <p> Maneuvers (i.e. velocity increments defined by <b>tman,dvman</b>)
// are taken into account. </p>
// <p><b>tman</b> defines the time of the maneuver.</p>
// <p><b>dvman</b> = [delta-Vx;delta-Vy;delta-Vz] defines the delta-V 
// in the QSW local frame at the time of maneuver.</p>
// <p>If one maneuver occurs at one of the final times, the output orbital elements mean_t2 and osc_t2 
// include the effect of the maneuver at that time.</p>
// <p>The final times as well as the maneuver time should be sorted (increasing order).</p>
// <p>Zonals coefficients up to J5 are taken into account.</p>
// <p></p></listitem>
// <listitem> 
// <p>The type of orbital elements is keplerian.</p> 
// <p></p></listitem>
// <listitem>
// <p>Warnings :</p>
// <p> - This function does not work for inclinations close to the critical inclinations (63.43494882 deg and 116.5650512 deg)</p>
// <p> - This function does not work for eccentricities higher than 0.9</p>
// <p></p></listitem>
// <listitem>
// <p>Warning :</p>
// <p> - The input argument "zonals" is deprecated as of CelestLab v3.0.0. It has been replaced by "j1jn".</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// t1: Initial time [days] (1x1)
// mean_t1: Mean Keplerian elements at time t1 [sma;ecc;inc;pom;raan;anm] (6x1)
// t2: Final time [days] (1xN)
// tman: (optional) Maneuver times [days] (1xP)
// dvman: (optional) Delta-V in QSW local frame [dvx;dvy;dvz] [m/s] (3xP)
// er: (optional) Equatorial radius [m] (default is %CL_eqRad)
// mu: (optional) Gravitational constant [m^3/s^2] (default value is %CL_mu)
// j1jn: (optional) Vector of zonal coefficients J1 to Jn (troncated to J5) to be used (default is %CL_j1jn(1:5)) (1xM)
// mean_t2: Mean Keplerian elements at t2 [sma;ecc;inc;pom;raan;anm] (6 x N)
// osc_t2: Osculating Keplerian elements at t2 [sma;ecc;inc;pom;raan;anm] (6 x N)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_ex_lyddane
// CL_ex_meanLyddane
// CL_fr_qswMat
//
// Examples
// t1 = 21915;
// mean_t1 = [42166712 ; 7.9e-3 ; CL_deg2rad([97.2; 89; 125; 0])];
// t2 = [21916 21918 21920];
// tman = [21915.5 21916.5];
// dvman = [ [ 1 ; 1 ; 1] , [ 1 ; 0 ; 1]];
// [mean_t2,osc_t2] = CL_ex_lyddaneMan(t1,mean_t1,t2,tman,dvman,er=6378138)
//

// Declarations:

// Code:
CL__warnDeprecated(); // deprecated function

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", 1:5); end

// Ensure that j1jn has at least 5 elements (fill with zeros)
j1jn = [matrix(j1jn, 1, -1), zeros(1,5)];

// local function : propagates to t and apply dv at t
// NB : use local variables err, mu, j1jn

function [mean_t] = extrap_lyd_man(t1, mean_t1, t, dv)
  [mean_t,osc_t] = CL_ex_lyddane(t1,mean_t1,t,er,mu,j1jn);
  [osc_t_car_pos,osc_t_car_vel] = CL_oe_kep2car(osc_t,mu);
  dv_inert = CL_fr_qsw2inertial(osc_t_car_pos,osc_t_car_vel,dv);
  osc_t_car_vel = osc_t_car_vel + dv_inert;
  osc_t = CL_oe_car2kep(osc_t_car_pos,osc_t_car_vel,mu);
  mean_t = CL_ex_meanLyddane(osc_t,er,mu,j1jn);
endfunction


Nt1 = size(t1,2)            // should be 1
Nmean_t1 = size(mean_t1,2)  // should be 1
if (Nt1 ~= 1) then
   CL__error('t1 should be (1x1)');
end
if (Nmean_t1 ~= 1) then
   CL__error('mean_t1 should be (6x1)');
end

Nman = 0; 
Ndvman = 0;
if exists('tman', 'local'); Nman = size(tman,2); end      // number of maneuver times 
if exists('dvman', 'local'); Ndvman = size(dvman,2), end  // number of maneuver DV 

Nt2 = size(t2,2);  //extrapolation dates

if (Nman ~= Ndvman ) then
   CL__error('tman et dvman not of the same size. Should be (1xP) and (3xP)');
end


// Cas of no maneuvers :
if (Nman == 0) then 
   [mean_t2,osc_t2] = CL_ex_lyddane(t1,mean_t1,t2,er,mu,j1jn);
   return; // *** EXIT ***  
end

// initialization of output arguments
mean_t2 = zeros(6,Nt2);   
osc_t2 = zeros(6,Nt2);


// Propagatiopn to first maneuver (outputs for t < first man)
I = find(t2 < tman(1));
if (I ~= [])
   [mean_t2(:,I),osc_t2(:,I)] = CL_ex_lyddane(t1,mean_t1,t2(I),er,mu,j1jn);
end

// state after 1st maneuver
mean_tman = extrap_lyd_man(t1, mean_t1, tman(1), dvman(:,1)); 


// Extrapolation entre la kieme-1 et la kieme maneuvre :
for kman = 2 : Nman

   I = find (t2 >= tman(kman-1) & t2 < tman(kman));
  
   if (I ~= [])
      [mean_t2(:,I),osc_t2(:,I)] = CL_ex_lyddane(tman(kman-1),mean_tman,t2(I),er,mu,j1jn);
   end
  
   // state after kman_th maneuver:
   mean_tman = extrap_lyd_man(tman(kman-1), mean_tman, tman(kman), dvman(:,kman-1)); 

end

// Extrapolation jusqu'a la derniere date :
I = find(t2 >= tman($));
if (I ~= [])
   [mean_t2(:,I),osc_t2(:,I)] = CL_ex_lyddane(tman($),mean_tman,t2(I),er,mu,j1jn);
end

endfunction