//  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 [pv] = CL_cw_propagate(t0,pv0,t,alt,acc,er,mu)
// Propagation using Clohessy-Whiltshire model
//
// Calling Sequence
// pv = CL_cw_propagate(t0,pv0,t,alt [,acc,er,mu])
//
// Description
// <itemizedlist><listitem>
// <p>Propagates using clohessy-Whiltshire equations (including a constant relative acceleration).</p> 
// <p></p>
// <p>The reference frame used is the target's LVLH local frame (origin = target). </p>
// <p>The reference altitude <b>alt</b> is the altitude of the target = semi major axis of target's (circular) orbit
// minus equatorial radius. </p>
// <p></p>
// <p><b>acc</b> is a (constant) differential inertial acceleration between t0 and t with components 
// in the reference frame. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// t0: Initial time [days] (1x1 or 1xN)
// pv0: Initial (relative) position and velocity vectors in reference frame [m,m/s] (6x1 or 6xN)
// t: Final time [days] (1xN or 1x1)
// alt: Reference altitude (= altitude of target) [m] (1x1 or 1xN)
// acc: (optional) Differential (constant) acceleration in reference frame (m/s^2). Default is []. (3x1 or 3xN)
// er: (optional) Equatorial radius [m]. Default is %CL_eqRad
// mu: (optional) Gravitational constant [m^3/s^2]. Default value is %CL_mu
// pv: Relative position and velocity vectors at time t in reference frame [m,m/s] (6xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_cw_Mmatrix
// CL_cw_Nmatrix
//
// Bibliography
// 1) Mecanique spatiale, CNES - Cepadues 1995, Tome II
//
// Examples
// alt = 450.e3;
// t0 = 0;
// pv0 = [1;0;0;1;0;0];  
// t = (100:100:500)/86400; 
// pv = CL_cw_propagate(t0, pv0, t, alt)


// Declarations:

// Code:
if (~exists('acc','local')); acc=[]; end
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

// check / resize
// Note: it should not be necessary to resize everything, but doing so leads to 
// simpler code. 
[t0, pv0, t, alt, acc, N] = CL__checkInputs(t0,1,pv0,6,t,1,alt,1,acc,3);

// no input dates
if (t == [])
  pv = []; 
  return; 
end

delta_t = (t - t0) * 86400; 

// ballistic propagation
mat_M = CL_cw_Mmatrix(alt, delta_t, er, mu);
pv = mat_M * pv0; 

// add effect of acceleration if not empty
if (acc <> [])
  mat_N = CL_cw_Nmatrix(alt, delta_t, er, mu);
  pv = pv + mat_N * acc; 
end


endfunction