//  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 [incdot] = CL_op_incdotSsoCirc(mltan, sma, inc, mu_b, r_b, inc_b, er, mu, j2, rotr_pla_sun)
// Long-term time derivative of inclination (thirdbody perturbation) 
//
// Calling Sequence
// [incdot] = CL_op_incdotSsoCirc(mltan, sma [, inc, mu_b, r_b, inc_b, er, mu, j2, rotr_pla_sun])
//
// Description
// <itemizedlist><listitem>
// <p>Long-term (that is doubly averaged) time derivative of inclination due to thirdbody perturbation 
// for a circular and Sun-synchronous orbit. </p>
// <p>The formula is the result of an analytical calculation assuming a simplified model for the orbit of the 
// perturbing body (circular orbit, constant inclination). </p>
// <inlinemediaobject><imageobject><imagedata fileref="incdotSsoCirc.png"/></imageobject></inlinemediaobject>
// <p></p>
// <p><b>Notes</b>:</p>
// <p>- Inclination can be [] or not specified in which case it is recomputed. </p>
// <p>- The function could theoretically be used for another body than the Sun, 
// provided the orbit plane (RAAN) rotates as the same angular rate as
// the perturbing body. </p>
// <p>- The default value used for the apparent inclination of the orbit of the 
// perturbing body is the value at J2000 in ECI (23.44 deg). In reality the value is not
// constant: there is a long-term drift of about -0.014 deg/century and
// variations of maximum amplitude about 0.005 deg. </p> 
// </listitem>
// </itemizedlist>
//
// Parameters
// mltan: Mean local time of the ascending node [h] (1xN or 1x1)
// sma: Mean semi-major axis [m] (1xN or 1x1)
// inc: (optional) Mean inclination [rad]. Default is []. (1xN or 1x1)
// mu_b: (optional) Gravitational constant of perturbing body [m^3/s^2]. Default is %CL_body.Sun.mu (1x1)
// r_b: (optional) Distance from central body to perturbing body [m]. Default is %CL_au. (1x1 or 1xN)
// inc_b: (optional) Inclination of orbit of perturbing body [rad]. Default value is 23.44 deg. (1x1 or 1xN)
// er: (optional) Equatorial radius [m]. Default is %CL_eqRad. 
// mu: (optional) Gravitational constant of central body [m^3/s^2]. Default is %CL_mu 
// j2: (optional) Second zonal harmonic term. Default is %CL_j1jn(2)
// rotr_pla_sun : (optional) Mean apparent rotation rate of the Sun around the planet. Default is %CL_rotrBodySun
// incdot: Time derivative of inclination [rad/s] (1xN)
//
// Examples
// mltan = 0 : 0.1 : 24; // hours
// sma = 6378.e3 + 700.e3; // m
//
// // Inclination average time derivative (rad/s) 
// incdot = CL_op_incdotSsoCirc(mltan, sma); 
// 
// // Plot result (deg/year)
// scf(); 
// plot(mltan, incdot * (180/%pi) * 86400 * 365.25); 

// Initialize optional inputs
if (~exists("inc","local")); inc = []; end;
if (~exists("mu_b","local")); mu_b = CL__dataGetEnv(["body", "Sun", "mu"]); end
if (~exists("r_b","local")); r_b = CL__dataGetEnv("au"); end
if (~exists("inc_b","local")); inc_b = 23.44 * %pi / 180; end
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
if (~exists("rotr_pla_sun", "local")); rotr_pla_sun = CL__dataGetEnv("rotrBodySun"); end

// Check number of input arguments
if (argn(2) < 2)
  CL__error("Too few arguments: at least 2 expected"); 
end

// Eccentricity assumed to be 0
ecc = 0; 

// Compute inclinaison (SSO) if not defined 
if (inc == [])
  inc = CL_op_ssoJ2("i", sma, ecc, er, mu, j2, rotr_pla_sun); 
end

// Check/resize input arguments
[sma, inc, mltan, inc_b, r_b] = CL__checkInputs(sma, 1, inc, 1, mltan, 1, inc_b, 1, r_b, 1);

// Check inputs validity
I = find (sma <= 0); 
if (I <> []); CL__error("Semi major axis out of range"); end

I = find (inc < 0 | inc > %pi + %eps); // margin in case !!!
if (I <> []); CL__error("Inclination out of range"); end

I = find (inc_b < 0 | inc_b > %pi + %eps); // margin in case !!!
if (I <> []); CL__error("Inclination of perturbing body out of range"); end

I = find (r_b <= 0); 
if (I <> []); CL__error("Distance to perturbing body out of range"); end

// Mean motion (Keplerian value)
// as not checked that the "J2" value gives a more accurate result   
// mm = CL_op_paramsJ2("mm", sma, ecc, inc, er, mu, j2, rotr_pla_sun);
mm = CL_kp_params("mm", sma, mu);

// Time derivative of inclination
h = (mltan - 12) * %pi/12; // local time angle
incdot = 0.75 * mu_b * (cos(inc_b/2)).^4 .* sin(inc) .* sin(2*h) ./ (mm .* r_b.^3);

endfunction