//  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 [res] = CL_op_repeatGroundTracks(N,P,Q, option)
// Longitudes of repeating ground tracks
//
// Calling Sequence
// res = CL_op_repeatGroundTracks(N,P,Q [,option])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the longitudes of the ground tracks (at a reference latitude) of a repeat orbit.</p>
// <p><b>res</b> contains :</p>
// <p> - first column : res(:,1) : orbit number</p>
// <p> - second column : res(:,2) : number of orbital days in the cycle</p>
// <p> - third column : res(:,3) : longitude gap between first ground track and current ground track (rad)</p>
// <p> - fourth column : res(:,4) : longitude gap between first ground track and current ground track expressed 
// in reference longitude gaps, where reference longitude gap = 2*%pi/(N*Q+P)</p>
// </listitem></itemizedlist>
//
// Parameters
// N: Value of N (1x1)
// P: Value of P (1x1)
// Q: Value of Q (1x1)
// option: (optional) 0 = all longitudes, 1 = successive longitudes in [0; DL[, 2 = successive longitudes in [-DL/2; DL/2[ such that gap with initial longitude is decreasing.  Default is 1. (1x1)
// res: see description (Mx4)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_op_searchRepeatOrbits
// CL_op_repeat2smaInc
//
// Examples
// // number of longitude reference gaps 
// // versus number of orbital days :
// N = 15;
// P = 3;
// Q = 16;
// res = CL_op_repeatGroundTracks(N,P,Q);
// plot(res(:,2),res(:,4),'s');
//
// // number of longitude reference gaps versus
// // number of orbital days (decreasing values) :
// res = CL_op_repeatGroundTracks(N,P,Q,option=2);
// plot(res(:,2),res(:,4),'s');

// Declarations:


// Code:

if (~exists('option','local')); option=1; end


n = N*Q+P;

DL = 2*%pi * Q / (N*Q+P); // intertrace
DL0 = 2*%pi / (N*Q+P);    // maille elementaire


tab_k = (0:n)';
// longitudes entre 0 et intertrace suivante = DL
if (option == 1)
   tab_l = CL_rMod(tab_k * DL, -DL0/2, 2*%pi-DL0/2);
   ind = find (tab_l < DL);
// longitudes entre -DL/2 et DL/2, seulement quand l'ecart a diminue.
elseif (option == 2)
   tab_l = CL_rMod(tab_k(1:$) * DL, -%pi, %pi);
   ind=[];
   ecart_max=DL;
   jj = 0;
   last_ind = 1 ;

   while(jj~=[])
    ii = find(abs(tab_l) <= ecart_max + DL0/2);
    jj = find(ii > last_ind);
    if(jj~=[])
        last_ind = ii(jj(1));
      ind = [ind, last_ind];
      ecart_max = abs(tab_l(last_ind));
      tab_tmp = tab_l(last_ind);
    end
   end

// toutes les longitudes
else // option == 0
   tab_l = CL_rMod(tab_k * DL, -DL0/2, 2*%pi-DL0/2);
   ind = find (tab_l < 2*%pi); // tous!
end

tab_ki = tab_k(ind);

res = [ tab_ki, tab_ki*Q/n, tab_l(ind), round(tab_l(ind)/DL0) ];

endfunction