//  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 [kep_t2, jac] = CL_ex_secularJ2(t1,kep_t1,t2, er,mu,j2)
// Orbit propagation using "secular J2" model
//
// Calling Sequence
// [kep_t2, jac] = CL_ex_secularJ2(t1,kep_t1,t2 [,er,mu,j2])
//
// Description
// <itemizedlist><listitem>
// <p>Propagates orbital elements considering the mean secular effects due to J2 only. </p>
// <p>The jacobian of the transformation (i.e. the transition matrix) is also 
// optionally computed. </p>
// <p></p></listitem>
// <listitem>
// <p>The type of orbital elements is the following:</p>  
// <p><inlinemediaobject><imageobject><imagedata fileref="kep_par.gif"/></imageobject></inlinemediaobject></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>
// </listitem>
// </itemizedlist>
//
// Parameters
// 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: (optional) Equatorial radius [m] (default is %CL_eqRad)
// mu: (optional) Gravitational constant [m^3/s^2] (default value is %CL_mu)
// j2: (optional) (Unnormalized) zonal coefficient (second zonal harmonic) (default is %CL_j1jn(2))
// kep_t2: Orbital elements propagated to time t2. (6xN)
// jac: Transition matrix = d(kep_t2)/d(kep_t1) (6x6xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_ex_propagate
// CL_op_driftJ2
//
// Examples
// // propagation of one satellite to several times:
// sma = 7.e6;
// ecc = 0.001;
// inc = CL_deg2rad(98);
// pom = CL_deg2rad(90);
// gom = 0;
// anm = 0;
// t1 = 21915;
// t2 = t1:0.1:t1+1;
// kep_t1 = [sma; ecc; inc; pom; gom; anm];
// [kep_t2] = CL_ex_secularJ2(t1, kep_t1, t2);
//
// // propagation of 2 element sets to one final time:
// t2 = t1+1;
// t1 = [t1, t1+0.5];
// kep_t1 = [kep_t1, [sma+1; ecc+0.1; inc; pom; gom; anm]];
// [kep_t2] = CL_ex_secularJ2(t1, kep_t1, t2);

// Gestion arguments optionnels

// Declarations:

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

// Inputs checking
// adjust sizes : same for t1 and kep_t1
[t1, kep_t1, N1] = CL__checkInputs(t1, 1, kep_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);

// compute_jac: true if transition matrix is computed
compute_jac = %t; 
if (argn(1) == 1); compute_jac= %f; end

// drifts due to J2
if (compute_jac)
  [pomdot, gomdot, anmdot, dpomdotdaei, dgomdotdaei, danmdotdaei] = ..
         CL_op_driftJ2(kep_t1(1,:), kep_t1(2,:), kep_t1(3,:), er=er, mu=mu, j2=j2);
else
  [pomdot, gomdot, anmdot] = CL_op_driftJ2(kep_t1(1,:), kep_t1(2,:), kep_t1(3,:), er=er, mu=mu, j2=j2);
end

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

// orbital elements at final times
kep_t2 = [ kep_t1(1,:) .* ones(1,N);
           kep_t1(2,:) .* ones(1,N);
           kep_t1(3,:) .* ones(1,N);
           kep_t1(4,:) + pomdot .* delta_t;
           kep_t1(5,:) + gomdot .* delta_t;
           kep_t1(6,:) + anmdot .* delta_t ];


// transition matrix if there are at least 2 output arguments
// jac(i,j,:) = d(kep_t2(i,:)) / d(kep_t1(j,:))
if (compute_jac)
  jac = repmat(eye(6,6),[1,1,N]); // identity hypermatrix
  for (k = 1:3)
    jac(4,k,:) = dpomdotdaei(k,:) .* delta_t; 
    jac(5,k,:) = dgomdotdaei(k,:) .* delta_t; 
    jac(6,k,:) = danmdotdaei(k,:) .* delta_t; 
  end
end

           
endfunction