//  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_orbSphTriangle(inc, angname, angval, output)
// Angle conversion in "orbit" spherical triangle 
//
// Calling Sequence
// [val1, val2,...] = CL_gm_orbSphTriangle(inc, angname, angval, output=[name1, name2, ...])
// res = CL_gm_orbSphTriangle(inc, angname, angval [,output="all"])
//
// Description
// <itemizedlist><listitem>
// <p>Converts angles in the orbit spherical triangle. </p>
// <p>Knowing the orbit inclination and given one angle of the triangle that can be: 
// declination (or latitude), right ascension (measured from the orbit's ascending node), 
// argument of latitude, or bearing (angle between North and course directions, clockwise), 
// this function computes the other ones. </p>
// <p></p>
// <p>Angles are referred to by name: </p>
// <p>- <b>aol</b>: argument of latitude = angle in the orbit plane from the ascending node. </p>
// <p>- <b>decl</b>: declination (or latitude). Specifying the declination only is not sufficient to 
// determine the other angles. So one may specify if the argument of latitude is increasing 
// (decl_a) or decreasing (decl_d). </p>
// <p>- <b>rra</b>: relative right ascension = angle in the equator measured from the ascending node. </p>
// <p>- <b>bear</b>: dihedral angle between the North direction and the course direction (tangent to the 
// path on the unit sphere). Specifying the bearing only is not sufficient to determine the other 
// angles. So one may specify if the argument of latitude correspond to a position in the Northern 
// hemisphere (bear_n) or in the Southern hemisphere (bear_s). </p>
// <p><inlinemediaobject><imageobject><imagedata fileref="orb_sph_triangle.gif"/></imageobject></inlinemediaobject></p>
// <p>NB: </p>
// <p>- In the figure: "phi", "delta", "alpha", "A" stand for <b>aol</b>, <b>decl</b>, <b>rra</b>, <b>bear</b>, respectively. </p>
// <p></p>
// </listitem>
// <listitem>
// <p>Notes: </p>
// <p>If output is "all", all the results are gathered in a structure containing the fields: 
// "aol", "decl", "rra", "bear". </p>
// <p>- All output values are recomputed values (including the one given as input). </p>
// <p>- Output values are silently set to %nan if input values are not in the correct range, 
// but an error occurs if the value is not in the definition range: inclination not in  [0, pi], 
// declination not in [-pi/2, pi/2], etc... </p>
// <p>- Beware of the range of output angles: aol, rra, bear in [-pi, pi], decl in [-pi/2, pi/2]. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// inc: Inclination [rad] (1xN or 1x1)
// angname: (string) Name of the input angle among: "aol", "rra", "decl_a", "decl_d", "bear_n", "bear_s" (1x1)
// angval: Value of the input angle [rad] (1xN or 1x1)
// output: (string) Names of the output angles among: "decl", "rra", "aol", "bear" or "all". Default is "all". (1x1)
// val1, val2... : Output data values (1xN)
// res : Structure containing all the results if output is "all" (containing 1xN arrays)
//
// Authors
// CNES - DCT/SB
//
// Examples
// inc = 1;
// aol = linspace(0, %pi/2, 10);
// decl = CL_gm_orbSphTriangle(inc, "aol", aol, output=["decl"]); 
// res = CL_gm_orbSphTriangle(inc, "aol", aol, output="all"); 


// ---------------------------------
// internal functions
// ---------------------------------

// Asin = asin
// + checks that argument is in [-1,1], otherwise: %nan
function [y] = Asin(x)
  I = find(abs(x) > 1); 
  y = real(asin(x)); 
  y(I) = %nan; 
endfunction

// Acos = acos
// + checks that argument is in [-1,1], otherwise: %nan
function [y] = Acos(x)
  I = find(abs(x) > 1); 
  y = real(acos(x)); 
  y(I) = %nan; 
endfunction

// Atan = atan
// + returns 0 or nan (if ret_nan=%t) if x and y are close to 0
function [z] = Atan(y,x,ret_nan)
  eps = 1.e-15;
  z = atan(y,x);
  I = find(abs(x) < eps & abs(y) < eps); 
  if (ret_nan)
    z(I) = %nan; 
  else
    z(I) = 0; 
  end
endfunction

// Atan = atan
// + returns 0 if x and y are close to 0
function [z] = Atan2(y,x)
  eps = 1.e-15;
  z = atan(y,x);
  I = find(abs(x) < eps & abs(y) < eps); 
  z(I) = %nan; 
endfunction

// ---------------------------------
// Main function to compute results
// ---------------------------------
function [res] = orbSphTriangle_compute(inc, angname, angval)
  // In all cases: first computes the argument of latitude
  // Then computes all other angles 
  // And puts everything in the structure

  res = struct();  

  if (angname == "aol") 
    res.aol = angval; 

  elseif (angname == "rra")
    // Atan(0,0) => pso undefined => return %nan
    res.aol = Atan(sin(angval).*cos(inc), cos(angval).*(cos(inc)).^2, %t); 

  elseif (angname == "decl_a" ) 
    if (find(angval < -%pi/2 - 2*%eps | angval > %pi/2 + 2*%eps) <> [])
      CL__error("Declination out of range"); 
    end
    angval = max(-%pi/2, min(%pi/2,angval)); 
    res.aol = Asin(sin(angval) .* CL__csc(inc)); 

  elseif (angname == "decl_d" )
    if (find(angval < -%pi/2 - 2*%eps | angval > %pi/2 + 2*%eps) <> [])
      CL__error("Declination out of range"); 
    end
    angval = max(-%pi/2, min(%pi/2,angval)); 
    res.aol = %pi - Asin(sin(angval) .* CL__csc(inc)); 

  elseif (angname == "bear_n")
    I = find(cos(inc) .* sin(angval) < 0); // invalid range
    res.aol = Acos(cos(inc).*cos(angval).*CL__csc(inc).*CL__csc(angval)); 
    res.aol(I) = %nan; 

  elseif (angname == "bear_s") 
    I = find(cos(inc) .* sin(angval) < 0); // invalid range
    res.aol = -Acos(cos(inc).*cos(angval).*CL__csc(inc).*CL__csc(angval)); 
    res.aol(I) = %nan; 

  else
    CL__error("Invalid angle name"); 

  end


  // computes remaining angles
  // no modulo except for aol, for consistency (various formulas)
  // Atan(0,0) => rra and bear = 0 (only possible if inc = pi/2 and pso = +/-pi/2, 
  //                                => rra and bear can be 0 or pi)

  res.aol = CL_rMod(res.aol, -%pi, %pi);  
  res.decl = Asin(sin(inc).*sin(res.aol));
  res.rra = Atan(sin(res.aol).*cos(inc), cos(res.aol), %f);
  res.bear = Atan(cos(inc), sin(inc).*cos(res.rra).*cos(res.decl), %f); 
endfunction


// ---------------------------------
// Code:
// ---------------------------------
if (~exists("output", "local")); output = "all"; end

// Checks / updates sizes
[inc, angval] = CL__checkInputs(inc, 1, angval, 1); 

// Checks number of output arguments (not too many!)
nbargout = argn(1); 
if (nbargout > size(output,2))
  CL__error("Too many output arguments"); 
end

if (find(inc < -2*%eps | inc > %pi+2*%eps) <> [])
  CL__error("Inclination out of range"); 
end
inc = max(0, min(%pi,inc)); 


// Special code if input is []
if (angval == [])
  res = struct("aol", [], "decl", [], "rra", [], "bear", []);
  
// Main computation
else
  res = orbSphTriangle_compute(inc, angname, angval);
end


// Copy results in output arguments
// ("all" is considered as a regular name)
for (k = 1:size(output,2))
  name = output(k); 
  if (name == "all")
    varargout(k) = res;
  else
    varargout(k) = res(name);
  end
end


endfunction