//  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 [t,rel_pos_vel] = CL_cw_contPropa(t_ini,rel_pos_vel_ini,gama_diff,alt,t_end, delta_t,ballistic_coef_chaser,ballistic_coef_target,er,mu)
// Continuous thrust trajectory - DEPRECATED
//
// Calling Sequence
// [t,rel_pos_vel] = CL_cw_contPropa(t_ini,rel_pos_vel_ini,gama_diff,...
//    alt,t_end,delta_t,[ballistic_coef_chaser,ballistic_coef_target,er,mu])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the position and velocity vector of a chaser relative to a target when 
// the motion of the chaser is affected by continuous thrust. </p>
// <p>The output argument <b>rel_pos_vel</b> contains
// the positions and velocities of the chaser relative to the target in the target's LVLH 
// reference frame at the times <b>t</b>. </p>
// <p>If no ballistic coefficients are given, acceleration due to drag is not taken into account. </p>
// <p><inlinemediaobject><imageobject><imagedata fileref="contPropa.gif"/></imageobject></inlinemediaobject></p>
// </listitem>
// </itemizedlist>
//
// Parameters
// t_ini: Initial time [seconds] (1x1)
// rel_pos_vel_ini: Initial position and velocity vector of the chaser relative to the target in the target's LVLH frame [X;Y;Z;Vx;Vy;Vz] [m;m/s] (6x1)
// gamma_diff: Continuous thrust differential acceleration [ax, ay, az] [m/s**2] (3xM)
// alt: Altitude of target [m]
// t_end: Final time [seconds] (1x1)
// delta_t: (optional) Time step [s] (default is 100)
// ballistic_coef_chaser: (optional) Ballistic coefficient of chaser : S*cx/m S=equivalent drag surface, cx = drag coefficient, m = mass [m^2/kg] (default is 0)
// ballistic_coef_target: (optional) Ballistic coefficient of target : S*cx/m S=equivalent drag surface, cx = drag coefficient, m = mass [m^2/kg] (default is 0)
// er: (optional) Equatorial radius [m] (default is %CL_eqRad)
// mu: (optional) Gravitational constant [m^3/s^2] (default value is %CL_mu)
// t: Computation times [seconds] (1xN)
// rel_pos_vel: Position and velocity vectors of chaser in target's LVLH frame at corresponding times [X;Y;Z;Vx;Vy;Vz] (6xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_cw_impulPropa
// CL_cw_diffDrag
// CL_cw_Mmatrix
// CL_cw_Nmatrix
//
// Bibliography
// 1) Mecanique spatiale, CNES - Cepadues 1995, Tome II
//
// Examples
// t_ini = 11550.677;
// rel_pos_vel_ini = [-250; 0. ; 0 ; -3e-6 ; 0. ; 3e-6];
// gama_diff = [0. ; 0. ; 0.0003238];
// alt = 450.e3;
// ballistic_coef_chaser = 100 * 1 / 20.e3;
// ballistic_coef_target = 200 * 1 / 400.e3;
// t_end = 4*3600; // 4 hours of ballistic propagation
// delta_t = 100; // 100 seconds time step
// [t,rel_pos_vel] = CL_cw_contPropa(t_ini, ..
// rel_pos_vel_ini,gama_diff,alt,t_end,delta_t, ..
// ballistic_coef_chaser,ballistic_coef_target)
//

// Declarations:


// Code:

CL__warnDeprecated(); // deprecated function

if (~exists('ballistic_coef_chaser','local')); ballistic_coef_chaser=0; end
if (~exists('ballistic_coef_target','local')); ballistic_coef_target=0; end
if (~exists('delta_t','local')); delta_t=100; end
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

nbpas = int( (t_end-t_ini)/delta_t );     //number of steps
pasfinal = t_end - nbpas*delta_t - t_ini; //last step to perform. non zero if nbpas was not int.

acc_diff = CL_cw_diffDrag(alt,ballistic_coef_chaser,ballistic_coef_target,er,mu)  //differential acceleration between target and chaser

gama = [acc_diff;0;0];   //differential acceleration vector
gama = gama + gama_diff;  //add differential acceleration to continuous thrust acceleration vector.

mat_M = CL_cw_Mmatrix(alt,delta_t,er,mu);
mat_N = CL_cw_Nmatrix(alt,delta_t,er,mu);

t = zeros(1,nbpas+2);
rel_pos_vel = zeros(6,nbpas+2);
t(1) = t_ini;
rel_pos_vel(:,1) = rel_pos_vel_ini;

for i=1:nbpas
  V = rel_pos_vel(:,i);
  W = mat_M*V + mat_N*gama;
  rel_pos_vel(:,i+1) = W;
  t(i+1) = t(i) + delta_t
end

//last step
if pasfinal~=0
  mat_M = CL_cw_Mmatrix(alt,pasfinal,er,mu);
  mat_N = CL_cw_Nmatrix(alt,pasfinal,er,mu);

  V = rel_pos_vel(:,nbpas+1);
  W = mat_M*V + mat_N*gama;
  rel_pos_vel(:,nbpas+2) = W;
  t(nbpas+2) = t(nbpas+1) + pasfinal
else
  rel_pos_vel(:,$) = [];
  t(:,$) = [];
end

endfunction