//  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_kepler(t1,kep_t1,t2, mu)
// Keplerian orbit propagation
//
// Calling Sequence
// [kep_t2, jac] = CL_ex_kepler(t1,kep_t1,t2 [,mu])
//
// Description
// <itemizedlist><listitem>
// <p>Propagates orbital elements considering the central force 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>
// <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 any type of orbit and any type of orbital elements
// ("kep", "cir", "cireq" or "equin"). </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).  
// mu : (optional) Gravitational constant. [m^3/s^2]. Default value is %CL_mu
// kep_t2: Orbital elements propagated to time t2. (6xN)
// jac: Transition matrix = d(kep_t2)/d(kep_t1) (6x6xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Mecanique spatiale, CNES - Cepadues 1995, Tome I
//
// See also
// CL_ex_secularJ2
//
// Examples
// // Example 1 : 
// t1 = 20000 // initial time (days)
// kep_t1 = [ 10000.e3; 0.7; 0; 0; 0; 0]
// T = CL_kp_params('per', kep_t1(1,:)) // orbital period (s)
// t2 = t1 + (1:5) * (T/2) / 86400 // final times (days)
// kep_t2 = CL_ex_kepler(t1, kep_t1, t2)
//
// // Example 2 : 
// t1 = 0  // initial time (arbitrary origin)
// kep_t1 = [ [7.e6; 0.;  %pi/2; 0; 0; 0 ], ..
//            [8.e6; 0.1; %pi/2; 0; 0; 0 ] ] 
// t2 = [1,2] // final times (1 or 2 days later)
// kep_t2 = CL_ex_kepler(t1, kep_t1, t2)
//

// Declarations:

// Code:
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

// 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); // size of result

sma = kep_t1(1,:);  // semi major axis

if (find(sma <= 0) <> []) 
  CL__error("Invalid orbital elements"); 
end

// mean motion
mm = CL_kp_params('mm', sma, mu=mu);  

// 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,:) .* ones(1,N); 
           kep_t1(5,:) .* ones(1,N); 
           kep_t1(6,:) + mm .* delta_t ];  
           
// transition matrix if there are at least 2 output arguments
if (argn(1) > 1)
  jac = repmat(eye(6,6),[1,1,N]); // identity hypermatrix
  jac(6,1,1:N) = -(3/2) * (mm./sma) .* delta_t;
end
           
endfunction