//  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 [res] = CL_gm_eclipse(sma,ecc,inc,argp,raan,alpha_sun,delta_sun, er,mu)
// Analytical eclipse calculation for elliptical orbits
//
// Calling Sequence
// res = CL_gm_eclipse(sma,ecc,inc,argp,raan,alpha_sun,delta_sun [,er,mu])
//
// Description
// <itemizedlist><listitem>
// <p>Computes various results that characterize the portion of the orbit path where the satellite is 
// in the shadow of the planet. </p>
// <p>The eclipsed region is a half cylinder of diameter the planet's diameter, and axis the Sun direction. 
// The calculation is purely geometrical: the Sun direction is supposed constant with respect to the orbit plane. </p>
// <p><b>res</b> is a structure with the following fields :</p>
// <p><b>res.start</b> and <b>res.end</b> contain quantities that define 
// the eclipse start and end positions: 'pso': (true) argument of latitude (w+v), 'ra': right ascension, 'decl': declination, tlt: true local time. </p>
// <p><b>res.sun_orb</b> contains quantities that define the Sun's position 
// ('alpha' and 'delta': spherical coordinates) in a frame tied to the orbit.</p>
// <p>The frame tied to the orbit is defined as follows:</p>
// <p>X-axis: Towards the ascending node</p>
// <p>Z-axis: Parallel to (and same direction as) the angular momentum vector</p>
// <p>Y-axis: Such that the frame is direct. </p>
// <p><b>res.angle</b> is the total eclipse's length (res.end.pso - res.start.pso)</p>
// <p><b>res.duration</b> is the duration, computed assuming a keplerian motion. </p>
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p>
// <p> - The planet is assumed spherical</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// sma: Semi major axis [m] (1xN or 1x1)
// ecc: Eccentricity (1xN or 1x1)
// inc: Inclination [rad] (1xN or 1x1)
// argp: Argument of periapsis [rad] (1xN or 1x1)
// raan: Right ascension of ascending node [rad] (1xN or 1x1)
// alpha_sun: Sun right ascension [rad] (1xN or 1x1)
// delta_sun: Sun declination [rad] (1xN or 1x1)
// er: (optional) Equatorial radius [m] (default is %CL_eqRad)
// mu: (optional) Gravitational constant [m^3/s^2] (default value is %CL_mu)
// res: (structure) Various quantities that define the eclipse [rad,sec] (each field is 1xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_gm_betaEclipse
//
// Examples
// cjd = CL_dat_cal2cjd(2009,03,21,6,0,0);
// pos_sun = CL_eph_sun(cjd);
// pos_sun_sph = CL_co_car2sph(pos_sun);
// alpha_sun = pos_sun_sph(1);
// delta_sun = pos_sun_sph(2);
// sma = 8000.e3;
// ecc = 0.1; 
// inc = 1.8;
// argp = 0; 
// raan = 0;
// res = CL_gm_eclipse(sma,ecc,inc,argp,raan,alpha_sun,delta_sun);
// res.start
// res.end
// res.sun_orb
// res.angle
// res.duration


// Declarations:

// Code:
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

s = [size(sma); size(ecc); size(inc); size(argp); size(raan); size(alpha_sun); size(delta_sun)];
smin = min(s, "r"); // min nb of rows, min nb of columns
smax = max(s, "r"); // max nb of rows, max nb of columns
N = smax(2); // max number of columns
  
// Check columns (must be 1 or N)
I = find(s(:,2) <> 1 & s(:,2) <> N);
if (I <> [])
  CL__error("Invalid argument sizes (number of columns)");
end

if (s(1,2) < N); sma = sma * ones(1,N); end
if (s(2,2) < N); ecc = ecc * ones(1,N); end
if (s(3,2) < N); inc = inc * ones(1,N); end
if (s(4,2) < N); argp = argp * ones(1,N); end
if (s(5,2) < N); raan = raan * ones(1,N); end
if (s(6,2) < N); alpha_sun = alpha_sun * ones(1,N); end
if (s(7,2) < N); delta_sun = delta_sun * ones(1,N); end

// -------------------------------------------------------------

// Keplerian mean motion
mm = CL_kp_params('mm',sma,mu); 

// Sun in orbit's frame
[alpha_sun_orb,delta_sun_orb] = CL_gm_inertial2orbitSph(inc,raan,alpha_sun,delta_sun); 

// eclipse computation 
betaa = delta_sun_orb; // beta angle
argsp = argp - alpha_sun_orb; // Sun -> Periapsis
[half_span, mid_pos] = CL_gm_betaEclipse(sma,betaa,ecc,argsp,er); 

// Final ress
// NB:
// index 1 <-> beginning of eclipse
// index 2 <-> end of eclipse

// (true) PSO (from asc. node) start/end of eclipse
// Note: pso1 <= pso2
pso1 = alpha_sun_orb + mid_pos - half_span;
pso2 = alpha_sun_orb + mid_pos + half_span;

// Right ascension/declination start/end of eclipse
[asc1, decl1] = CL_gm_orbit2inertialSph(inc,raan,pso1,0);
[asc2, decl2] = CL_gm_orbit2inertialSph(inc,raan,pso2,0);

// (true) local time start/end of eclipse
tlt1 = CL_rMod((asc1-alpha_sun)*12/%pi+12,0,24); 
tlt2 = CL_rMod((asc2-alpha_sun)*12/%pi+12,0,24); 

// Total angle (true anomaly)
ang = 2*half_span; 

// Mean anomalies and duration
anm1 = CL_kp_v2M(ecc, pso1-argp); 
anm2 = CL_kp_v2M(ecc, pso2-argp); 
dursec = CL_rMod(anm2-anm1,0,2*%pi) ./ mm; // seconds


// -------------------------------------------------------------
res = struct();

res.start = struct();
res.start.pso = CL_rMod(pso1, 0, 2*%pi);
res.start.ra = asc1;
res.start.decl = decl1;
res.start.tlt = tlt1;

res.end = struct();
res.end.pso = CL_rMod(pso2, res.start.pso, res.start.pso + 2*%pi);
res.end.ra = asc2;
res.end.decl = decl2;
res.end.tlt = tlt2;

res.sun_orb = struct();
res.sun_orb.alpha = alpha_sun_orb;
res.sun_orb.delta = delta_sun_orb;

res.angle = ang;
res.duration = dursec;


endfunction