//  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 [deccdt,dpomdt] = CL_op_frozenOrbitDer(sma,ecc,inc,pom, er,mu,j1jn)
// Derivatives of eccentricity and argument of periapsis with respect to time
//
// Calling Sequence
// [deccdt,dpomdt] = CL_op_frozenOrbitDer(sma,ecc,inc,pom [,er,mu,j1jn])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the time derivatives of eccentricity and argument of periapsis 
// resulting from the gravitational effects due to J2 and J3.</p>
// <p>The time derivative of the argument of periapsis is undefined (%nan) if the orbit is circular (ecc = 0) or equatorial
// (inc = 0 or pi).</p>
// <p>(See the formulas, where: R = equatorial radius, n = keplerian mean motion)</p>
// <p><inlinemediaobject><imageobject><imagedata fileref="driftJ2J3.gif"/></imageobject></inlinemediaobject></p>
// <p></p></listitem>
// <listitem>
// <p>Warning :</p>
// <p> - The input argument "zonals" is deprecated as of CelestLab v3.0.0. It has been replaced by "j1jn".</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// sma: Semi-major axis [m] (1x1 or 1xN)
// ecc: Eccentricity (1x1 or 1xN)
// inc: Inclination [rad] (1x1 or 1xN)
// pom: Argument of pariapsis [rad] (1x1 or 1xN)
// er: (optional) Equatorial radius [m] (default is %CL_eqRad)
// mu: (optional) Gravitational constant [m^3/s^2] (default value is %CL_mu)
// j1jn: (optional) Vector of zonal coefficients J1 to Jn, troncated to J3. Default is %CL_j1jn(1:3)). (1xNz)
// deccdt: Time derivative of eccentricity [s^-1] (1xN)
// dpomdt: Time derivative of argument of periapsis [rad/s] (1xN)
//
// Authors
// CNES - DCT/SB
// 
// Bibliography
// 1) "Frozen orbits in the J2+J3 problem", Krystyna Kiedron and Richard Cook, AAS 91-426. 
//
// Examples
// sma = [7000.e3, 73000.e3];
// inc = CL_deg2rad([51.6, 91.6]);
// [ecc,pom] = CL_op_frozenOrbit(sma,inc);
// ecc = 0.999*ecc;
// pom = 0.999*pom;
// [deccdt,dpomdt] = CL_op_frozenOrbitDer(sma,ecc,inc,pom)

// Declarations:

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

[sma, ecc, inc, pom] = CL__checkInputs(sma, 1, ecc, 1, inc, 1, pom, 1); 

if (find(ecc < 0 | ecc >= 1 | sma <= 0) <> []); 
  CL__error("Invalid arguments"); 
end

if (length(j1jn) < 3)
  CL__error('j1jn must be a vector of size 3 (or more)');
end

j2 = j1jn(2);
j3 = j1jn(3);

sini = sin(inc);
// n = mean motion
n = CL_kp_params('mm',sma,mu); 

deccdt = -0.375 * (er./sma).^3 .* n*j3 ./ (1-ecc.^2).^2 .* sini.*cos(pom) .* (4-5*sini.^2) ;

// check possible division by 0 (result is then %nan)
I = find(abs(sini) < %eps | abs(ecc) < %eps);
ecc(I) = %nan; // makes the calculation possible
sini(I) = %nan;  

dpomdt = 0.75 * (er./sma).^2 .* n*j2 ./ (1-ecc.^2).^2 .*(4-5*sini.^2) + ...
        -0.375 * (er./sma).^3 .* n*j3 ./ (ecc.*(1-ecc.^2).^3) .* sin(pom)./sini .* ...
        ((5*sini.^2-4).*sini.^2 + ecc.^2 .* (4-35*sini.^2 .* (1-sini.^2)));
// NB : the only modification of this formula compared to the one in the documentation is
//      the last cosi^2 replaced by (1-sini^2) (in order to improve performance)


endfunction