//  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 [sma,inc] = CL_op_repeat2smaInc(N,P,Q,ecc,sso, incInput,er,mu,j2,rotr_pla,rotr_pla_sun)
// Semi major axis and inclination corresponding to (N,P,Q)  
//
// Calling Sequence
// [sma,inc] = CL_op_repeat2smaInc(N,P,Q,ecc,sso [,incInput,er,mu,j2,rotr_pla,rotr_pla_sun])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the semi major axis (and inclination for sunsynchronous orbits) of repeat (also called phased) orbits 
// with given repeat parameters N,P,Q.</p>
// <p>The 3 integers N,P,Q are such that: </p>
// <p>Tc (duration of repeat cycle) = (N*Q+P)*T = Q * Tr, where T is the orbit period (more exactly the mean nodal period), 
// and Tr is the revolution period of the planet with respect to the orbit's (ascending or descending) node. </p>
// <p>Then: </p>
// <p>Q (=Tc/Tr) is the number of planet revolutions per repeat cycle.</p>
// <p>N+P/Q (=Tr/T) is the number of orbits per planet revolutions with respect to the orbit's (ascending or descending) node.</p>
// <p>N*Q+P (= Tc/T) is the (whole) number of orbits per repeat cycle.</p>
// <p>Note that Tr depends on the orbit parameters because the node drifts under the effect of J2.</p>
// <p> If sso = 1, the orbits are sun-synchronous (the inclination is then computed)</p>
// <p> If sso = 0, the function uses incInput as the inclination to be used. </p>
// </listitem></itemizedlist>
//
// Parameters
// N: Values of N (1x1 or 1xNt)
// P: Values of P (1x1 or 1xNt)
// Q: Values of Q (1x1 or 1xNt)
// ecc: Eccentricity (1x1 or 1xNt)
// sso: 1 = sun synchronous orbits, 0 = fixed inclination (1x1)
// incInput : (mandatory in case sso=0, unused otherwise) inclinations (rad) (1x1 or 1xNt)
// 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 coefficient (default is %CL_j1jn(2))
// rotr_pla : (optional) rotation rate of the planet (default is %CL_rotrBody) (1x1)
// rotr_pla_sun : (optional) mean apparent rotation rate of the Sun around the planet (default is %CL_rotrBodySun) (1x1)
// sma: Semi major axis (m). sma is %nan if value cannot be computed.  (1xNt)
// inc: Inclination (rad). inc is %nan if value cannot be computed. (1xNt)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_op_searchRepeatOrbits
// CL_op_repeatGroundTracks
// CL_op_ssoJ2
//
// Examples
// N = 15;
// P = 3;
// Q = 16;
// ecc = 0.01;
// sso = 0;
// incInput = CL_deg2rad(98);
// [sma,inc] = CL_op_repeat2smaInc(N,P,Q,ecc,sso,incInput)
//
// N = [15,15];
// P = [3,5];
// Q = [16,6];
// ecc = [0.01, 0.001];
// sso = 1 ;
// [sma,inc] = CL_op_repeat2smaInc(N,P,Q,ecc,sso)


  function [f,dfda] = fct_xnorb(sma,ecc,sso,inc,er,mu,j2,rotr_pla,rotr_pla_sun)
  // fct_xnorb: 
  // f = number of orbits / planet revolution 
  // dfda = derivative of f / semi major axis
  // NOTE : inc must be initialized in all cases
  
    [pomp,gomp,anmp,dpomp,dgomp,danmp] = CL_op_driftJ2(sma,ecc,inc,er,mu,j2);
    alphap = pomp+anmp; 
    dalphap = dpomp+danmp; 

    if (sso == 1)
      // sun-synchronous => gomp == rotr_pla_sun (constant)
      // derivative : df/da + df/di * di/da
      // with : di/da computed assuming Sun-synchronicity (cosi = K(e) * a^7/2 )
      // --> -sini*di = K(e) * 7/2 * a^5/2 * da and replace replace K(e) by cosi / a^7/2
      // --> di/da = -7/2 / (tani * a)
      f = alphap / (rotr_pla - rotr_pla_sun); 
      dfda = ( dalphap(1,:) + dalphap(3,:) .* ((-7/2) ./ (tan(inc) .* sma)) ) ..
                                                         / (rotr_pla-rotr_pla_sun);

    else
      f = alphap ./ (rotr_pla - gomp);
      dfda = ( dalphap(1,:) .* (rotr_pla-gomp) + alphap.*dgomp(1,:) ) ..
                                                        ./ (rotr_pla-gomp).^2;
    end
  
  endfunction


// 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
if (~exists("rotr_pla", "local")); rotr_pla = CL__dataGetEnv("rotrBody"); end
if (~exists("rotr_pla_sun", "local")); rotr_pla_sun = CL__dataGetEnv("rotrBodySun"); end

// Verification des entrees
I = find(ecc < 0 | ecc >= 1);
if (~isempty(I)); CL__error("Eccentricity out of range"); end
if (sso <> 0 & sso <> 1); CL__error("sso: wrong value"); end

if (sso == 0)
  if ~exists('incInput','local')
    CL__error('Inclination missing (sso == 0)');
  end
  I = find(incInput < 0 | incInput > %pi+%eps);
  if ~isempty(I); CL__error("Inclination out of range"); end
else
  incInput = [];
end


dimN = length(N);
dimP = length(P);
dimQ = length(Q);
dimEcc = length(ecc);
dimInc = length(incInput); // == 0 if sso == 1
Nmax = max(dimN,dimEcc,dimInc);

if ~(dimN == dimP & dimN == dimQ); CL__error("Wrong lengths for N,P or Q"); end

// empty input interval => returns empty
if (Nmax == 0) 
  sma = []; 
  inc = []; 
  return; 
end

if ~((dimEcc == 1 | dimEcc == Nmax) & (dimInc == 1 | dimInc == Nmax | sso == 1) & (dimN == 1 | dimN == Nmax))
  CL__error("Invalid input arguments sizes")
end

I = find (Q <= 0 | (P < 0 | P >= Q) | N < 0);
if (~isempty(I)); CL__error("N,P,Q out of range"); end


// adjust sizes of N, P, Q, ecc and inc
ecc = ecc .* ones(1,Nmax);
inc = incInput .* ones(1,Nmax);
N = N .* ones(1,Nmax);
P = P .* ones(1,Nmax);
Q = Q .* ones(1,Nmax);


// more initializations
fref = (N + P./Q);  // target number of orbits / planet revolution
f = zeros(N);   
fprime = zeros(N);

// initializes sma (aproximate value)
nmoy = fref .* rotr_pla;   // approximate orbit's mean motion
sma = (mu./nmoy.^2).^(1/3);

iterMax = 30;
iter = 1;
prec = 1.e-13;  // accuracy on N+P/Q

// convergence process
I = 1 : Nmax;   // indices to be computed 

while( ~isempty(I) & iter <= iterMax)

  if (sso == 1)
    inc = CL_op_ssoJ2('i',sma,ecc,er,mu,j2,rotr_pla_sun);
    ind_nan = find(isnan(inc));  // if nan present : remove index from I 
    sma(I(ind_nan)) = %nan;
    I(ind_nan) = [];
  end

  [f(I),fprime(I)] = fct_xnorb(sma(I),ecc(I),sso,inc(I),er,mu,j2,rotr_pla,rotr_pla_sun);

  // better sma value
  sma(I) = sma(I) + (fref(I) - f(I)) ./ fprime(I);

  I = find( abs(f - fref) > prec & sma > 0); // discard negative values of sma
  iter = iter+1;

end

if (iter > iterMax)
  CL__error("No convergence (max number of iterations reached)")
end

I = find( sma <= 0);
inc(I) = %nan; 
sma(I) = %nan;

endfunction