//  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 [Traj_trace] = CL__3b_traceManifold(Traj_dir, mod, dt, MU, Traj_type)
// 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 whether it's convergent or divergent
//
// Inputs: 
// Traj_dir: Matrix(19xn), n=number of points of the orbit
//   1:7: position, velocity and date
//   8:13: divergent direction on position and velocity 
//   14:19: Convergent direction on position and velocity
//   Traj_dir is usually generated by CL__3b_dirDirConv
// mod: perturbation module,(POSITIF ou NEGATIF)(real)
// dt: Extrapolation times relative to orbit times (>= 0)
// Traj_type: 0 Convergent manifold, 1 Divergent manifold (???)
//
// Outputs: 
// Traj_trace:  HyperMatrix. dimension (6 x nbpts x norb)
//   nbpts: number of points of each trajectory
//   norb: number of orbit points (input)
//
// Author:
// A. BLAZQUEZ (CNES DCT/SB/MO) - initial version 
// AL - new version (simpler)

// Declarations:


// Code:

// norb = number of input points 
norb = size(Traj_dir, 2); 

// number of points of each trajectory
nbpts = size(dt, 2); 
Traj_trace = hypermat([6,nbpts,norb]);


for k = 1 : norb
  t0 = Traj_dir(7,k); 

  if (Traj_type == 0)
    X0 = Traj_dir(1:6,k) + mod*Traj_dir(8:13,k);
    t = t0 + abs(dt);
  else
    X0 = Traj_dir(1:6,k) + mod*Traj_dir(14:19,k);
    t = t0 - abs(dt);
  end

  X = CL__3b_integratorFix(X0, t0, t, CL__3b_RHS, MU);
  Traj_trace(:,1:nbpts,k) = X; 

end

endfunction