//  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 [visi_dates] = CL_ev_visibility(cjd,mean_kep,stations,stations_masks,sim_period, visi_min,prec,propag_model,er,mu,zonals,obla)
// Geometrical visibility start and end times - DEPRECATED
//
// Calling Sequence
// [visi_dates] = CL_ev_visibility(cjd,mean_kep,stations,stations_masks, ...
//                sim_period, [visi_min,prec,propag_model,er,mu,zonals,obla])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p>Replacement function: <link linkend="CL_ev_visibilityExtr">CL_ev_visibilityExtr</link></p>
// <p></p></listitem>
// <listitem>
// <p>Computes the periods of time (start and end times) when a satellite is visible from a given set of ground stations.</p> 
// <p>The satellite is visible from a ground station when its elevation is over a given threshold (<b>stations_masks</b>). </p>
// <p>The satellite trajectory is computed using an analytical model (specified by <b>propag_model</b>). </p>
// <p>The results are computed in the simulation period defined by <b>sim_period</b>. </p>
// <p>The intervals returned in <b>visi_dates</b> define the visibility periods start and end times 
// of the satellite by <b>at least one ground station</b>. 
// It means that the visibility intervals for ground stations considered independently are concatenated. </p>
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p> 
// <p> - Visibility intervals with length less than <b>visi_min</b> are not computed. 
// This parameter is also used for the detection of visibility intervals. So that choosing a small value for visi_min will increase the computation time. </p>
// <p> - Setting planet oblateness (<b>obla</b>) to 0 (so that computation are 
// done assuming a spherical planet) results in faster computation. </p>
// <p> - The inertial reference frame where the orbit parameters are defined is necessarily Gamma50 (Veis) as CL_mod_SidTimeG50 is used to convert positions to the rotating (planet fixed) reference frame. </p>
// <p> - The orbit parameters are the "classical" orbital elements. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// cjd: Modified julian days from 1950.0 (TUC)  (1x1)
// mean_kep: Satellite's keplerian mean orbital elements at time cjd in Gamma50 (Veis) frame [sma;ecc;inc;pom;raan;anm] (6x1)
// stations: Stations positions in the rotating (planet fixed) reference frame in elliptical (geodetic) coordinates [long,lat,alt] [rad,rad,m] (3xN)
// stations_masks: Station minimum elevations (above which there can be visibility) [rad] (1xN or 1x1)
// sim_period: Simulation time interval ([cjd_start; cjd_end]) in modified julian days from 1950.0 (TUC) (2x1)
// visi_min: (optional) Minimum visibility duration (default is 60 seconds) [sec] (1x1)
// prec: (optional) Computation accuracy on start/end visibility times (default is 1 second) [sec] (1x1)
// propag_model: (optional) Propagation model: 'kep' for keplerian, 'j2' for secular J2, 'lyd' for lyddane or 'eh' for Eckstein Hechler (default is lyddane) (1x1)
// er : (optional) Planet equatorial radius (default is %CL_eqRad) [m] (1x1)
// mu: (optional) Gravitational constant [m^3/s^2] (default value is %CL_mu)
// zonals: (optional) Vector of zonals coefficients J1 to Jn (troncated to J5) to be used (default is %CL_j1jn(1:5)) (1 x N)
// obla : (optional) Planet oblateness (default is %CL_obla) (1x1)
// visi_dates: Visibility start and end times: [cjd_visi_start ; cjd_visi_end] in modified julian days from 1950.0 (TUC) (2xM)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_ev_visibilityExtr
// CL_ev_visibilityEph
//
// Examples
// t0 = 21915;
// mean_kep0 = [7070.e3 ; 0.001 ; CL_deg2rad([98;90;10;15])];
//
// // Definition of ground stations
// sta1 = [CL_deg2rad(2);CL_deg2rad(70);200]; // high latitude
// sta2 = [CL_deg2rad(20);CL_deg2rad(0);400]; // equator
// stations = [sta1,sta2];
// stations_masks = [ CL_deg2rad(10) , CL_deg2rad(2) ];
//
// sim_period = [21915 ; 21918 ];  // 3 days
//
// // Visibility computation
// [visi_dates] = CL_ev_visibility(t0, mean_kep0, stations, ..
//                         stations_masks, sim_period);
//
// // Plot visibility duration (mn) as function of time 
// scf();
// plot2d3(visi_dates(1,:) - t0, .. 
//        (visi_dates(2,:) - visi_dates(1,:)) * 1440, style=2); 
//
// // Same computation with obla=0 (faster)
// [visi_dates] = CL_ev_visibility(t0, mean_kep0, stations, ..
//                 stations_masks, sim_period, obla=0);
//
// // Plot visibility duration (mn) as function of time 
// plot2d3(visi_dates(1,:) - t0, .. 
//        (visi_dates(2,:) - visi_dates(1,:)) * 1440, style=5); 


// Declarations:

// Code:

CL__warnDeprecated(); // deprecated function

  // -------------------------------------------------
  // computation of max(elev - elevmin)
  // -------------------------------------------------
  function [z] = f_elev(t)

  if(propag_model =='lyd')
    [moy,osc] = CL_ex_lyddane(cjd,mean_kep,t,er,mu,zonals);
    [pos_sat_G50,vit_sat_G50] = CL_oe_kep2car(osc);
  elseif(propag_model =='eh')
    [pos,vel] = CL_oe_kep2car(mean_kep);
    mean_cir = CL_oe_car2cir(pos,vel);
    [moy_cir,osc_cir] = CL_ex_eckHech(cjd,mean_cir,t,er,mu,zonals);
    [pos_sat_G50,vit_sat_G50] = CL_oe_cir2car(osc_cir);
  elseif(propag_model =='j2')
    [osc] = CL_ex_secularJ2(cjd,mean_kep,t,er,mu,zonals(2));
    [pos_sat_G50,vit_sat_G50] = CL_oe_kep2car(osc);
  elseif(propag_model =='kep')
    [osc] = CL_ex_kepler(cjd,mean_kep,t,mu);
    [pos_sat_G50,vit_sat_G50] = CL_oe_kep2car(osc);
  else
    CL__error('propag_model unknown');
  end

  // positions => terretrial frame
  pos_sat_ter = CL_fr_G502ter(t,pos_sat_G50);

  // compute elevations minus respective minimum values
  elev = CL_gm_stationElevation(pos_sat_ter,stations,er,obla);
  elev =  elev - stations_masks' * ones(1,size(elev,2));
  z = max(elev,'r');

  endfunction


  // -------------------------------------------------
  // generic computation of zero(f) by secant method
  // f is such that y = f(t) (t and y can be 1-D vectors)
  // solution is in ]y1, y2]
  // eps: threshold on abscissa variation 
  // -------------------------------------------------
  function [t, y] = zero_sec(f, t1, y1, t2, y2, eps)

  itermax = 20;
  
  K = 1:length(t1); 
  t = zeros(K); 
  y = zeros(K);
  dt = ones(K) * %inf; 
  iter = 1; 

  while (iter <= itermax & K ~= []) 
    t(K) = (y2(K) .* t1(K) - y1(K) .* t2(K)) ./ (y2(K) - y1(K)); 
    y(K) = f(t(K)); 

    i = find(y(K) .* y2(K) >= 0);
    I = K(i); 
    dt(I) = abs(t2(I)-t(I));    
    t2(I) = t(I); 
    y2(I) = y(I); 

    i = find(y(K) .* y1(K) > 0);
    I = K(i); 
    dt(I) = abs(t1(I)-t(I));    
    t1(I) = t(I);
    y1(I) = y(I);

    K = find (dt > eps); 
    iter = iter + 1; 
  end

  if (K ~= []) 
     CL__error("No convergence in zero computation");  
  end

  endfunction


  // -------------------------------------------------
  // computation of start/end visi times
  // -------------------------------------------------
  function [visi_dates] = calc_visi(f,tdeb,tfin,pas,eps)
  // eps : accuracy on t

  n = round((tfin-tdeb)/pas) + 2; // n >= 2
  t = linspace(tdeb,tfin,n); 

  z = f(t); // elevation minus 'min elevation for visibility'

  // --- Start of visibility ---
  start_dates = [];
  I = find( z(1:$-1) < 0 & z(2:$) >= 0 ); // not selected if z(1) >= 0 => added later 
  if (~isempty(I))
     start_dates = zero_sec(f, t(I), z(I), t(I+1), z(I+1), eps); 
  end

  // visibility from the beginning 
  if (z(1) >= 0)
     start_dates = [ tdeb, start_dates ];
  end

  // --- End of visibility ---
  end_dates = [];
  I = find( z(1:$-1) >= 0 & z(2:$) < 0 ); // not selected if z($) >= 0 => added later
  if (~isempty(I))
    end_dates = zero_sec(f, t(I), z(I), t(I+1), z(I+1), eps); 
  end
  
  // visibility at the end
  if (z($) >= 0)
    end_dates = [ end_dates, tfin ];
  end

  // --- Controls ---
  if (length(start_dates) <> length(end_dates))
    CL__error("Problem with visibility algorithm..."); 
  end
  
  if (length(start_dates) <> 0)
    if (find(start_dates > end_dates) ~= [])
      CL__error("Problem with visibility algorithm..."); 
    end
  end

  // outputs : beginning/end intervals
  visi_dates = [ start_dates ; end_dates ]; 

  endfunction

  // -------------------------------------------------
  // MAIN 
  // -------------------------------------------------

  if (~exists('visi_min','local')); visi_min=60.0; end
  if (~exists('prec','local')); prec=1; end
  if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
  if (~exists('propag_model','local')); propag_model='lyd'; end
  if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end
  if (~exists("zonals", "local")); zonals = CL__dataGetEnv("j1jn", 1:5); end
  if (~exists("obla", "local")); obla = CL__dataGetEnv("obla"); end

  // number of ground stations
  Nsta = size(stations , 'c');      
  if (Nsta == 0) 
     CL__error('At least one station expected'); 
  end

  // number of ground station masks (min elevations)
  Nmask = size(stations_masks , 'c');      
  if ~(Nsta == Nmask | Nmask ==1) 
     CL__error('Invalid size for stations_masks'); 
  end

  // number de dates
  Mt = size(cjd , 'c');          
  if (Mt ~= 1)
     CL__error('Invalid size for initial time (cjd)'); 
  end

  // number of satellites
  Morb = size(mean_kep , 'c');   
  if (Morb ~= 1)
     CL__error('Invalid size for mean orbital elements (mean_kep)'); 
  end

  if (sim_period(2) <= sim_period(1))
     CL__error('Invalid simulation period'); 
  end

  if (visi_min <= 0)
     CL__error('Invalid value for visi_min'); 
  end

  if (Nmask == 1); stations_masks = stations_masks*ones(1,Nsta); end


  // If the number of propagation dates is too big, 
  // => loop over shorter intervals 
  // => computation of the maximal number of propagation dates (approximate value)
  sz = stacksize() ; 
  n_tab_max = 0.7*(sz(1) - sz(2));  
  
  // Propagation times
  tdeb = sim_period(1); 
  tfin = sim_period(2); 
  pas = visi_min/86400.0;  // time step in days
  eps = prec/86400; // accuracy in days

  n_tab_tot = max(round((tfin-tdeb)/pas)*3*Nsta*3 , round((tfin-tdeb)/pas)*100); 
  nb_boucles = floor(n_tab_tot / n_tab_max) + 1; // value >= 1
  
  if (nb_boucles == 1)
    visi_dates = calc_visi(f_elev, tdeb, tfin, pas, eps); 

  else 
    visi_dates = [];
    tk = linspace(tdeb, tfin, nb_boucles+1); // start/end of sub-segments
    
    for k = 1 : nb_boucles
  
      visi_dates1 = calc_visi(f_elev, tk(k), tk(k+1), pas, eps); 

      // merging of visibility intervals (if they overlap)
      if ( ~isempty(visi_dates) & ~isempty(visi_dates1))
         if(visi_dates(2,$) == visi_dates1(1,1))
            visi_dates(2,$) = visi_dates1(2,1); 
            visi_dates = [ visi_dates, visi_dates1(:,2:$) ];
         else
       visi_dates = [ visi_dates, visi_dates1 ];
     end
      else
         visi_dates = [ visi_dates, visi_dates1 ];
      end
    end

  end

  // intervals with length <= minimum are discarded
  I = find(visi_dates(2,:) - visi_dates(1,:) > pas);
  visi_dates = visi_dates(:,I);


endfunction