//  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_eclipse(t, pos, pos1, pos2, sr1, sr2, typ, tperiods, opts)
// Eclipse intervals calculation
// 
// Calling Sequence
// [interv] = CL_ev_eclipse(t, pos, pos1 [, pos2, sr1, sr2, typ, tperiods, opts]); 
// 
// Description
// <itemizedlist><listitem>
// <p>Computes eclipse entry and exit times. </p> 
// <p>The type of event is defined by <b>typ</b>: </p>
// <p>- <b>umb</b>: umbra, </p>
// <p>- <b>pen</b>: penumbra, </p>
// <p>- <b>umbc</b>: umbra, whereby the shadowed region is a cyclinder with radius equal to 
// that of the occulting body. </p>
// <p></p>
// <p>The positions are interpolated using order 8 Lagrange interpolation (using t and pos). </p>
// <p>The results are computed in the simulation time intervals defined by <b>tperiods</b>. </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>: Accuracy on the characteristic angle. Default value is 1.e-3 deg.</p>
// <p>- <b>dmin</b>: Minimum length of eclipse 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>The intervals returned in <b>interv</b> define the eclipse entry and exit times. </p>
// <p></p>
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p> 
// <p> - Eclipse intervals with length less than <b>dmin</b> are discarded. </p>
// <p> - The times <b>t</b> should be well chosen so that interpolation errors are small enough. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// t: Time at which the observer positions are defined. (1xN)
// pos: Position of observer [m]. (3xN or 3x1)
// pos1: Position of body1 (= potentially eclipsed body) [m]. (3xN or 3x1)
// pos2: (optional) Position of body2 (= occulting body) [m]. Default is [0;0;0]. (3xN or 3x1)
// sr1: (optional) Radius of body1 (sphere) [m]. Default is %CL_body.Sun.eqRad. (1x1)
// sr2: (optional) Radius of body2 (sphere) [m]. Default is %CL_eqRad. (1x1)
// typ: (string, optional) Type of event: "umb", "pen", "umbc". Default is "umbc". (1x1)
// tperiods: (optional) Simulation periods (time intervals). Default is [-%inf; %inf] (2xQ)
// opts: (structure, optional) Additional options (see above). Default value = empty structure. 
// interv: Eclipse entry and exit times: [t_start; t_end]. (2xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_gm_eclipseCheck
// CL_detectSign
//
// Examples
// // Start/end of eclipse periods (umbra) for an observer in space:
// t = 0 : 5/1440 : 3; 
// pos_sat = 42164.e3 * [cos(2*%pi*t); sin(2*%pi*t); zeros(t)];
// au = CL_dataGet("au");
// pos_sun = au * [1; 0; 0]; // Sun
// CL_ev_eclipse(t, pos_sat, pos_sun, typ="umb") 
//

// Declarations:

// -------------------------------------------------
// Computation of characteristic angle
// (eclipse <=> angle is positive)
// args must contain: t, p = [pos; pos1; pos2], ninterp, sr1, sr2, ang_type
// NB: ind: not used
// -------------------------------------------------
function [val] = ecl_fct(t, ind, args)
  // interpolates all positions at the same time for efficiency
  P = CL_interpLagrange(args.t, args.p, t, args.ninterp);
  // eclipse angle (change sign)
  val = CL_gm_eclipseCheck(P(1:3,:), P(4:6,:), P(7:9,:), args.sr1, args.sr2, args.ang_type); 
endfunction

// -------------------------------------------------
// Checks
// -------------------------------------------------
if (argn(2) < 3)
  CL__error("Invalid number of input arguments (at least 3 expected)"); 
end

if (~exists("pos2", "local")); pos2 = [0;0;0]; end
if (~exists("sr1", "local")); sr1 = CL__dataGetEnv(["body", "Sun", "eqRad"]); end
if (~exists("sr2", "local")); sr2 = CL__dataGetEnv("eqRad"); end
if (~exists("typ","local")); typ = "umbc"; end 
if (~exists("tperiods", "local")); tperiods = [-%inf; %inf]; end
if (~exists("opts", "local")); opts = struct(); 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

if (size(sr1, "*") <> 1 | size(sr2, "*") <> 1)
  CL__error("Invalid size for arguments sr1 or sr2"); 
end

// Check that radii are positive
if (sr1 < 0 | sr2 < 0)
  CL__error("Invalid values for arguments sr1 or sr2 (must be positive)"); 
end

// check size / resize
[pos, pos1, pos2, N] = CL__checkInputs(pos,3, pos1,3, pos2,3);

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

if (typeof(typ) <> "string" | size(typ, "*") <> 1)
  CL__error("Invalid type or size for argument typ"); 
end

Types = ["umb", "pen", "umbc"]; 
AngTypes = ["angu", "angp", "angc"]; 

I = find(typ == Types); 
if (I == [])
  CL__error("Invalid value for argument: typ"); 
end

// arg for eclipseCheck
ang_type = AngTypes(I); 

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

// -------------------------------------------------
// Computation
// -------------------------------------------------

// data for interval computation
args = struct(); 
args.t = t; 
args.p = [pos; pos1; pos2]; 
args.sr1 = sr1; 
args.sr2 = sr2; 
args.ang_type = ang_type; 
args.ninterp = opts.ninterp; // interpolation order

// interval computation
meth = "ds"; 
nsub = 0; // no intermediate points
isgn = 1; // positive

interv = CL__detectSign(t, [], ecl_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