//  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 [ecc,pom] = CL_op_frozenOrbit(sma,inc, er,j1jn)
// Eccentricity and argument of periapsis of a frozen orbit
//
// Calling Sequence
// [ecc,pom] = CL_op_frozenOrbit(sma,inc [,er,j1jn])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the eccentricity and argument of periapsis of a frozen
// orbit, so that the mean value of eccentricity and argument of periapsis remain 
// constant over time. </p>
// <p></p>
// <p> Freezing the orbit is possible by balancing the effects of harmonic even 
// and odd terms of the potential (J1=0, J2, J3, ...). </p>  
// <p> The argument of periapis returned is 90 degrees or -90 degrees. The eccentricity 
// is usually small: of the order of 1.e-3 for the Earth. </p>
// <p> If j1jn are limited to J2 and J3, the frozen orbit eccentricity is given by: </p>
// <p> eccg = -0.5 * sin(inc) * (Req/sma) * (J3/J2)  </p>
// <p></p></listitem>
// <listitem>
// <p><b>Note</b>: </p>
// <p>The relative accuracy decreases as the number of zonal terms increases because of 
// numerical errors. The result can be very inaccurate for 40 terms or more. </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] (1xN or 1x1) 
// inc: Inclination [rad] (1xN or 1x1)
// er: (optional) Equatorial radius [m]. Default is %CL_eqRad.
// j1jn: (optional) Vector of zonal coefficients J1 to Jn (1xNz or Nzx1). Default is %CL_j1jn(1:3).  
// ecc: Mean eccentricity (1xN)
// pom: Mean argument of periapsis [rad] (1xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// sma = [7000.e3, 7300.e3];
// inc = CL_deg2rad(98);
// [ecc,pom] = CL_op_frozenOrbit(sma,inc)
// j1jn = CL_dataGet("j1jn");
// [ecc,pom] = CL_op_frozenOrbit(sma,inc,j1jn=j1jn(1:20))


// Declarations:

// ----------------------------------------------------------
// internal function: eccentricity/argument of perigee
// of frozen orbit
// sma, inc: semi major axis and inclination (same size = 1xN)
// er: equatorial radius
// j1jn: zonal terms to be used (at least: J1, J2, J3 required)
// ----------------------------------------------------------

function [eccg,pomg] = frozenOrbit(sma,inc,er,j1jn)  

nzmax = length(j1jn);

// ancillary variables  
sini = sin(inc);
cosi2 = 1 - sini.^2; // cos(inc).^2
ra = er ./ sma;

// sin(i) .^ k , k=0..nzmax (nzmax+1 rows) 
MSIN = (ones(nzmax+1,1) * sini) .^ ((0:nzmax)' * ones(sini));  

// precomputation of factorial terms
MFACT = factorial(0:2*nzmax);

// computation of: sum( c(k) * sin(i) .^ L(k) ), k=1,2...
// c and L: same size (not empty)
function [s] = sum_sini(c,L)
  s = sum((c' * ones(sini)) .* MSIN(L+1,:),'r');
endfunction

// computation of: factorial(L), L not empty
function [fact] = f(L)
  fact = MFACT(L+1);
endfunction


den = zeros(sini); 
num = zeros(sini); 

for i = 2 : nzmax
  if (modulo(i,2) == 0)

    k = i/2; 
    L = 0:k;

    X = ((-1).^L).* f(2*(i-L)) ./ (f(L).*f(i-L).*(4.^(i-L)));

    ta = X ./ f(k - L) .^ 2;
    sy = -k * (2*k+1) * sum_sini(ta, i-2*L);

    L1 = L(1:$-1); // not empty
    sx = 2 * cosi2 .* sum_sini((k-L1) .* ta(1:$-1), i-2*L1-2);

    sz = zeros(sini); 
    if (k > 1)
      te = X(1:$-1) ./ (f(k-1-L1) .* f(k+1-L1));
      sz = -(k-1) * (2*k-1) * sum_sini(te, i-2*L1);
    end

    den = den + (sx+sy+sz) .* (ra.^i) .* j1jn(i);

  else
 
    k = (i-1)/2; 
    L = 0:k;

    td = (((-1).^L).*f(2*(i-L))) ./ (f(L).*f(i-L).*f(k-L).*f(k+1-L).*(4.^(i-L)));

    s = sum_sini(td, i-2*L);
    num = num + 2 * k * s .* (ra.^i) .* j1jn(i);

  end
end

eccg = num ./ den;  
pomg = (%pi/2) * ones(eccg); 

// if eccg < 0, the argument of perigee is in fact -90 deg
// (the value of : e * sin(pom) remains the same) 

I = find(eccg < 0); 
eccg(I) = -eccg(I); // becomes positive
pomg(I) = -pomg(I); // becomes -90 deg

endfunction


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

Nsma = size(sma,2); 
Ninc = size(inc,2); 

N = max(Nsma, Ninc); // max number of columns

if ((Nsma <> 1 & Nsma <> N) | (Ninc <> 1 & Ninc <> N))
   CL__error('Invalid arguments sizes'); 
end

if (Nsma < N); sma = sma * ones(1,N); end
if (Ninc < N); inc = inc * ones(1,N); end


if (length(j1jn) < 3)
  CL__error('Invalid number of zonal terms (3 or more required)');
end

// calculation (internal function)
[ecc,pom] = frozenOrbit(sma,inc,er,j1jn); 

endfunction