//  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 [orb,omegahalo] = CL_3b_halo(env, Az, direction, t_orb)
// Halo orbit
//
// Calling Sequence
// [orb, omegahalo] = CL_3b_halo(env, Az, direction, t_orb)
//
// Description
// <itemizedlist><listitem>
// <p>Computes a halo orbit around a Lagrangian point.</p>
// <p></p>
// <p>The approximate amplitude along the Z axis if given by <b>Az</b>. </p>
// <p>The motion around +X axis can specified by <b>direction</b>: </p>
// <p>- direction = "pro": prograde motion around +X axis</p>
// <p>- direction = "retro": retrograde motion around +X axis</p>
// <p>The orbit is such that the Y component at the initial time (t_orb = 0) is 0. </p>
// <p></p></listitem>
// <listitem>
// <p><b>Notes: </b></p> 
// <p>- Before using this function, it is needed to create an "environment" (<b>env</b>) for 
// the chosen libration point and the chosen system (see <link linkend="CL_3b_environment">CL_3b_environment</link>).</p> 
// <p>- For the definition of adimensional quantities, see <link linkend="CL_3b_environment">CL_3b_environment</link>.</p> 
 // </listitem>
// </itemizedlist>
//
// Parameters
// env: (struct) Lagrangian point structure 
// Az: Estimate of amplitude along the Z axis (adimensional)
// direction: (string) Motion direction around +X axis: "pro" = prograde, "retro" = retrograde 
// t_orb: Time instants at which the orbit is computed (adimensional) (1xN)
// orb: Position and velocity [x;y;z;vx;vy;vz] defining the orbit (adimensional) (6xN)
// omegahalo: Orbit's angular rate (adimensional)
//
// See also
// CL_3b_environment
// CL_3b_lissajous
//
// Authors
// CNES - DCT/SB
//
// Examples
// env = CL_3b_environment("S-EM", "L2");
// Az = 150.e6 / env.D; // adimensional
// direction = "pro";
// t_orb = linspace(0,360*86400,100) * env.OMEGA; // 360 days -> adimensional
// [orb,omega] = CL_3b_halo(env,Az,direction,t_orb);
//
// // Plot orbit (normalized coordinates)
// scf(); 
// param3d((orb(1,:)-env.gl), orb(2,:), orb(3,:)); 
// xtitle("Orbit (3D), origin = Lagrange point"); 

// Bibliography
// 1) Introduction au probleme a trois corps et dynamique linearisee autour des points de Lagrange, G. Collange, Note Technique CCT Mecanique Orbitale num.7, CNES 2006
// 2) Estimation numerique des varietes stables et instables des orbites quasi-periodiques de Lissajous autour des points d'Euler (Lagrange L1, L2, L3), R. Alacevich, CNES septembre 2006
// 3) Exploration numerique d'orbites homoclines et heteroclines autour de L1 et L2 dans le probleme restreint a trois corps, rapport de stage, A. Martinez Maida, DCT/SB/MO 2007.0029301, CNES septembre 2007
//

// Declarations:

// Code:

// Arguments checking

if (Az <= 0)
  CL__error("Invalid amplitude value"); 
end

if (and([typeof(direction) <> "string", typeof(direction) <> "constant"]) | size(direction,"*") <> 1)
  CL__error("Invalid direction type or size"); 
end

// Compatibility with old version
// 0 <=> "retro"
// 1 <=> "pro"
if (typeof(direction) == "constant")
  I = find(direction == [0, 1]);
  if (I == []); CL__error("Invalid direction value"); end
  dirnames = ["retro", "pro"];
  direction = dirnames(I);
end

I = find(direction == ["pro", "retro"]); 
if (I == []); CL__error("Invalid direction value"); end
VALS = [1, -1]; 
idirection = VALS(I); 

if (find(t_orb < 0) <> [])
  CL__error("Invalid value(s) for t_orb"); 
end

// Main Code

// Calculating Initial Conditions
[X0,omegacorr] = CL__3b_condInitHalo(Az, idirection, env);

// Adjusting every point
[correction, Xt, omegahalo] = CL__3b_genHaloMat(X0, omegacorr, t_orb(1), t_orb($), env.pas, env.MU);

// Generating orbit
orb = CL__3b_formatOrbits(Xt, t_orb, env.MU);


endfunction