//  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 [pos_sat] = CL_gm_stationVisiLocus(station, azim, elev, rsat, er, obla)
// Satellite position for visibility in specified local direction 
//
// Calling Sequence
// [pos_sat] = CL_gm_stationVisiLocus(station, azim, elev, rsat [, er, obla])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the position of the satellite on a sphere (with a given radius) given the direction 
// of the satellite in the local frame of some location. </p>
// <p></p>
// <p><b>station</b> is a location defined by its elliptical (geodetic) coordinates. </p> 
// <p><b>azim</b> (azimuth, posive towards West) and <b>elev</b> (elevation) define the direction 
// of the satellite in the local (topocentric North) frame. The corresponding unit vector in 
// spherical coordinates is [azim; elev; 1]. </p>
// <p><b>pos_sat</b> is the position on the sphere of radius <b>rsat</b> in cartesian coordinates. </p>
// <p></p>
// </listitem>
// <listitem>
// <p><b>Notes:</b></p>
// <p>This function can be used to compute "visibility circles". </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// station: Station position, in elliptical (geodetic) coordinates  [long,lat,alt] [rad,m] (3x1 or 3xN)
// azim: First spherical coordinate of the pointing direction in the topocentric North frame [rad] (1x1 or 1xN)
// elev: Second spherical coordinate of the pointing direction in the topocentric North frame [rad] (1x1 or 1xN)
// rsat: Satellite radius (= distance from body center) [m] (1x1 or 1xN)
// er : (optional) Equatorial radius (default is %CL_eqRad) [m] (1x1)
// obla : (optional) Oblateness (default is %CL_obla) (1x1)
// pos: Satellite position (cartesian coordinates) [m] (3xN)
//
// See also
// CL_gm_stationPointing
//
// Authors
// CNES - DCT/SB
//
// Examples
// station = [10 * %pi/180; 45 * %pi/180; 0]; 
// azim = linspace(0, 2*%pi, 100); 
// rsat = 8000.e3; 
//
// // Elevation = constant value
// elev = 10 * %pi/180; 
// pos_sat = CL_gm_stationVisiLocus(station, azim, elev, rsat); 
//
// // Elevation = f(azimuth)
// elev = 10 * (%pi/180) * (1+sin(4*azim)); 
// pos_sat2 = CL_gm_stationVisiLocus(station, azim, elev, rsat); 
//
// scf(); 
// CL_plot_earthMap(color_id=color("grey60"));  
// CL_plot_ephem(pos_sat, color_id=2); 
// CL_plot_ephem(pos_sat2, color_id=5); 

// Declarations:

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


// check station separately (no unnecesary resizing)  
[azim,elev,rsat,N] = CL__checkInputs(azim,1,elev,1,rsat,1); 
[station,N2] = CL__checkInputs(station,3); 

if (N2 <> 1 & N2 <> N)
  CL__error("Invalid size for ''station''"); 
end 


// cartesian coordinates of "station"
pos_sta = CL_co_ell2car(station, er=er, obla=obla); 

// transf matrix to topoN frame (Z=vertical, X=north)
// (X_topo_frame = mat_topoN * X_initial_frame)
mat_topoN = CL_fr_topoNMat(station); 

// direction in initial frame
u = mat_topoN' * CL_co_sph2car([azim; elev; ones(azim)]); 

// intersection with sphere of radius "rsat"
pos_sat = CL_gm_inters3dLineEllips(pos_sta, u, er=rsat, obla=0); 


endfunction