//  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 [interv] = CL_ev_stationVisibility(t, pos, stations, elevmin, tperiods, opts, er, obla)
// Geometrical visibility intervals
// 
// Calling Sequence
// [interv] = CL_ev_stationVisibility(t, pos, stations, elevmin [, tperiods, opts, er, obla])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the periods of time (start and end times) when a spacecraft is visible from 
// ground stations. </p> 
// <p></p>
// <p>A spacecraft is visible from a ground station when its elevation (angle above horizon) 
// is over a given threshold (<b>elevmin</b>). The station coordinates must be elliptical 
// (that is geodetic) coordinates. There can be one minimum elevation value for all the stations, 
// or one value per station. </p>
// <p>The results are computed in the simulation time intervals defined by <b>tperiods</b>. </p>
// <p>The intervals returned in <b>interv</b> define the visibility periods start and end times 
// when the spacecraft is in view of <b>at least one ground station</b>. </p> 
// <p></p>
// <p>The position of the spacecraft (<b>pos</b>) must be defined in the same coordinates system 
// as the stations and can be: </p> 
// <p>1) A 3xN array of positions (one position for each time t). The values at intermediate times 
// are obtained by Lagrange interpolation of order <b>ninterp</b> (using t and pos).</p>
// <p>2) A list containing the function that computes the positions. The list must have the form:  </p> 
// <p> pos = list(fun, args), where the function fun has the interface: position = fun(t, args). </p>
// <p></p>
// <p>Other options can be defined in <b>opts</b>, which is a structure (optionally) containing 
// the following fields: </p>
// <p>- <b>prec</b>: Computation accuracy on the elevation angle. Default value is 1.e-3 deg.</p>
// <p>- <b>dmin</b>: Minimum length of visibility intervals. Default value is 60 s.</p>
// <p>- <b>ninterp</b>: Number of points for interpolating the positions. Default value is 8.</p>
// <p>- <b>impr</b>: True (%t) to activate the improved detection mechanism. Default value is %t.</p>
// <p></p>
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p> 
// <p>- The time step should be well chosen so that interpolation errors are small enough. </p>
// <p>- The results may be degraded if the position or the velocity (position time derivative) 
// is not continuous. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// t: Times at which the spacecraft positions are defined. (1xN)
// pos: Position vector in cartesian coordinates [m] or list. See above for details. (3xN or list) 
// stations: [lon;lat;alt] Stations elliptical (geodetic) coordinates [rad,rad,m]. (3xM)
// elevmin: Stations minimum elevations [rad]. (1xM or 1x1)
// tperiods: (optional) Simulation periods (time intervals). Default is [-%inf; %inf] (2xQ)
// opts: (structure, optional) Additional options (see above). Default value = empty structure. 
// er: (optional) Planet equatorial radius [m]. Default is %CL_eqRad.(1x1)
// obla: (optional) Planet oblateness. Default is %CL_obla. (1x1)
// interv: Visibility start and end times: [tvisi_start; tvisi_end]. (2xP)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_gm_stationPointing
// CL_detectSign
//
// Examples
// // Orbital elements and time instants
// cjd0 = 20000; 
// kep0 = [1.e7; 0; 1; 0; 0; 0]; 
// cjd = cjd0 + (0 : 5/1440 : 1); // days
// 
// // Ground stations definition (ECF)
// stations = [0; 0; 0]; // lon = lat = alt = 0
// elevmin = 5 * %pi/180;
// 
// // Computation of spacecraft position
// args = struct(); 
// args.cjd0 = cjd0; 
// args.kep0 = kep0; 
// 
// function [pos_ecf] = computePosition(cjd, args)
//   pos_eci = CL_oe_kep2car(CL_ex_kepler(args.cjd0, args.kep0, cjd)); 
//   pos_ecf = CL_fr_convert("ECI", "ECF", cjd, pos_eci); 
// endfunction
//
// // Case 1: Data interpolation
// pos = computePosition(cjd, args); 
// interv = CL_ev_stationVisibility(cjd, pos, stations, elevmin)
//
// // Case 2: Direct computation (no interpolation)
// interv = CL_ev_stationVisibility(cjd, list(computePosition, args), stations, elevmin)


// Declarations:

// -------------------------------------------------
// Intermediate function for elevation computation (data interpolation) 
// args contains: t, pos, ninterp, stations, elevmin, er, obla
// z = max(elev - elevmin), size = 1xN
// NB: 
// - args.elevmin: as many as stations (1xP)
// - ind unused ! 
// -------------------------------------------------
function [z] = stavisi_fct1(t, ind, args)
  pos  = CL_interpLagrange(args.t, args.pos, t, args.ninterp);
  // elev: each row <=> one different station
  elev = CL_gm_stationPointing(args.stations, pos, "elev", er=args.er, obla=args.obla); 
  z = max(elev - args.elevmin' * ones(t), "r");
endfunction

// -------------------------------------------------
// Intermediate function for elevation computation (user-defined function) 
// args contains: f, a, stations, elevmin, er, obla
// z = max(elev - elevmin), size = 1xN
// NB: 
// - args.elevmin: as many as stations (1xP)
// - ind unused ! 
// -------------------------------------------------
function [z] = stavisi_fct2(t, ind, args)
  pos  = args.f(t, args.a);
  // elev: each row <=> one different station
  elev = CL_gm_stationPointing(args.stations, pos, "elev", er=args.er, obla=args.obla); 
  z = max(elev - args.elevmin' * ones(t), "r");
endfunction

// -------------------------------------------------
// Argument Checking
// Only the specific arguments are checked. Other checks done in called functions
// -------------------------------------------------
if (argn(2) < 4)
  CL__error("Invalid number of input arguments (at least 4 expected)"); 
end

if (~exists("tperiods", "local")); tperiods = [-%inf; %inf]; end
if (~exists("opts", "local")); opts = struct(); end
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
if (~exists("obla", "local")); obla = CL__dataGetEnv("obla"); end

// check that all field names in opts are correct
Names = ["prec", "dmin", "ninterp", "impr"]; 
optsfields = fieldnames(opts); 
if (setdiff(optsfields, Names) <> [])
  CL__error("Invalid fields in opts structure"); 
end

if (~isfield(opts, "prec")); opts.prec = 1.e-3 * %pi/180; end  
if (~isfield(opts, "dmin")); opts.dmin = 60.0; end // sec
if (~isfield(opts, "ninterp")); opts.ninterp = 8; end
if (~isfield(opts, "impr")); opts.impr = %t; end

Nsta  = size(stations, 2);
Nelevmin = size(elevmin , 2);

if (size(stations,1) <> 3) 
  CL__error("Invalid size for argument stations"); 
end

if (size(elevmin,1) <> 1 | (Nelevmin <> 1 & Nelevmin <> Nsta)) 
  CL__error("Invalid size for argument elevmin"); 
end

if (size(t,1) <> 1)
  CL__error("Invalid size for argument t"); 
end

typeofpos = typeof(pos); 
if (typeofpos <> "constant" & typeofpos <> "list")
  CL__error("Invalid type for argument pos"); 
end

if (typeofpos == "constant")
  if (size(pos, 1) <> 3 | size(pos, 2) <> size(t, 2))
    CL__error("Invalid size for argument pos"); 
  end

else // if (typeofpos == "list")
  if (lstsize(pos) <> 2)
    CL__error("Invalid size for argument pos (list of 2 elements expected)"); 
  end
  if (typeof(pos(1)) <> "function")
    CL__error("Invalid values in argument pos (function expected)"); 
  end
end

if (opts.dmin < 0)
  CL__error("Invalid value for argument opts.dmin"); 
end

// ---------------------------------
// Computation 
// ---------------------------------
// Build arguments for intermediate functions
Args = struct();
Fct = []; // will be a function ! 

// common fields
Args.stations = stations;
Args.elevmin = elevmin .* ones(stations(1,:)); // size = 1xP, P = number of stations
Args.er = er;
Args.obla = obla; 

if (typeofpos == "constant")
  // for data interpolation
  Fct = stavisi_fct1; 
  Args.t = t; 
  Args.pos = pos; 
  Args.ninterp = opts.ninterp; 

else // if (typeofpos == "list")
  // for user-defined function 
  Fct = stavisi_fct2; 
  Args.f = pos(1); // function
  Args.a = pos(2); // 2nd argument

end


// visibility intervals
meth = "ds"; 
nsub = 0; // no intermediate points
isgn = 1; 

interv = CL__detectSign(t, [], Fct, Args, isgn, opts.prec, tperiods, nsub, meth, opts.impr);
			
// removes intervals of length <= dmin
if (interv <> [])
  I = find((interv(2,:) - interv(1,:)) * 86400 < opts.dmin); 			
  interv(:,I) = [];
end 
 
endfunction