//  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 [pso1,pso2,inters] = CL_gm_intersectPlanes(inc1,raan1,inc2,raan2)
// Intersection of 2 orbit planes
//
// Calling Sequence
// [pso1,pso2,inters] = CL_gm_intersectPlanes(inc1,raan1,inc2,raan2)
//
// Description
// <itemizedlist><listitem>
// <p>Computes the argument of latitude at the intersection of 2 orbit planes.</p> 
// <p>The orbit planes are defined by their inclinations and right ascensions of ascending node </p>
// <p></p></listitem>
// <listitem>
// <p><b>Notes:</b></p>
// <p> - Only one solution is computed (such that pso1 and pso2 
// belong to [0, pi]). The arguments of latitude for the second solution are pso1+pi and pso2+pi respectively.</p>
// <p> - If the orbits are circular, then <b>pso1</b> and 
// <b>pso2</b> define the positions in the orbits where the orbit paths intersect. </p>
// <p> - If the planes are identical, the number of intersections is infinite. 
// The output argument (<b>inters</b>) is then set to 0, otherwise it is equal to 1.</p>
// </listitem>
// </itemizedlist>
//
// Parameters
// inc1: Inclination of orbit 1 [rad] (1xN or 1x1)
// raan1: Right ascension of ascending node of orbit 1 [rad] (1xN or 1x1)
// inc2: Inclination of orbit 2 [rad] (1xN or 1x1)
// raan2: Right ascension of ascending node of orbit 2 [rad] (1xN or 1x1)
// pso1: Argument of latitude in orbit 1 where the 2 planes intersect [rad] (1xN)
// pso2: Argument of latitude in orbit 2 where the 2 planes intersect [rad] (1xN)
// inters: Flag indicating if the planes intersect (inters = 1) or is there exists an infinity of intersections (inters=0) (1xN)
//
// Authors
// CNES - DCT/SB
//
// Examples
// inc1 = 1; // rad
// raan1 = 2;
// inc2 = 0.1;
// raan2 = 1.9;
// [pso1,pso2] = CL_gm_intersectPlanes(inc1,raan1,inc2,raan2)
//

// Declarations:


// Code:

// checks arguments sizes are OK / resizes
[inc1,raan1,inc2,raan2] = CL__checkInputs(inc1,1, raan1,1, inc2,1, raan2,1); 


// Angular momentum vectors
w1 = [sin(raan1).*sin(inc1); -cos(raan1).*sin(inc1); cos(inc1)];
w2 = [sin(raan2).*sin(inc2); -cos(raan2).*sin(inc2); cos(inc2)];

// w = +/- w1 ^w2, with:
// w(3,:) >= 0
// NB: w may be %nan if w1 // w2

// eps: threshold on ||w||
eps = 1.e-10; 

[w, nw] = CL_unitVector(CL_cross(w1,w2));
I = find(nw < eps); 
w(:,I) = %nan; 
I = find(w(3,:) < 0); 
w(:,I) = -w(:,I); 

cos_pso1 = cos(raan1) .* w(1,:) + sin(raan1) .* w(2,:);
cos_pso2 = cos(raan2) .* w(1,:) + sin(raan2) .* w(2,:);

pso1 = real(acos(cos_pso1));
pso2 = real(acos(cos_pso2));

inters = ones(pso1);

// remove %nan
I = find(isnan(CL_dot(w)));
inters(I) = 0;
pso1(I) = 0; 
pso2(I) = 0; 


endfunction