//  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, omega, nu] = CL_3b_lissajous(env, Ax, phix, Az, phiz, epsilon, t_orb)
// Lissajous orbit
//
// Calling Sequence
// [orb, omega, nu] = CL_3b_lissajous(env, Ax, phix, Az, phiz, epsilon, t_orb)
//
// Description
// <itemizedlist><listitem>
// <p>Computes a Lissajous orbit around a Lagrangian point.</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 
// Ax: Amplitude in the "x" direction (adimensional)
// phix: Phase in the "x" direction [rad]
// Az: Amplitude in the "z" direction (adimensional)
// phiz: Phase in the "z" direction [rad]
// epsilon: Accuracy for convergence
// 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)
// omega: Orbit's angular rate in "xy" plane (adimensional)
// nu: Orbit's angular rate along "z" axis (adimensional)
//
// See also
// CL_3b_environment
// CL_3b_halo
//
// Authors
// CNES - DCT/SB
//
// Examples
// env = CL_3b_environment("S-EM", "L2");
// Ax = 30.e6 / env.D; // adimensional
// phix = 0;
// Az = 100.e6 / env.D; // adimensional
// phiz = 0;
// epsilon = 1.e-10;
// t_orb = linspace(0,360*86400,100) * env.OMEGA; // 360 days -> adimensional
// orb = CL_3b_lissajous(env, Ax, phix, Az, phiz, epsilon, 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
//

// calcul des conditions initiales lineaires

// Declarations:


// Code:

// Calculating Initial Conditions
points_per_orbit = 5; 

// In order to have enough accuracy on the first and last points 
// we must add/subtract pi/omega at the beginning/end
t0 = t_orb(1) - %pi/env.omega_init; 
tf = t_orb($) + %pi/env.omega_init; 
t_init = linspace(t0, tf, int(tf-t0/env.omega_init)*points_per_orbit);
[X0,omega,nu] = CL__3b_condInitLiss(Ax,phix,Az,phiz,env,t_init);

// Adjusting every point
Xt = CL__3b_shooting(X0, t_init, env.MU, epsilon);

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

endfunction