//  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 kepculated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [oe_t2] = CL__ex_propagJ2sec(type_oe, t1, oe_t1, t2, er, mu, j2)
// Orbit propagation using "secular J2" model ( no conversion to kep ty
//
// Calling Sequence
// [oe_t2] = CL__ex_propagJ2sec(type_oe, t1, oe_t1, t2, er, mu, j2)
//
// Description
// <itemizedlist><listitem>
// <p>Propagates orbital elements considering the mean secular effects due to J2 only. </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> - This function works for elliptical orbits only. </p>
// <p> - Values of equatorial radius, gravitational constant and second zonal harmonic are Celestlab default ones (%CL_eqRad, %CL_mu and %CL_j1jn(2), repectively) </p>
// </listitem>
// <listitem>
// <p>The type of orbital elements can be : "kep", "cir", "cireq" or "equin" 
// (see <link linkend="Orbital elements">Orbital elements</link>) </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// type_oe: (string) Type of orbital elements: "kep", "cir", "cireq" or "equin".
// t1: Initial time [days] (1x1 or 1xN).
// kep_t1: Orbital elements at time t1 (6x1 or 6xN).  
// t2: Final time [days] (1xN or 1x1).  
// er: Equatorial radius [m] 
// mu: Gravitational constant [m^3/s^2] 
// j2: (Unnormalized) zonal coefficient (second zonal harmonic) 
// kep_t2: Orbital elements propagated to time t2 (6xN).

// Inputs checking
// adjust sizes : same for t1 and kep_t1
[t1, oe_t1, N1] = CL__checkInputs(t1, 1, oe_t1, 6); 
[t2, N2] = CL__checkInputs(t2, 1); // check nb of rows == 1

if ~(N1 == 1 | N2 == 1 | N1 == N2)
   CL__error('Wrong size of input arguments'); 
end
N = max(N1, N2);
    
// Computation of the drifts (mean time derivatives) of orbital elements due to J2
// Conversion to kep is used to compute sma, ecc and inc only. 
if (type_oe == "kep")
  kep_t1 = oe_t1; 
else
  kep_t1 = CL_oe_convert(type_oe, "kep", oe_t1);  
end  

[dpomdt, dgomdt, danmdt] = CL_op_driftJ2(kep_t1(1,:), kep_t1(2,:), kep_t1(3,:), er, mu, j2);


// propagation time (seconds)
delta_t = (t2-t1) * 86400; 

// Orbital elements at final times
if (type_oe == 'kep')
  oe_t2 = [ oe_t1(1,:) .* ones(1,N); ..
            oe_t1(2,:) .* ones(1,N); ..
            oe_t1(3,:) .* ones(1,N); ..
            oe_t1(4,:) + dpomdt .* delta_t; ..
            oe_t1(5,:) + dgomdt .* delta_t; ..
            oe_t1(6,:) + danmdt .* delta_t ];
            
elseif (type_oe == 'cir')
  oe_t2 = [ oe_t1(1,:) .* ones(1,N); ..
            oe_t1(2,:) .* cos(dpomdt .* delta_t) - oe_t1(3,:) .* sin(dpomdt .* delta_t); ..
            oe_t1(3,:) .* cos(dpomdt .* delta_t) + oe_t1(2,:) .* sin(dpomdt .* delta_t); ..
            oe_t1(4,:) .* ones(1,N); ..
            oe_t1(5,:) + dgomdt .* delta_t; ..
            oe_t1(6,:) + (dpomdt + danmdt) .* delta_t ];

elseif (type_oe == 'cireq' | type_oe == 'equin')
  oe_t2 = [ oe_t1(1,:) .* ones(1,N);
            oe_t1(2,:) .* cos((dpomdt + dgomdt) .* delta_t) - oe_t1(3,:) .* sin((dpomdt + dgomdt) .* delta_t); ..
            oe_t1(3,:) .* cos((dpomdt + dgomdt) .* delta_t) + oe_t1(2,:) .* sin((dpomdt + dgomdt) .* delta_t); ..
            oe_t1(4,:) .* cos(dgomdt .* delta_t) - oe_t1(5,:) .* sin(dgomdt .* delta_t); ..
            oe_t1(5,:) .* cos(dgomdt .* delta_t) + oe_t1(4,:) .* sin(dgomdt .* delta_t); ..
            oe_t1(6,:) + (dpomdt + dgomdt + danmdt) .* delta_t ];

else
  CL__error("Invalid type of orbital elements"); 
end

endfunction