//  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 [varargout] = CL_gm_eclipseCheck(pos_obs, pos1, pos2, sr1, sr2, res)
// Characteristic eclipse quantities
//
// Calling Sequence
// ratio = CL_gm_eclipseCheck(pos_obs, pos1 [, pos2, sr1, sr2])
// [val1, val2, ...] = CL_gm_eclipseCheck(pos_obs, pos1, [, pos2, sr1, sr2, res]); with res = [name1, name2, ...]
//
// Description
// <itemizedlist><listitem>
// <p>Computes quantities that characterize the eclipse conditions: ratio of solid angle that is eclipsed, 
// or characteristic angles defining umbra or penumbra. </p>
// <p></p>
// <p> - <b>Eclipse ratio</b> (<b>"eclrat"</b>):</p>
// <p> It is the hidden solid angle fraction of a body (body1) eclipsed by the second body (body2), as seen 
// from the observer. </p>
// <p>- ratio = 1.0 : 100% of body1 is eclipsed (umbra),</p>
// <p>- ratio = 0.7 : 70% of body1 is eclipsed (penumbra),</p>
// <p>- ratio = 0.0 : 0% of body1 is eclipsed. </p>
// <p></p>
// <p><inlinemediaobject><imageobject><imagedata fileref="eclipse_check.gif"/></imageobject></inlinemediaobject></p>
// <p></p>
// <p> - <b>Eclipse angles</b> (<b>"angu"</b>,<b>"angp"</b>,<b>"angc"</b>):</p>
// <p> Characteristic angles (positive or negative) that define the type of event: </p>
// <p> - angu: umbra </p>
// <p> - angp: penumbra </p>
// <p> - angc: umbra or penumbra, the eclipse region being defined by a cylinder of radius that of the occulting body (i.e. body2). The radius of body1 is not used in this case. </p>
// <p></p>
// <p>All the angles are: </p>
// <p>- continuous quantities (of the input data),  </p>
// <p>- positive if the eclipse condition is met and negative otherwise. </p>
// <p></p>
// </listitem>
// <listitem>
// <p><b>Notes</b>: </p>
// <p>The bodies are assumed spherical. The spheres should have positive radii, and should not intersect. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// pos_obs: 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)
// res: (string, optional) Names of computed quantities: "eclrat", "angu", "angp", "angc". Default is "eclrat" (1xP)
// val1, val2...: Computed quantities [rad or -] (1xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_gm_inters3dConesSa
//
// Examples
// // Eclipse of the Sun by Earth from an observer in space:
// theta = linspace(2.88, 3.05, 1000); 
// pos_sat = 42164.e3 * [cos(theta); sin(theta); zeros(theta)];
// au = CL_dataGet("au");
// pos_sun = au * [1; 0; 0]; // Sun
// eclrat = CL_gm_eclipseCheck(pos_sat, pos_sun); 
// angp = CL_gm_eclipseCheck(pos_sat, pos_sun, res="angp");  
// angu = CL_gm_eclipseCheck(pos_sat, pos_sun, res="angu"); 
//
// // Plot the results
// scf();
// plot(theta*180/%pi, eclrat, "k"); 
// plot(theta*180/%pi, angp, "r"); 
// plot(theta*180/%pi, angu+1, "b"); 
// CL_g_legend(gca(), ["eclrat", "angp", "angu + 1"]); 

// Declarations: 

// --------------------------------------------------------
// Computation of eclipse ratio. 
// same meaning of arguments as in main function 
// (but here: have all same size)
// No arguments checking done. 
// --------------------------------------------------------
function [rat_ecl] = ecl_ratio(pos_obs, pos1, pos2, sr1, sr2)

  // Check that observer is outside both spheres :
  norm1 = CL_norm(pos_obs-pos1);
  norm2 = CL_norm(pos_obs-pos2);
  I = find(norm1 < sr1 | norm2 < sr2);
  if (I <> [])
    CL__error("Observer must be outside both spheres");   
  end;
  
  alpha1 = asin( sr1 ./ norm1 ); // apparent radius of sphere1
  alpha2 = asin( sr2 ./ norm2 ); // apparent radius of sphere2

  // separation angle between the two spheres (as seen from obs)
  alpha = CL_vectAngle( pos1 - pos_obs , pos2 - pos_obs ); 
  
  // Solid angle of body1 occulted by body2
  ang_inters = CL_gm_inters3dConesSa(alpha,alpha1,alpha2);
  
  // If body2 was behind body1 : body1 is not occulted
  I = find( norm2 > norm1 );
  ang_inters(I) = 0;
  
  // Solid angle of body1 : 
  ang_tot = 2*%pi * (1 - cos(alpha1));
 
  // ratio of body1 eclipsed
  rat_ecl = %nan * ones(ang_inters); 
  I = find(ang_tot > 0); 
  rat_ecl(I) = ang_inters(I) ./ ang_tot(I);

endfunction


// --------------------------------------------------------
// Computation of eclipse angles. 
// same meaning of arguments as in main function 
// (but here: have all same size)
// typ = "angu", "angp", "angc"
// No arguments checking done. 
// --------------------------------------------------------
function [ang] = ecl_angles(pos, pos1, pos2, sr1, sr2, typ)

  // u: pos1 -> pos2 (e.g. Earth -> Sun)
  [u, D] = CL_unitVector(pos1-pos2);

  // theta: angle between axis (pos1, pos2) and line defining
  // the eclipse condition (less than 90 deg)
  theta = zeros(D); 

  if (typ == "angp")
    C = CL_dMult((sr2 ./ (sr2+sr1)) .* D, u); 
    D2 = CL_norm(C); 
    theta = asin(sr2 ./ D2); 
  
  elseif (typ == "angu")
    // If radii are ~equal => same as "angc" (C at infinity)
    cond_ok = (abs(sr2-sr1) > 1.e-8 * sr2); 
    I = find(cond_ok);  
    C = CL_dMult((sr2(I) ./ (sr2(I)-sr1(I))) .* D(I), u(:,I)); 
    D2 = CL_norm(C); 
    theta(I) = asin(sr2(I) ./ D2); 
  end

  // phi: angle between line of bodies and pos2->pos
  phi = CL_vectAngle(pos - pos2, -u); 

  d = CL_norm(pos - pos2); 
  I = find(d == 0 | d <= sr2); 
  d(I) = %nan; 
 
  // half view angle of occulting body from observer
  alpha = asin(sr2 ./ d);

  if (typ == "angp" | typ == "angc")
    ang = phi - alpha - theta; 
  else
    ang = phi - alpha + theta; 
    I = find(sr2 > sr1); 
    ang(I) = phi(I) - alpha(I) - theta(I); 
  end

  // opposite sign
  ang = -ang; 

endfunction


// Code :
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("res", "local")); res = "eclrat"; end

if (typeof(res) <> "string" | size(res,"*") <> 1)
  CL__error("Invalid type or size for argument: res"); 
end
  
Names = ["eclrat", "angu", "angp", "angc"]; 
if (setdiff(res, Names) <> [])
  CL__error("Invalid value for argument: res"); 
end

[pos_obs, pos1, pos2, sr1, sr2, N] = CL__checkInputs(pos_obs,3, pos1,3, pos2,3, sr1,1 ,sr2,1);

// Check that radii are positive (or 0)
I = find(sr1 < 0 | sr2 < 0);
if (I <> [])
  CL__error("Radius of spheres must be positive"); 
end
  
// Same number of output args as there are names in "res"
[lhs,rhs]=argn();
  
if (lhs <> size(res, "*"))
  CL__error("Wrong number of output arguments"); 
end


// Check that spheres don't intersect
// I = indices for which there is no intersection
I = find(CL_dot(pos1-pos2) > (sr1+sr2).^2);


for (k = 1:size(res, "*"))

  typ = res(k); 
  val = %nan * ones(pos_obs(1,:)); 

  if (typ == "eclrat" & I <> [])
    val(I) = ecl_ratio(pos_obs(:,I), pos1(:,I), pos2(:,I), sr1(I), sr2(I)); 

  elseif (I <> [])
    val(I) = ecl_angles(pos_obs(:,I), pos1(:,I), pos2(:,I), sr1(I), sr2(I), typ); 

  end

  varargout(k) = val; 
end


endfunction