//  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_manifolds(env,orb,t_orb,period,epsilon,tint,pars)
// Manifolds (divergent and convergent) from Halo and Lissajous - DEPRECATED
//
// Calling Sequence
// [manifold1, ..., manifoldN] = CL_3b_manifolds(env, orb, t_orb, epsilon, tint, pars)
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p>Replacement function: <link linkend="CL_3b_manifold">CL_3b_manifold</link></p>
// <p></p></listitem>
// <listitem>
// <p>Computes <b>manifolds</b> for Halo or Lissajous orbits. </p>
// <p>It can compute all four branches depending on 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 then 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 
// tint: Integration time (adimensional) 
// pars: (string) Type of manifold to be computed: "conv", "-conv", "div", "-div" (1xP)
// manifolds: Generated manifolds (6 x n x nb_points,  6 corresponds to position and velocity, 
// n is the size of each trajectory making up the manifold, nb_points is the number of trajectories)
//
// 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
//
// 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 = 0;
// 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;
// tint = 150*86400 * env.OMEGA // adimensional
// [conv_out, div_out] =  ..
// CL_3b_manifolds(env, orb, t_orb, omega, epsilon, tint, ["-conv","-div"]);
// 
// // Plot orbit and manifolds (conv_out only)
// 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; 
//

// Auxiliary function
// NB: KEPT HERE FOR ASCENDING COMPATIBILITY PURPOSES

function [Traj_trace] = manifolds_trace(Traj_dir,mod,duree,MU,Traj_type)
  //This function calculates all the trajectories of the manifolds defined by de Traj_dir(1:6,:)
  //For each point of the orbit. The point is perturbated by mod*Traj_dir(8:13,:)
  //or mod*Traj_dir(14:18,:), it depends on if it's convergent or divergent
  //
  
  //Imputs: Traj_dir:      Matrix. Dimension (19,n=number of points of the orbit)
  //                        position, velocity and date,(1:7)
  //                        divergent direction on position and velocity (8:13)
  //                        Convergent direction on position and velocity (14:18)
  //                        Traj_dir is usually generated by manifolds_dirconv.sci
  //          mod:          perturbation module,(POSITIF ou NEGATIF)(real)
  //      duree:        Extrapolation duration for each manifold (real POSITIVE)
  //          Traj_type:    0 Convergent manifold
  //                        1 Divergent manifold
  //
  //Outputs: Traj_trace:    HyperMatrix. dimension (6,n,duree/env.pas)
  //                        Represents for each point (n) position and velocity (6)
  //                        during the trajectory.
  //                        if a trajectory ends before the others, it is completed by NaN
  //
  // Author:
  // A. BLAZQUEZ (CNES DCT/SB/MO)
  
  // Declarations:
  
  
  // Code:
  
  sizeTd = size(Traj_dir);
  
  i = 1;
  j = 1;
  
  //Initialisation (convergent or divergent)
  if Traj_type==0
    X0 = Traj_dir(1:6,i) + mod*Traj_dir(8:13,i);
    duree = abs(duree);
    f=CL__3b_RHS
  else
    duree = -abs(duree);
    X0 = Traj_dir(1:6,i) + mod*Traj_dir(14:19,i);
    f=CL__3b_RHSReb
  end
  
  tinit = Traj_dir(7,i);
  [X,t] = manifolds_3bIntegrVar(X0,tinit,duree,f,MU);
  
  sx = size(X);
  Traj_trace = hypermat([6,sx(2),1]);
  Traj_trace(:,:,1) = X;
  
  i =i+1;
  j =j+1;
  
  
  // loop for each point of Traj_dir
  
  while i < sizeTd(2)
    if Traj_type==0
      X0 = Traj_dir(1:6,i) + mod*Traj_dir(8:13,i);
    else
      X0 = Traj_dir(1:6,i) + mod*Traj_dir(14:19,i);
    end
  
    tinit = Traj_dir(7,i);
    [X,t] = manifolds_3bIntegrVar(X0,tinit,duree,f,MU);
  
  // Plot for testing
  // plot(X(1,:),X(2,:),'k');
  
    strace = size(Traj_trace);
    sx = size(X);
    Traj_tracestock = Traj_trace;
    if (strace(2) < sx(2))
      Traj_trace = hypermat([6,sx(2),j]);
      Traj_trace(:,1:strace(2),1:j-1) = Traj_tracestock;
      Traj_trace(:,strace(2)+1:sx(2),1:j-1) = ones((j-1)*6*(sx(2)-strace(2)))*%nan;
      Traj_trace(:,:,j) = X;
    else
      Traj_trace = hypermat([6,strace(2),j]);
      Traj_trace(:,:,1:j-1) = Traj_tracestock;
      Traj_trace(:,:,j) = [X %nan*ones(6,strace(2)-sx(2))];
    end
  
    i = i+1;
    j = j+1;
  
  end
  
endfunction


// NB: KEPT HERE FOR ASCENDING COMPATIBILITY PURPOSES

function [Traj_dir] = manifolds_dirconv(orb,T,env)
//Author
// 12/08/07 R.Alacevich et B. Meyssignac ( CNES, DCT/SB/MO)
// 12/02/09 A. Blazquez( CNES, DCT/SB/MO)

//fonction qui donne, a chaque point de la trajectoire regeneree, le vecteur propre des varietes convergentes et divergentes a partir de la matrice de monodromie calculee au bout du temps T.
//  entree: orb: 7xn: sortie de la regeneration,
//                             chaque colonne est un point de l'orbite: position vitesse et date
//              T: 1 : temps au bout duquel est estime la matrice de monodromie
//  sortie: matrice Traj_dir:
//    premieres 7 lignes: les memes que orb (point et date)
//    ligne 8-13: direction divergeante normee du point correspondant
//    ligne 14-19: direction convergeante normee du point correspondant


// Declarations:


// Code:

sizeTL = size(orb,2);
Traj_dir=[]
div_dir = [];
conv_dir = [];
CL__message("Directions computation for " + msprintf("%d",sizeTL) + " points\n");

for i=1:sizeTL
  CL__message("%d ", i);
  //calcul de la monodromie a T pour la direction divergente
  //et a -T pour la direction convergente
  [XTpos,monodTpos] = CL__3b_monodromy(orb(1:6,i),T,env.MU);
  [XTneg,monodTneg] = CL__3b_monodromy(orb(1:6,i),-T,env.MU);

  //diagonalisation des matrices de monodromie
  // AB A modifier apartir d'ici il a besoin que de vecteur propre de Lambda
  VectPropPos = CL__3b_monodroVectProp(monodTpos);
  VectPropNeg = CL__3b_monodroVectProp(monodTneg);

  div=real(VectPropPos)
  signe_div = sign(div(1));
  div_direction = signe_div*div/norm(div);

  convv=real(VectPropNeg)
  signe_conv = sign(convv(1));
  conv_direction = signe_conv*convv/norm(convv);

  div_dir = [div_dir,div_direction];
  conv_dir = [conv_dir,conv_direction];

end
CL__message("\n");

Traj_dir = [orb;div_dir;conv_dir];
endfunction


// NB: KEPT HERE FOR ASCENDING COMPATIBILITY PURPOSES

function Xdot = CL__3b_RHSReb(t,X,MU)
  //  Inputs: t: time
  //      X: vector of 6 lines (position+velocity)
  //          MU: from tb_environnement
  //
  //  Outputs: XCL_dot: vector of 6 lines (position+velocity)
  //
  // Author:
  // B. MEYSSIGNAG (CNES DCT/SB/MO)
  
  // Declarations:
  
  
  // Code:
  
  r0 = sqrt((X(1)+MU)^2+X(2)^2+X(3)^2);
  r1 = sqrt((X(1)-1+MU)^2+X(2)^2+X(3)^2);
  
  Xdot = [X(4);
          X(5);
          X(6);
         -2*X(5)+X(1)-(1-MU)*(X(1)+MU)/(r0^3)-MU*(X(1)-1+MU)/(r1^3);
          2*X(4)+X(2)-(1-MU)*X(2)/(r0^3)-MU*X(2)/(r1^3);
         -(1-MU)*X(3)/(r0^3)-MU*X(3)/(r1^3)];
  
endfunction

// initially : CL__3b_integratorVar
function [Y,t] = manifolds_3bIntegrVar(Y0,t0,tint,f,MU)


// Author:
// B. Meyssignac (CNES DCT/SB/MO)

//fonction qui utilise l'integration "root", pour le seul calcul de trajectoire
//  entrees: Y0: vecteur colonne des conditions initiales (six composantes: pos - vitesses)
//     t0: temps initial
//     tint: duree de l'integration
//               f: second membre: YCL_dot = f(t,Y). il n est pas gere ici en interne l'integration a rebours:
//Scilab dans sa version actuelle n'integre pas a rebours quand il fonctionne avec un pas variables.
//pour integrer a rebours il faut changer t en -t dans f
//Dans ce cas on integre non pas avec la fonction f mais avec f_reb
//qui doit etre definie dans la pile de variable. f_reb est  f dans laquelle t a ete change en -t.
//  sorties: Y: matrice resultat de l'integration: tous les points de la trajectoire
//     t: vecteur des dates associees a chaque point (colonne de Y)


// Declarations:


// Code:

%ODEOPTIONS = [2,0,0,%inf,0,2,500000,12,5,0,-1,-1];
sizeY = size(Y0);
rtol = 1.e-8*ones(sizeY(1),1);
atol = 1.e-9*ones(sizeY(1),1);
ng = 1;

if  tint>0

  deff("[u]=r(t,X)","u=t0+tint-t");
  [Yt,rd] = ode("root",Y0,t0,t0+tint,rtol,atol,f,ng,r);
  Y = Yt(2:$,:);
  t = Yt(1,:);

else

  deff("[u]=r(t,X)","u=t0+tint+t");
  Y0_reb = [Y0(1:3);-Y0(4:6)];
  [Yt,rd] = ode("root",Y0_reb,-t0,-(t0+tint),rtol,atol,f,ng,r);
  Y = [Yt(2:4,:);-Yt(5:7,:)];
  t = -Yt(1,:);

end

endfunction

// Declarations:


// Code:

CL__warnDeprecated(); // deprecated function

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

orb = [orb ; t_orb];  // compatibilite avec l'ancienne interface

CL__message("-> Convergent and divergent directions computation\n");

[orb_dir] = manifolds_dirconv(orb,period/2,env);
CL__message("-> Manifolds generation\n");


manifolds = list();
Np = size(pars,'*');

for i=1:Np

  par = pars(i);
  select par
    case "div"  //calcul de la variete divergente avec epsilon
      manifold = manifolds_trace(orb_dir,epsilon,tint,env.MU,0);
    case "-div" //calcul de la variete divergente avec -epsilon
      manifold = manifolds_trace(orb_dir,-epsilon,tint,env.MU,0);
    case "conv" //calcul de la variete convergente avec epsilon
      manifold = manifolds_trace(orb_dir,epsilon,tint,env.MU,1);
    case "-conv" //calcul de la variete convergente avec -epsilon
      manifold = manifolds_trace(orb_dir,-epsilon,tint,env.MU,1);
    else
      CL__error('Unknown name of manifolds (pars)');
  end

  manifolds($+1) = manifold;

end

lhs = argn(1);

if lhs==1 & Np <> 1
  varargout = list(manifolds);
elseif lhs==Np
  varargout = manifolds;
else
  CL__error('Bad size on input/output arguments');
end

endfunction