//  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 [alpha2,delta2] = CL_gm_orbit2inertialSph(inc,raan,alpha1,delta1)
// Spherical coordinates conversion from an orbit frame to an inertial frame
//
// Calling Sequence
// [alpha2,delta2] = CL_gm_orbit2inertialSph(inc,raan,alpha1,delta1)
//
// Description
// <itemizedlist><listitem>
// <p>Performs a conversion of spherical coordinates from an "orbit" reference frame to an inertial reference frame.</p> 
// <p>The "orbit" reference frame is defined as follows:</p>
// <p>- X-axis: Towards the ascending node</p>
// <p>- Z-axis: Parallel to the angular momentum vector</p>
// <p>- Y-axis: Such that the frame is direct. </p>
// <p><inlinemediaobject><imageobject><imagedata fileref="orbit_inertial_sph.gif"/></imageobject></inlinemediaobject></p>
// </listitem>
// </itemizedlist>
//
// Parameters
// inc: Orbit's inclination [rad] (1xN)
// raan: Orbit's right ascension of ascending node [rad] (1xN)
// alpha1: Right ascension in orbit frame [rad] (1xN)
// delta1: Declination in orbit frame [rad] (1xN)
// alpha2: Right ascension in inertial frame [rad] (1xN)
// delta2: Declination in inertial frame [rad] (1xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// inc = CL_deg2rad(98);
// raan = CL_deg2rad(15);
// alpha1 = CL_deg2rad(50);
// delta1 = CL_deg2rad(60);
// 
// [alpha2,delta2] = CL_gm_orbit2inertialSph(inc,raan,alpha1,delta1)


// Declarations:


// Code:

Ng = size(raan,2)
Ni = size(inc,2)
Na = size(alpha1,2)
Nd = size(delta1,2)
N = max(Ng,Ni,Na,Nd)
coherence = (Ng==1|Ng==N)&(Ni==1|Ni==N)&(Na==1|Na==N)&(Nd==1|Nd==N)
if ~coherence then CL__error('bad size of input arguments'); end
if N~=1
  if Ng==1 then raan=raan*ones(1,N); end
  if Ni==1 then inc=inc*ones(1,N); end
  if Na==1 then alpha1=alpha1*ones(1,N); end
  if Nd==1 then delta1=delta1*ones(1,N); end
end

cosi = cos(inc);
sini = sin(inc);
ca1 = cos(alpha1);
sa1 = sin(alpha1);
cd1 = cos(delta1);
sd1 = sin(delta1);
cg = cos(raan);
sg = sin(raan);


//delta2 = elevation on the orbit plan
sindelta2 = sa1.*cd1.*sini + cosi.*sd1 ;
sindelta2(find(sindelta2<-1)) = -1 ;
sindelta2(find(sindelta2>1)) = 1 ;
delta2 = asin(sindelta2) ;

//alpha2: angle on the orbit plan / node
cd2ca2 = cg.*ca1.*cd1 - sg.*cosi.*sa1.*cd1 + sg.*sini.*sd1 ; //cos delta2 cos alpha2
cd2sa2 = sg.*ca1.*cd1 + cg.*cosi.*sa1.*cd1 - cg.*sini.*sd1 ; //cos delta2 sin alpha2
alpha2 = zeros(cd2ca2) ;
alpha2(find(cd2ca2==0 & cd2sa2==0)) = 0 ;
ii = find(~(cd2ca2==0 & cd2sa2==0)) ;
alpha2(ii) = atan(cd2sa2(ii),cd2ca2(ii)) ;

endfunction