//  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 [half_span,mid_pos] = CL_gm_betaEclipse(sma,betaa, ecc,argsp,er)
// Eclipse interval for circular or elliptical orbits.
//
// Calling Sequence
// [half_span,mid_pos] = CL_gm_betaEclipse(sma,betaa [,ecc,argsp,er])
//
// Description
// <itemizedlist><listitem>
// <p>Computes quantities that characterize the portion of the orbit 
// path where the Sun is eclipsed by the planet as seen from the satellite.</p> 
// <p> The shadow region is a cylinder whose axis is the (fixed) Sun direction and whose radius is that of the planet (assumed spherical). </p>
// <p> The computed quantities are: </p>
// <p> - <b>half_span</b>: eclipse interval half-length (in true anomaly) </p>
// <p> - <b>mid_pos</b>: anomaly of mid-eclipse, measured from the projection of the Sun direction onto the orbit plane. </p>
// <p> The true anomalies that correspond to the beginning and end of eclipse (resp. tan_beg and tan_end) can then be computed by: </p>
// <p> tan_beg = mid_pos - argsp - half_span </p>
// <p> tan_end = mid_pos - argsp + half_span </p>
// <p> where argsp is the angle between the projection of the Sun direction onto the orbit plane and the direction of the periapsis. </p>
// <p>If the Sun is not eclipsed, half_span is set to 0, and mid_pos to pi (arbitrarily). </p>
// <inlinemediaobject><imageobject><imagedata fileref="beta_eclipse.gif"/></imageobject></inlinemediaobject>
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p> 
// <p> - The trajectory of the satellite is an ellipse (no perturbations assumed). </p>
// <p> - If the orbit is circular (ecc = 0), the value of argsp has no impact on the results. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// sma: Semi major axis [m] (1xN or 1x1)
// betaa: Sun beta angle (can be positive or negative) [rad] (1xN or 1x1)
// ecc: (optional) Eccentricity. Default value is 0. (1xN or 1x1)
// argsp: (optional) Angle between the Sun direction and the periapsis in the orbit plane. Same sign convention as for the argument of periapsis. Default value is 0. (1xN or 1x1)
// er: (optional) Planet radius (default is %CL_eqRad) [m] (1xN or 1x1)
// half_span: Eclipse interval half-length (in true anomaly) [rad] (1xN)
// mid_pos: Angle from the Sun direction in the orbit plane defining the middle of the eclipse interval [rad] (1xN)
//
// See also
// CL_gm_eclipse
// CL_gm_raan2beta
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Example 1
// betaa = CL_deg2rad([10,20]);
// sma = 7.e6;
// CL_gm_betaEclipse(sma,betaa)
// CL_gm_betaEclipse(sma,betaa,0) // same
//
// // Example 2
// cjd = 20050;
// 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);
// raan = %pi/4;
// inc = CL_deg2rad(98.7);
// betaa = CL_gm_raan2beta(inc,raan,alpha_sun,delta_sun)
// sma = 10000.e3;
// ecc = [0, 0.01]; 
// argsp = 0;  
// [half_span,mid_pos] = CL_gm_betaEclipse(sma,betaa,ecc,argsp)

// Declarations:

// circular orbit
function [half_span, mid_pos] = eclipse_cir(sma, betaa, er)

  half_span = zeros(sma); 
  mid_pos = %pi * ones(sma); 

  K = sqrt(1 - (er./sma).^2); // sma assumed < er
  I = find(K < cos(betaa)); 
  half_span(I) = acos(K(I) ./ cos(betaa(I))); 

endfunction


// elliptical orbit (also works if ecc == 0)
function [half_span, mid_pos] = eclipse_ell(sma, betaa, ecc, argsp, er)

  // Solution to: 
  // 1 - cos(betaa)^2*cos(alpha)^2 = K*(1+e*cos(alpha-argsp))^2
  // t = tan(%pi/2 - alpha/2) 
  // => polynomial in t of degree 4
  // (alpha in [%pi/2, 3*%pi/2] => t in [-1,1]) 

  p = sma .* (1-ecc.^2); 
  K = (er ./ p).^2; 
  ex = ecc .* cos(argsp); 
  ey = ecc .* sin(argsp); 
  sb2 = sin(betaa).^2; 

  // polynomial coefficients (degree 0,1,...) 
  a0 = sb2 - K .* (1-ex).^2; 
  a1 = -4 * K .* ey .* (1-ex); 
  a2 = 2 * ( 2 - sb2 - K .* (1 + 2*ey.^2 - ex.^2) ); 
  a3 = -4 * K .* ey .* (1+ex); 
  a4 = sb2 - K .* (1+ex).^2;

  s = poly(0, "s"); 
  pols = a0 + a1*s^1 + a2*s^2 + a3*s^3 + a4*s^4; 

  N = size(sma,2); 
  t = %nan * ones(4,N);  // t: roots of polynomials (max: 4)

  // loop necessary because "roots" is not vectorized
  // NB : degree of polynomial may be less than 4
  for k = 1:N
    r = roots(pols(k)); // column vector
    t(1:length(r),k) = r; 
  end

  // select real roots belonging to [%pi/2, 3*%pi/2]
  // (%nan if condition not met)
  // note : "real(atan)" : "real" needed because t is complex
  I = find(imag(t) <> 0);
  t(find(imag(t) <> 0)) = %nan; 
  alpha = CL_rMod(%pi - 2 * real(atan(t)), 0, 2*%pi); 
  I = find(alpha <= %pi/2 | alpha >= 3*%pi/2); 
  alpha(I) = %nan; 

  // final results

  // nbr = number of roots <> %nan
  nbr = sum(bool2s(~isnan(alpha)), 'r'); 

  // check the number of roots (in case)
  if (find(nbr <> 0 & nbr <> 2) <> []) 
    CL__warning("CL_gm_betaEclipse: Number of roots not 0 nor 2"); 
  end

  alpha1 = %pi * ones(1:N); // beginning of interval by default
  alpha2 = alpha1; // end of interval by default  

  I = find(nbr > 0); 
  if (I <> []) 
    alpha1(I) = nanmin(alpha(:,I), 'r'); // min excluding %nan 
    alpha2(I) = nanmax(alpha(:,I), 'r'); // max excluding %nan 
  end

  half_span = (alpha2 - alpha1)/2; 
  mid_pos = (alpha1 + alpha2)/2;   
 
endfunction

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

// Code:
if (~exists('ecc','local')); ecc = 0; end
if (~exists('argsp','local')); argsp = 0; end
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end

s = [size(sma); size(betaa); size(ecc); size(argsp); size(er)];
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)
if (find(s(:,2) <> 1 & s(:,2) <> N) <> [])
  CL__error("Invalid argument sizes (number of columns)");
end

// adjust number of columns
if (s(1,2) < N); sma = sma * ones(1,N); end
if (s(2,2) < N); betaa = betaa * ones(1,N); end
if (s(3,2) < N); ecc = ecc * ones(1,N); end
if (s(4,2) < N); argsp = argsp * ones(1,N); end
if (s(5,2) < N); er = er * ones(1,N); end

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

// condition for abnormal inputs
cond_err = ecc < 0 | ecc >= 1 | sma.*(1-ecc) <= er; 

// (sufficient) condition for no eclipse (for efficiency)
cond_noecl = sma.*(1-ecc).*abs(sin(betaa)) >= er; 


// default results (no eclipse)

half_span = zeros(1:N);
mid_pos = %pi * ones(1:N);


// case: abnormal inputs => %nan

I = find(cond_err); 
half_span(I) = %nan; 
mid_pos(I) = %nan; 


// case: circular orbit (faster)

I = find(~cond_err & ~cond_noecl & ecc == 0); 
if (I <> [])
  [half_span(I), mid_pos(I)] = eclipse_cir(sma(I),betaa(I),er(I)); 
end


// case: elliptical orbit (general case)

I = find(~cond_err & ~cond_noecl & ecc <> 0); 
if (I <> [])
  [half_span(I), mid_pos(I)] = eclipse_ell(sma(I),betaa(I),ecc(I),argsp(I),er(I)); 
end


endfunction