//  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_visibilityExtr(cjd,type_oe,mean_oe,stations,stations_masks,sim_period, visi_min,prec,propag_model,ut1_tref,er,mu,j1jn,obla,fun_ci2cf)
// Geometrical visibility calculation - DEPRECATED
//
// Calling Sequence
// [visi_dates] = CL_ev_visibilityExtr(cjd,type_oe,mean_oe,stations,stations_masks, sim_period, ..
//                [visi_min,prec,propag_model,ut1_tref,er,mu,j1jn,obla,fun_ci2cf])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p>Replacement function: <link linkend="CL_ev_stationVisibility">CL_ev_stationVisibility</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> - Beware that the default propagation model is "j2sec" (J2 secular effects only). </p>
// <p> - Visibility intervals with a duration less than <b>visi_min</b> are discarded. </p>
// <p> - 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></p></listitem>
// <listitem>
// <p>Advanced users can specify an external function (<b>fun_ci2cf</b>) that converts positions 
//    from the inertial frame to the body fixed frame.</p>
// <p>This is particularly useful when one wants to compute visibilities for a planet that is not the Earth.</p>
// <p>The calling sequence of the function must be the following:</p>
// <p><literal>pos_cf = fun_ci2cf(cjd_tref,pos_ci,ut1_tref)</literal></p>
// <p>- <literal>cjd_tref</literal>: Modified (1950.0) julian date (Time scale: TREF)</p>
// <p>- <literal>pos_ci</literal>: Position vector in the inertial frame.</p>
// <p>- <literal>pos_cf</literal>: Position vector in the (body) fixed frame.</p>
// <p>- <literal>ut1_tref</literal>: UT1-TREF [seconds]</p>
// <p></p>
// <p>If <b>fun_ci2cf</b> is omitted, the frame conversion from ECI to ECF is used. 
// (See <link linkend="Reference frames">Reference frames</link> for more details)</p>
// <p></p></listitem>
// <listitem>
// <p>Warning:</p>
// <p> - The default propagation model is "j2sec": osculating elements = mean elements. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// cjd: Modified (1950.0) Julian date (Time scale: TREF) (1x1)
// type_oe: (string) Type of orbital elements used: "kep" or "cir" (1x1)
// oe: Satellite's mean orbital elements at date cjd in inertial frame (6x1)
// stations: Stations positions in 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]). Modified (1950.0) Julian date (Time scale: TREF) (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). See <link linkend="CL_ex_propagate">CL_ex_propagate</link> for details. Default is "j2sec" (1x1)
// ut1_tref: (optional) UT1-TREF [seconds]. Default is %CL_UT1_TREF (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)
// j1jn: (optional) Vector of j1jn coefficients J1 to Jn (troncated to J6) to be used (default is %CL_j1jn(1:6)) (1xN)
// obla: (optional) Planet oblateness (default is %CL_obla) (1x1)
// fun_ci2cf: (optional) See description for more details.
// visi_dates: Visibility start and end times ([cjd_visi_start ; cjd_visi_end]). Modified (1950.0) Julian date (Time scale: TREF) (2xM)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_ev_visibilityEph
// CL_gm_stationPointing
// CL_ex_propagate
//
// Examples
// t0 = 21915;
// mean_kep0 = [7070.e3; 1.e-3; 1.7; 0.1; 0.2; 0.3];
//
// // 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_visibilityExtr(t0, "kep", 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_visibilityExtr(t0, "kep", 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:

// -------------------------------------------------
// Default function to convert positions from inertial to rotating frame
// -------------------------------------------------
function [pos_ECF] = fun_ci2cf_ref(cjd_tref, pos_ECI, ut1_tref)
  pos_ECF = CL_fr_convert("ECI", "ECF", cjd_tref, pos_ECI, ut1_tref=ut1_tref);
endfunction


// -------------------------------------------------
// computation of max(elev - elevmin)
// NB: ind,args: not used but necessary for CL_fsolveb
// WARNING : This function uses variables defined in calling functions
// -------------------------------------------------
function [z] = f_elev_extr(t,ind,args)

  // Positions in inertial frame
  osc = CL_ex_propagate(propag_model,type_oe,cjd,mean_oe,t,"o",er,mu,j1jn);
  posvel_sat_ECI = CL_oe_convert(type_oe, "pv", osc, mu);

  // Positions in rotating frame
  pos_sat_ECF = fun_ci2cf(t, posvel_sat_ECI(1:3,:), ut1_tref);

  // Compute elevations minus respective minimum values
  elev = CL_gm_stationPointing(stations, pos_sat_ECF, "elev", er, obla);
  elev =  elev - stations_masks' * ones(1,size(elev,2));
  z = max(elev,'r');

endfunction


// -------------------------------------------------
// computation of start/end visi times
// -------------------------------------------------
function [visi_dates] = calc_visi_extr(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 = CL_fsolveb(f, t(I), t(I+1), dxtol = eps, meth="s", y1=z(I), y2=z(I+1));
  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 = CL_fsolveb(f, t(I), t(I+1), dxtol = eps, meth="s", y1=z(I), y2=z(I+1));
  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 
// -------------------------------------------------
CL__warnDeprecated(); // deprecated function

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 = "j2sec"; end
if (~exists("ut1_tref", "local")); ut1_tref = CL__dataGetEnv("UT1_TREF"); end
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end
if (~exists("j1jn", "local")); j1jn = CL__dataGetEnv("j1jn", 1:6); end
if (~exists("obla", "local")); obla = CL__dataGetEnv("obla"); end
if (~exists('fun_ci2cf','local')); fun_ci2cf = fun_ci2cf_ref; end
   
// Check/resize number of stations and masks
[stations, stations_masks, Nsta] = CL__checkInputs(stations,3, stations_masks,1);

if (Nsta == 0) 
   CL__error("At least one station expected"); 
end
     
if (size(cjd,2) <> 1)
   CL__error("Invalid size for initial time (cjd)"); 
end

if (size(mean_oe,2) <> 1 | size(mean_oe,1) <> 6)
   CL__error("Invalid size for mean orbital elements (mean_kep)"); 
end

if (typeof(fun_ci2cf) <> "function")
  CL__error("Invalid type of argument (fun_ci2cf)");
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 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_extr(f_elev_extr, 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_extr(f_elev_extr, 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 duration <= minimum are discarded
I = find(visi_dates(2,:) - visi_dates(1,:) > pas);
visi_dates = visi_dates(:,I);


endfunction