//  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 [params] = CL_op_searchRepeatOrbits(smaMin,smaMax,QMin,QMax,ecc,sso, incInput,er,mu,j2,rotr_pla,rotr_pla_sun)
// Search for repeat orbits
//
// Calling Sequence
// params = CL_op_searchRepeatOrbits(smaMin,smaMax,QMin,QMax,ecc,sso [,incInput,er,mu,j2,rotr_pla,rotr_pla_sun])
//
// Description
// <itemizedlist><listitem>
// <p>Computes all repeat (also called "phased") orbits for a given range of semi major axis and number of planet revolutions per cycle.</p>
// <p>An orbit is considered as "phased" if the ground tracks (with respect to the rotating planet) repeat exactly, 
// which means that after a certain whole number of orbits, the planet must has rotated a whole number of revolutions with respect to the orbit's node.</p>
// <p>Let Tc be the duration of the repeat cycle. There must exist 2 integers M and Q such that:</p> 
// <p>Tc = M * T = Q * Tr </p>
// <p>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 (Tr is called "orbital day" in the following). </p>
// <p>If M is expressed as: M=N*Q+P (with P < Q), the problem comes down to determining the 3 integers N,P,Q such that: </p>
// <p>Tc = (N*Q+P)*T = Q * Tr (with P < Q)) </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>Two orbit cases are handled: </p>
// <p>- Orbits with a fixed inclination (incInput). There is only one unknown (sma), 
// which is determined using the equation: (N*Q+P)*T(sma) = Q * Tr(sma) </p>
// <p>- Sun-synchronous orbits: the inclination is computed in the process. 
// There are 2 unknowns (sma and inc). Two equations are then used: (N*Q+P)*T(sma) = Q * Tr(sma), and sun-sync(sma,inc) (condition for Sun-synchronicity). 
// For a Sun-synchronous orbit, an orbital day is the same as a usual day (86400 s).</p>
// <p>Also note that the computation uses the mean drifts, so that the true ground tracks may not repeat exactly if the argument of perigee drifts. </p>
// <p></p></listitem>
// <listitem>
// <p><b>params</b> is a 11-column matrix with as many rows as there are solutions.</p>
// <p>- params(:,1): Semi major axis [m]</p>
// <p>- params(:,2): Eccentricity</p>
// <p>- params(:,3): Inclination [rad]</p>
// <p>- params(:,4:6): N, P, Q </p>
// <p>- params(:,7): Number of orbits per cycle (=N*Q+P)</p>
// <p>- params(:,8): (mean) Nodal period [s]</p>
// <p>- params(:,9): Cycle duration [days]</p>
// <p>- params(:,10): Longitude gap between 2 closest ascending ground tracks at the equator [rad] = 2*pi/(N*Q+P)</p>
// <p>- params(:,11): Gap in longitude after one cycle modulo 2*pi [rad] (for checking purposes)</p>
// </listitem></itemizedlist>
//
// Parameters
// smaMin: Lower bound of semi major axis range (m) (1 x 1)
// smaMax: Upper bound of semi major axis range (m) (1 x 1)
// QMin: Minimum number of "orbital days" per repeat cycle (integer) (1 x 1)
// QMax: Maximum number of "orbital days" per repeat cycle (integer) (1 x 1)
// ecc : eccentricity (1 x 1)
// sso: 1 = search for sun-synchronous orbits, 0 = search for orbits with fixed inclination (1 x 1)
// incInput : (mandatory if sso==0, else optional) Inclination (rad) (1 x 1)
// 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) (1 x 1)
// rotr_pla_sun : (optional) mean apparent rotation rate of the Sun around the planet (default is %CL_rotrBodySun) (1 x 1)
// params : results (See description)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_op_repeatGroundTracks
// CL_op_repeat2smaInc
// CL_op_ssoJ2
//
// Examples
// smaMin = 7000.0 * 1000;
// smaMax = 7050.0 * 1000;
// QMin = 2;
// QMax = 6;
// ecc = 0.001;
// sso = 1;
// [params]=CL_op_searchRepeatOrbits(smaMin,smaMax,QMin,QMax,ecc,sso)
//
// sso = 0;
// incInput = CL_deg2rad(98.5);
// [params]=CL_op_searchRepeatOrbits(smaMin,smaMax,QMin,QMax,ecc,sso,incInput)


  function [N,P,Q] = CL__findNPQ(xNmin,xNmax,QMin,QMax)
  // Determines all (N,P,Q) such that: QMin<=Q<=QMax, 0<=P<Q, xNmin<=N+P/Q<=xNmax
  //
  // Calling Sequence
  // [N,P,Q] = CL__findNPQ(xNmin,xNmax,QMin,QMax)
  //
  // Parameters
  // xNmin: Nombre d'orbite par jour minimum (reel)
  // xNmax: Nombre d'orbite par jour maximum (reel)
  // QMin: Nombre de jours de phasage minimum (entier)
  // QMax: Nombre de jours de phasage maximum (entier)
  // N : valeurs de N (1 x Nt)
  // P : valeurs de P (1 x Nt)
  // Q : valeurs de Q (1 x Nt)
  //
  // Description
  // <itemizedlist><listitem>
  // Cette fonction renvoie toutes les combinaisons de N, P et Q qui respectent les conditions :
  // <p> xNmin <= N+P/Q <= xNmax </p>
  // <p> QMin <= Q <= QMax </p>
  // <p> P et Q premiers entre eux </p>
  // <p> 0<=P<Q (P=0 n'est sorti que si Q=1)</p>
  // <p> Note : la taille des vecteurs N,P,Q en sortie n'est pas connue a l'avance.</p>
  // </listitem></itemizedlist>
  //
  // Authors
  // A. LAMY - CNES - DCT/SB/MS
  // G. AZEMA - Thales Services pour le CNES - DCT/SB/MS
  //
  // See also
  //
  // Examples
  // xNmin = 14.25;
  // xNmax = 15.75;
  // QMin = 3;
  // QMax = 6;
  // [N,P,Q]=CL__findNPQ(xNmin,xNmax,QMin,QMax)

    Nmin = floor(xNmin);
    Nmax = ceil(xNmax);

    // Recherche des couples P,Q tels que P et Q soient premiers entre eux :
    P1=[];
    Q1=[];

    for iQ = QMin:QMax
      tabP = 1:iQ-1;
      I = find(CL_gcd(iQ*ones(tabP),tabP) <= 1); 
      tabP = tabP(I);
      if (iQ == 1); tabP = [0, tabP]; end    // add P=0 if Q == 1
      P1 = [P1, tabP];
      Q1 = [Q1, iQ * ones(tabP)];
    end

    N=[];
    P=[];
    Q=[];

    for iN = Nmin:Nmax
      N=[N, iN * ones(P1)];
      P=[P, P1]; // NB : P1 and Q1: same size
      Q=[Q, Q1];
    end

    // Keep N,P,Q if N+P/Q in [xNmin, xNmax] 
    I = find(N+P./Q >= xNmin & N+P./Q <= xNmax);
    N = N(I);
    P = P(I);
    Q = Q(I);

  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
if (smaMin < 0); CL__error("smaMin out of range"); end
if (smaMax <= smaMin); CL__error("smaMax out of range"); end
if (ecc < 0 | ecc >= 1); CL__error("eccentricity out of range"); end
if (QMin < 1); CL__error("QMin out of range"); end
if (QMax < QMin); CL__error("QMax out of range"); end
if (sso <> 0 & sso <> 1); CL__error("sso: wrong value"); end

if (sso == 0)
  if ~exists('incInput','local') then
    CL__error('You must give incInput when sso = 0');
  end
  I = find(incInput < 0 | incInput > %pi+%eps);
  if ~isempty(I); CL__error("inclination out of range"); end
else 
  incInput = 0; // not used ! 
end 



// computes xNmin, xNmax (range of number of orbits per planet rev) 
// NOTE :  
// it is assumed that : 
// - the number of orbit / planet rev is a 
//   monotonous function of the semi major axis (for sso == 0 or 1)
// - max sma of a sun-synchronous orbit corresponds to inc= %pi 

tabsma = [smaMin,smaMax]; 
if (sso == 1) 
  smaMax_sso = CL_op_ssoJ2('a', ecc, %pi); 
  tabsma(2) = min(smaMax, smaMax_sso) 
  tabinc =  CL_op_ssoJ2('i', tabsma, ecc);
else
  tabinc = [incInput, incInput];
end

[tabpomp, tabgomp, tabanmp] = CL_op_driftJ2(tabsma,ecc,tabinc,er,mu,j2);
tabTnod = 2*%pi ./ (tabpomp + tabanmp);  // nodal period
tabTrev = 2*%pi ./ (rotr_pla - tabgomp); // planet rev. period / node
tabxN = tabTrev ./ tabTnod; 
xNmin = min(tabxN); 
xNmax = max(tabxN); 

// (N,P,Q) such that N+P/Q in [xNmin, xNmax] 
[N,P,Q] = CL__findNPQ(xNmin,xNmax,QMin,QMax);

params = [];

if (and([N <> [], P <> [], Q <> []]))
  // computed sma and inc for each (N,P,Q)
  [sma,inc] = CL_op_repeat2smaInc(N,P,Q,ecc,sso,incInput,er,mu,j2,rotr_pla,rotr_pla_sun);
  
  // final results
  ecc = ecc .* ones(sma);
  [pomp, gomp, anmp] = CL_op_driftJ2(sma,ecc,inc,er,mu,j2);
  nbOrbCycle = N.*Q + P;
  maille = 2*%pi ./ nbOrbCycle;
  periodeNodale = 2*%pi ./ (pomp + anmp); 
  dureeCycle = nbOrbCycle .* periodeNodale;
  ecartLongitudeCycle = CL_rMod(dureeCycle .* (rotr_pla - gomp), -%pi, %pi);
  
  params = [sma; ecc; inc; N; P; Q; nbOrbCycle; periodeNodale; dureeCycle/86400.0; maille ; ecartLongitudeCycle]';
  
  // %nan removed if any 
  I = find(isnan(params(:,1)) | isnan(params(:,3))); // sma or inc is %nan
  if (~isempty(I)); CL__warning("%nan found in result: removed"); end
  params(I,:) = []; 
end
  
endfunction