//  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 [pomdot,gomdot,anmdot,dpomdotdaei,dgomdotdaei,danmdotdaei] = CL_op_driftJ2(sma,ecc,inc, er,mu,j2)
// Drifts (mean time derivatives) of orbital elements due to J2 + central force
//
// Calling Sequence
// [pomdot,gomdot,anmdot ,dpomdotdaei,dgomdotdaei,danmdotdaei] = CL_op_driftJ2(sma,ecc,inc [,er,mu,j2])
//
// Description
// <itemizedlist><listitem>
// <p>Computes: </p>
// <p>- the secular drifts (pomdot=d(pom)/dt, gomdot=d(gom)/dt, anmdot=d(anm)/dt) on the keplerian elements, 
// where pom=argument of periapsis, gom=right ascension of ascending node,anm=mean anomaly, 
// considering the effects of the central (Keplerian) force and the first zonal harmonic J2, </p>
// <p>- the partial derivatives of these drifts with respect to semi major axis, eccentricity and inclination.</p>
// <p>(See the formulas, where: R = equatorial radius, n = keplerian mean motion</p>
// <inlinemediaobject><imageobject><imagedata fileref="driftJ2.gif"/></imageobject></inlinemediaobject>
// </listitem>
// </itemizedlist>
//
// Parameters
// sma: Semi-major axis [m] (1x1 or 1xN)
// ecc: Eccentricity  (1x1 or 1xN)
// inc: Inclination [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)
// j2: (optional) Second zonal harmonic) (default is %CL_j1jn(2))
// pomdot: Mean time derivative of argument of periapsis [rad/s] (1xN)
// gomdot: Mean time derivative of right ascension of ascending node [rad/s] (1xN)
// anmdot: Mean time derivative of mean anomaly [rad/s] (1xN)
// dpomdotdaei: Partial derivative of pomdot with respect to sma, ecc and inc (3xN)
// dgomdotdaei: Partial derivative of gomdot with respect to sma, ecc and inc (3xN)
// danmdotdaei: Partial derivative of anmdot with respect to sma, ecc and inc (3xN)
//
// Bibliography
// 1) CNES - MSLIB FORTRAN 90, Volume E (me_deriv_secul_J2)
//
// Authors
// CNES - DCT/SB
//
// Examples
// sma = [42000.e3, 7000.e3];
// ecc = [0.001, 0.001];
// inc = CL_deg2rad([51.6, 51.6]);
// [pomdot,gomdot,anmdot] = CL_op_driftJ2(sma,ecc,inc)
//

// 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

lhs = argn(1) // number of output arguments
if (lhs > 6); CL__error('bad number of output arguments'); end 

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

I = find (sma <= 0); 
if (~isempty(I)); CL__error('semi major axis out of range'); end

I = find (ecc < 0 | ecc >= 1); 
if (~isempty(I)); CL__error('eccentricity out of range'); end

I = find (inc < 0 | inc > %pi+%eps); // margin in case !!!
if (~isempty(I)); CL__error('inclination out of range'); end


// -------- 
// formulas
// --------
K = 0.75 * j2 * er^2 * sqrt(mu);  // common coefficient

n = sqrt(mu ./ sma.^3); // mean motion 
cosi = cos(inc); 
f2 = 1 - ecc.^2;
f = sqrt(f2); 
D = sma.^(7.0 / 2) .* f2.^2; // denominator in following expressions

// 1) secular drifts of parameters
pomdot = K * (5 * cosi.^2 - 1) ./ D;        // arg of perigee
gomdot = -2 * K * cosi ./ D;                // RAAN
anmdotJ2 = K * (3 * cosi.^2 - 1) .* f ./ D;   // mean anomaly: effect of J2 only
anmdot = anmdotJ2 + n; // mean anomaly: effect of J2 + central body

// 2) derivatives
if lhs > 3

  dpomdotdaei = zeros(3,N);
  dgomdotdaei = zeros(3,N);
  danmdotdaei   = zeros(3,N);
  sini     = sin(inc); 

  // derivatives with respect to semi major axis
  dpomdotdaei(1,:) = -3.5 * pomdot ./ sma;                 // d(pomdot)/d(sma)
  dgomdotdaei(1,:) = -3.5 * gomdot ./ sma;                 // d(gomdot)/d(sma)
  danmdotdaei(1,:)   = (-3.5 * anmdotJ2 - 1.5 * n) ./ sma;  // d(anmdot)  /d(sma)

  // derivatives with respect to eccentricity
  dpomdotdaei(2,:) = 4 * pomdot .* ecc ./ f2;              // d(pomdot)/d(exc)
  dgomdotdaei(2,:) = 4 * gomdot .* ecc ./ f2;              // d(gomdot)/d(exc)
  danmdotdaei(2,:)   = 3 * anmdotJ2 .* ecc ./ f2;           // d(anmdot)  /d(exc)

  // derivatives with respect to inclination
  dpomdotdaei(3,:) = -10 * K * cosi .* sini ./ D;          // d(pomdot)/d(inc)
  dgomdotdaei(3,:) = 2 * K * sini ./ D;                    // d(gomdot)/d(inc)
  danmdotdaei(3,:)   = -6 * K * cosi .* sini .* f ./ D;    // d(anmdot)  /d(inc)

end

endfunction