//  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 [elev] = CL_gm_stationElevation(pos_ter,stations, er,obla)
// Elevation of any object as seen from ground stations - DEPRECATED
//
// Calling Sequence
// elev = CL_gm_stationElevation(pos_ter,stations [,er,obla])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the elevations (<b>elev</b>) of any object as seen from 
// one or several locations (<b>stations</b>).</p>
// <p>The position of the object (<b>pos_ter</b>) is given in cartesian coordinates 
// in a frame tied to the planet (that is in the same frame as <b>stations</b>). </p>
// <p>See below for detailed examples</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// pos_ter: positions in the rotating frame, in cartesian coordinates [X;Y;Z] [m] (3xM)
// stations: stations positions in the same rotating frame, in elliptical (geodetic) coordinates [long,lat,alt] [rad,m] (3xN)
// er : (optional) planet equatorial radius (default is %CL_eqRad) [m] (1x1)
// obla : (optional) planet oblateness (default is %CL_obla) (1x1)
// elev: elevations from each stations and each object position. elev(row=i,column=j) is the elevation of object j from station i. [rad] (NxM)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_ev_visibilityExtr
//
// Examples
// // Secular J2 propagation (in ECI frame)
// cjd0 = 21915;
// pas = 10.0 / 86400.0;
// cjd = cjd0 : pas : cjd0+1;
// kep0 = [7070.e3 ; 0.001 ; CL_deg2rad([98;90;10;15])];
// kep = CL_ex_secularJ2(cjd0, kep0, cjd);
//
// // Conversion to terrestrial frame (ECF frame)
// [pos_car,vel_car] = CL_oe_kep2car(kep);
// pos_ter = CL_fr_convert("ECI", "ECF", cjd, pos_car);
//
// // Stations definition
// sta1 = [CL_deg2rad(2);CL_deg2rad(70);200]; // high latitude
// sta2 = [CL_deg2rad(20);CL_deg2rad(0);400]; // equator
// stations = [sta1,sta2];
//
// // Elevations computation
// [elev] = CL_gm_stationElevation(pos_ter,stations);
//
// // Elevation in function of time
// scf();
// plot2d((cjd-cjd0)*24,CL_rad2deg(elev(1,:)),2) //station 1
// plot2d((cjd-cjd0)*24,CL_rad2deg(elev(2,:)),3) //station 2
//
// // Visibility duration if station's mask is 5 degrees :
// ind_1 = find(elev(1,:) > CL_deg2rad(5));  //station 1
// ind_2 = find(elev(2,:) > CL_deg2rad(5));  //station 2
// dur_1_minutes = pas*length(ind_1)*1440.0
// dur_2_minutes = pas*length(ind_2)*1440.0
//

// Declarations:


// Code:

CL__warnDeprecated(); // deprecated function

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

N = size(stations , 'c');      // nombre de stations
M = size(pos_ter , 'c');   // nombre de satellites

// On veut elev = pi/2 - acos(  [(verticale locale) scalaire (pos_sat-pos_sta)] / norme((pos_sat-pos_sta))  )
// On va decomposer le produit scalaire en 2 partie
pos_sta = CL_co_ell2car(stations,er,obla);
vert_loc = CL_co_sph2car( [stations(1:2,:) ; ones(1,N)]) ; // la verticale locale normee
  
// premiere partie : (verticale locale) scalaire (pos_sat)
prod_scal_1 = vert_loc' * pos_ter ;       // (NxM)
  
// deuxieme partie : (verticale locale) scalaire (pos_sta)
prod_scal_2 = CL_dot(vert_loc,pos_sta)' * ones(1,M) ;   // (NxM)
  
// la norme vaut : norme^2 = norme(pos_sat)^2 + norme(pos_sta)^2 - 2 * pos_sat scalaire pos_sta
norme = sqrt( ones(N,1) * (CL_norm(pos_ter).^2) + (CL_norm(pos_sta).^2)' * ones(1,M) - 2 * pos_sta'*pos_ter);
  
norme(find(norme == 0)) = %nan;
  
elev = %pi/2 - real(acos((prod_scal_1 - prod_scal_2) ./ norme)); // real pour eviter les erreurs numeriques
  
endfunction