//  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_visibilityEph(t_eph,pos_eph,stations,stations_masks,sim_period, visi_min,prec,er,obla)
// Geometrical visibility calculation - DEPRECATED
//
// Calling Sequence
// [visi_dates] = CL_ev_visibilityEph(t_eph, pos_eph, stations, stations_masks ...
//                [, sim_period, visi_min, prec, er, obla])
//
// 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 interpolated using order 8 Lagrange interpolation (using t_eph and pos_eph). </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> - The times <b>t_eph</b> should be well chosen so that intepolation errors are small enough. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// t_eph: Time at which the satellite positions are defined; Must be in increasing order and all different. [days] (1xNeph)
// pos_eph: Satellite's cartesian coordinates relative to the frame where the ground stations are defined (i.e. rotating frame). (3xNeph) 
// stations: Ground stations elliptical (geodetic) coordinates [long;lat;alt]. [rad,rad,m] (3xNsta)
// stations_masks: Station minimum elevations (above which there can be visibility). [rad] (1xNsta or 1x1)
// sim_period: (optional) Simulation time interval ([t_sim_start; t_sim_end]). Default is [t_eph(1); t_eph($)]. (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 sec. [sec] (1x1)
// er : (optional) Planet equatorial radius. Default is %CL_eqRad. [m] (1x1)
// obla : (optional) Planet oblateness. Default is %CL_obla. (1x1)
// visi_dates: Visibility start and end times: [t_visi_start; t_visi_end]. (2xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_ev_visibilityExtr
//
// Examples
// T = 1; // days
// t_eph = linspace(0,T,100);
// pos_eph = [cos(2*%pi*t_eph/T); zeros(t_eph); sin(2*%pi*t_eph/T); ] * 7.e6;  
//
// // ground stations definition
// sta1 = [0;0;0]; // equator
// sta2 = [0;%pi/2;0]; // North pole
// stations = [sta1,sta2];
// stations_masks = [CL_deg2rad(10), CL_deg2rad(0)];
//
// // visibility computation
// [visi_dates] = CL_ev_visibilityEph(t_eph,pos_eph,stations, ..
//                         stations_masks);
//
// [visi_dates] = CL_ev_visibilityEph(t_eph,pos_eph,stations, ..
//                         stations_masks, obla=0);
//

// Declarations:


// Code:

// -------------------------------------------------
// computation of max(elev - elevmin)
// NB: ind,args: not used but necessary for CL_fsolveb
// Warning: Uses "global" variables
// -------------------------------------------------
function [z] = f_elev_eph(t,ind,args)

  // positions aux dates voulues
  pos = CL_interpLagrange(t_eph, pos_eph, t, n=8);

  // stations elevations recalees par rapport au masque station
  elev = CL_gm_stationPointing(stations, pos, "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_eph(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

  // Tableau de sortie :
  visi_dates = [ start_dates ; end_dates ]; 

endfunction

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

CL__warnDeprecated(); // deprecated function

if (~exists('sim_period','local')); sim_period = [t_eph(1), t_eph($)]; end
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("obla", "local")); obla = CL__dataGetEnv("obla"); end

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

Nmask = size(stations_masks , 'c'); // nombre de stations masks
if ~(Nsta == Nmask | Nmask == 1) 
   CL__error('Invalid size for stations_masks'); 
end

if (size(t_eph, 1) <> 1)
   CL__error('Invalid size for t_eph'); 
end

if (size(pos_eph,1) <> 3 | size(pos_eph,2) <> size(t_eph,2))
   CL__error('Invalid size for pos_eph'); 
end

if (sim_period(2) <= sim_period(1) | sim_period(1) < t_eph(1) | ..
    sim_period(2) > t_eph($))
   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


// MAIN 
tdeb = sim_period(1); 
tfin = sim_period(2); 
pas = visi_min/86400.0;  // time step in days
eps = prec/86400; // accuracy in days

visi_dates = calc_visi_eph(f_elev_eph, tdeb, tfin, pas, eps); 

// Supprimer les visis de moins que le pas (visi_min)
I = find(visi_dates(2,:) - visi_dates(1,:) > pas);
visi_dates = visi_dates(:,I);


endfunction