//  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 [varargout] = CL_3b_manifold(env,orb,t_orb,period,epsilon,dtint,pars)
// Manifolds (divergent and convergent) for Halo and Lissajous orbits
//
// Calling Sequence
// [manifold1, ..., manifoldN] = CL_3b_manifold(env, orb, t_orb, period, epsilon, dtint, pars)
//
// Description
// <itemizedlist><listitem>
// <p>Computes <b>manifolds</b> for Halo or Lissajous orbits. </p>
// <p>It can compute all four branches specified by the value of <b>pars</b>: </p>
// <p>- pars = "conv": convergent, inward </p>
// <p>- pars = "-conv": convergent, outward </p>
// <p>- pars = "div": divergent, inward </p>
// <p>- pars = "-div": divergent, outward </p>
// <p>Manifolds in output are given in the same order as <b>pars</b>. </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>- A Halo or Lissajous orbit can be computed using <link linkend="CL_3b_halo">CL_3b_halo</link> or
// <link linkend="CL_3b_lissajous">CL_3b_lissajous</link>. </p>
// <p>- For the definition of adimensional quantities, see <link linkend="CL_3b_environment">CL_3b_environment</link>.</p> 
// <p>- In the literature it is said <b>epsilon</b> should be ~1.e-9, 
// but as the method is accurate enough, we recommend 1.e-5. </p> 
// </listitem>
// </itemizedlist>
//
// Parameters
// env: (struct) Lagrangian point structure 
// orb: Lissajous or Halo orbit (adimensional position and velocity) (6xN)
// t_orb: Time instants at which the orbit is given (adimensional) (1xN)
// epsilon: Tolerance. Recommended value: 1.e-5
// period: Period used to estimate the monodromy. It corresponds to omegahalo for the halo 
// orbits and nu for the Lissajous orbits 
// dtint: Relative (to t_orb) output times (>=0) (1xP) (adimensional) 
// pars: (string) Type of manifold to be computed: "conv", "-conv", "div", "-div" (1xP)
// manifolds: Generated manifolds (6 x P x N): 6 corresponds to position and velocity, P is the size of each trajectory making up the manifold, N is the number of trajectories)
//
// See also
// CL_3b_environment
// CL_3b_lissajous
// CL_3b_halo
//
// Authors
// CNES - DCT/SB
//
// Examples
// // Example with an Halo orbit:
// // Build structure for L2 and "Sun-EarthMoon" system:
// env = CL_3b_environment("S-EM","L2");
//
// // Generate a Halo orbit: 
// Az = 150.e6 / env.D;  // => adimensional
// direction = "pro";
// t_orb = linspace(0,180*86400,50) * env.OMEGA; // => adimensional
// [orb,omega] = CL_3b_halo(env, Az, direction, t_orb);
//
// // Compute manifolds:
// epsilon = 1.e-5;
// dtint = (0:0.5:200)*86400 * env.OMEGA; // adimensional
// conv_out = CL_3b_manifold(env, orb, t_orb, omega, epsilon, dtint, "-conv");
// 
// // Plot orbit and manifolds
// scf(); 
// for i = 1 : size(conv_out,3)
//   param3d(conv_out(1,:,i), conv_out(2,:,i), conv_out(3,:,i));
// end
// param3d(orb(1,:), orb(2,:), orb(3,:));
// h = CL_g_select(gce()); 
// h.foreground = 5; 

// 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:

// --- error checking
Norb = size(orb,2);
Nt_orb = size(t_orb,2);
if (Norb <> Nt_orb)
  CL__error("Invalid size for orb or t_orb"); 
end

tab_par = ["conv", "-conv", "div", "-div"]; 
if (setdiff(pars, tab_par) <> [])
  CL__error("Invalid types of manifolds (pars)");
end

Np = size(pars,"*");
lhs = argn(1);
if (lhs <> Np)
  CL__error("Invalid number of output arguments");
end

orb = [orb ; t_orb];  

// --- computation

// Computation of convergent and divergent directions
[orb_dir] = CL__3b_dirDirConv(orb, period/2, env.MU);

if (orb_dir == [])
  CL__error("Impossible to compute manifold"); 
end


// epsillon and type in the order: "conv", "-conv", "div", "-div" 
tab_eps = epsilon * [1, -1, 1, -1]; 
tab_typ = [1, 1, 0, 0]; 

// Manifolds generation
manifolds = list();

// loop on manifold types
for i = 1 : Np
  par = pars(i);
  ipar = find(par == tab_par); 

  // generate manifold and add to list
  manifold = CL__3b_traceManifold(orb_dir, tab_eps(ipar), dtint, env.MU, tab_typ(ipar));
  manifolds($+1) = manifold;
end


// --- output arguments
varargout = manifolds;


endfunction