//  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_visiParams(sat_radius,target_radius,type_par1,par1,type_par2)
// Visibility parameters (angles, distance, ...) for a spherical body
//
// Calling Sequence
// par2 = CL_gm_visiParams(sat_radius,target_radius,type_par1,par1,type_par2)
// [res1,..resN]=CL_gm_visiParams(sat_radius,target_radius,type_par1,par1,[type1,..,typeN])
// [result]=CL_gm_visiParams(sat_radius,target_radius,type_par1,par1,"all")
//
// Description
// <itemizedlist><listitem>
// <p>Computes various visibility parameters: </p>
// <p>- satellite's semi view-angle (<b>sat</b>): 
// angle between the descending vertical and the direction of the target location </p>
// <p>- satellite's elevation (<b>elev</b>): 
// elevation of the satellite as seen from the target location </p>
// <p>- satellite's incidence (<b>incid</b> = pi/2 - elev): 
// incidence of the target-satellite direction from the target location </p>
// <p>- distance (<b>dist</b>): 
// distance between the satellite and the target location </p>
// <p>- centre angle (<b>cen</b>): 
// angle between the (body centre -> satellite) direction and the (body centre -> target location) direction. </p>
// <p></p>
// <p>Given the distance from the body centre to the satellite (<b>sat_radius</b>),
// the distance from the body centre to the target (<b>target_radius</b>) 
// and one of the following parameters (<b>type_par1</b>):</p>
// <p>- type_par1 = 'sat' , par1 = satellite's semi view angle.</p>
// <p>- type_par1 = 'elev' , par1 = elevation from the target location.</p>
// <p>- type_par1 = 'incid' , par1 = incidence (=pi/2-elev) from the target location.</p>
// <p>- type_par1 = 'dist' , par1 = distance between the target location and the satellite.</p>
// <p>- type_par1 = 'cen' , par1 = centre angle: angle from the body centre between the target location and the satellite.</p>
// <p>the function computes <b>par2</b> whose type is defined by <b>type_par2</b>.</p> 
// <p><b>type_par2</b> can also be an array of strings (any of the 5 parameters above), 
// or 'all' and in that case the result is a structure whose fields are the 5 parameters above</p>
// <p><inlinemediaobject><imageobject><imagedata fileref="visiParams.gif"/></imageobject></inlinemediaobject></p>
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p>
// <p> - A spherical body is assumed.</p>
// <p> - The function works even if "target_radius" is greater then "sat_radius". </p>
// <p> - If the input parameter is 'sat', 'elev" or "incid", there may be two solutions. Only one  
// is computed, it corresponds to the smallest value of the centre angle. </p>
// </listitem></itemizedlist>
//
// Parameters
// sat_radius: Distance from the satellite to the body centre [m] (Px1 or PxN)
// target_radius: Distance from the target to the body centre [m] (Px1 or PxN)
// type_par1: (string) Type of input parameter 'par1'. It can be 'sat', 'elev', 'incid', 'dist', 'cen'
// par1: Satellite's semi view angle, elevation, indicence, distance or centre angle [rad,m] (Px1 or PxN)
// type_par2: (string) Type of output parameter: 'sat', 'elev', 'incid', 'dist', 'cen' or a vector containing any of them, or 'all'.  
// par2: output parameter(s). Structure if type_par2 == "all". [rad,m] (PxN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// eqRad = CL_dataGet("eqRad");
// sat_r = eqRad + 700.e3; // 700 km altitude
// target_r = eqRad + 0; // Ground
// 
// // Distance to Ground incidence: 
// [incid]=CL_gm_visiParams(sat_r,target_r,'dist',800.e3,'incid')
//
// // Satellite view angle to Ground elevation:
// [elev]=CL_gm_visiParams(sat_r,target_r,'sat',CL_deg2rad(10),'elev') 
//
// // Centre angle to distance:
// [dist]=CL_gm_visiParams(sat_r,target_r,'cen',CL_deg2rad(7),'dist') 
//
// // Ground incidence to satellite view angle and distance:
// [sat,dist]=CL_gm_visiParams(sat_r,target_r,'incid',CL_deg2rad(15),['sat','dist'])
//
// // Satellite view angle to everything:
// [result]=CL_gm_visiParams(sat_r,target_r,'sat',CL_deg2rad(37),'all');


// Declarations:


// Code:

[lhs,rhs] = argn();

if (rhs <> 5) 
   CL__error("This function requires 5 input arguments"); 
end

if (typeof(type_par2) <> "string") 
   CL__error('Wrong type for argument type_par2'); 
end

if ~(size(type_par2,2) == lhs | (type_par2 == 'all' & lhs == 1) ) 
   CL__error("Wrong number of output arguments"); 
end

// check sizes
n = size(sat_radius,1); 
[sat_radius, target_radius, par1] = CL__checkInputs(sat_radius,n, target_radius,n, par1,n); 

// Check that sat_r > 0, target_r > 0 (otherwise: error)  
I = find(sat_radius <= 0 | target_radius <= 0); 
if (I <> [])
  CL__error("Invalid input argument: sat_radius or target_radius"); 
end


// --- Step 1: calculation of center angle 
// NB: error raised only if input parameter out of definition domain 

if (type_par1 == "elev" | type_par1 == "incid") 
  // par1 = elevation or incidence
  
  if (type_par1 == "incid")
    elev = %pi/2 - par1; 
  else
    elev = par1; 
  end
  // %eps : as value may be approximate
  I = find(abs(elev) > %pi/2 + 2*%eps);  
  if (I <> []) 
    CL__error("Invalid value of input argument (elevation or incidence)"); 
  end
  cosang = (target_radius ./ sat_radius) .* cos(elev); 
  I = find(abs(cosang) > 1); 
  cosang(I) = %nan; 
  angcen = real(acos(cosang)) - elev; 
  I = find(target_radius >= sat_radius); 
  angcen(I) = -real(acos(cosang(I))) - elev(I); 


elseif (type_par1 == "sat")  
  // par1 = "satellite angle"

  angsat = par1; 
  I = find(angsat < 0 | angsat > %pi); 
  if (I <> []) 
    CL__error("Invalid value of input argument (sat. angle)"); 
  end
  cosang = (sat_radius ./ target_radius) .* sin(angsat); 
  I = find(abs(cosang) > 1); 
  cosang(I) = %nan; 
  angsit = real(acos(cosang)); 
  
  // 2 solutions, but only one computed (such that: elev > 0)
  angcen = %pi/2 - angsat - angsit
  I = find(target_radius >= sat_radius); 
  if (I <> []); angcen(I) = %pi/2 - angsat(I) + angsit(I); end

  
elseif (type_par1 == "dist") 
  // par1 = distance

  dist = par1; 
  I = find(dist < 0); 
  if (I <> [])
    CL__error("Invalid value of input argument (distance)"); 
  end
  cosang = (sat_radius.^2 + target_radius.^2 - dist.^2) ./ (2*sat_radius.*target_radius);
  I = find(abs(cosang) > 1); 
  cosang(I) = %nan; 
  angcen = acos(cosang); 

elseif (type_par1 == "cen")  
  // par1 = "center angle" (from body center)
  I = find(par1 < 0 | par1 > %pi+2*%eps); 
  if (I <> [])
    CL__error("Invalid value of input argument (center angle)"); 
  end

  angcen = par1; 

else
  CL__error('Unknown type_par1 value');
end

// --- Step 2: calculation of output parameters as function of centre angle
res = struct("sat", 0 , "elev", 0 , "incid", 0 , "dist", 0, "cen", 0);

res.elev = atan(cos(angcen) - target_radius./sat_radius, sin(angcen)); 
res.incid = %pi/2 - res.elev; 
res.sat = %pi/2 - angcen - res.elev;
res.dist = sqrt(sat_radius.^2 + target_radius.^2 - 2.*sat_radius.*target_radius.*cos(angcen));
res.cen = angcen;

// output
if (type_par2 == 'all') 
  varargout(1) = res;
else
  for k = 1:size(type_par2,2)
    varargout(k) = res(type_par2(k));
  end
end

endfunction