//  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_inertial2orbitSph(inc,raan,alpha1,delta1)
// Spherical coordinates conversion from an inertial frame to an orbit frame
//
// Calling Sequence
// [alpha2,delta2] = CL_gm_inertial2orbitSph(inc,raan,alpha1,delta1)
//
// Description
// <itemizedlist><listitem>
// <p>Performs a conversion of spherical coordinates from an inertial reference frame to an "orbit" 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 inertial frame [rad] (1xN)
// delta1: Declination in inertial frame [rad] (1xN)
// alpha2: Right ascension in orbit frame [rad] (1xN)
// delta2: Declination in orbit 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_inertial2orbitSph(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)
cd1 = cos(delta1)
sd1 = sin(delta1)
a1g = alpha1 - raan

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

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

endfunction