//  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'.

// -----------------------------------------------------------
//> Halo manifold around a given libration point.
//> Coordinates are adimensional 
// -----------------------------------------------------------

// Defaut values
ibodies = 1;
Lpoint = 1;
Az = 0.001;
idir = 1; // prograde
t_max = 3; // ~ 1 period
t_max_m = 2;

desc_param = list(..
   CL_defParam("Bodies (1=Sun-EarthMoon, 2=Earth-Moon)", ibodies, accv=1:2 ),...
   CL_defParam("Lagrange point (1=L1, 2=L2, 3=L3) ", Lpoint, accv=1:3 ),...
   CL_defParam("Amplitude along Z axis (adimensional)", Az, valid='$x > 0'),...
   CL_defParam("Motion direction / +X axis (1=prograde, 2=retrograde)",idir, accv=1:2),...
   CL_defParam("Propagation time for orbit (adimensional)", t_max, valid='$x > 0'),...
   CL_defParam("Propagation time for manifold (adimensional)", t_max_m, valid='$x > 0')...
);

[ibodies,Lpoint,Az,idir,t_max,t_max_m] = CL_inputParam(desc_param);

tab_env = ['S-EM' , 'E-M'];
tab_Lpoint = [ 'L1', 'L2', 'L3'];
env = CL_3b_environment(tab_env(ibodies),tab_Lpoint(Lpoint));

nb_dates = round(60 * t_max/(2*%pi));
t_orb = linspace(0,t_max,nb_dates);

nb_dates_m = max(10, round(1000 * t_max_m/(2*%pi)));
dt_int = linspace(0,t_max_m,nb_dates_m);

tab_direction = ["pro", "retro"]; 
[orb,omega_halo] = CL_3b_halo(env,Az,tab_direction(idir),t_orb);

epsilon = 3.e-4;
[div_in,conv_in,conv_out,div_out] = CL_3b_manifold(env,orb,t_orb,omega_halo,epsilon,dt_int,['div','conv','-conv','-div']);


// -----------------------------------------------------------
// plot
// -----------------------------------------------------------

f=scf();
f.immediate_drawing="off";
f.figure_size=[750,550];
nbcol = min(max(nb_dates,3),100); 
f.color_map = rainbowcolormap(nbcol) * 0.6;
f.visible="on";

Noir = addcolor([0,0,0]);
GrisF = addcolor([1,1,1]*0.4);
GrisC = addcolor([1,1,1]*0.6);

tab_system_names = ['(Sun-EarthMoon)' , '(Earth-Moon)'];
F = 1000; // unit factor

subplot(2,2,1);
a=gca();
a.title.text = "Halo manifold (conv in) - "+tab_Lpoint(Lpoint)+" - "+tab_env(ibodies) + " - unit = 0.001";
for i=1:nb_dates
  param3d((conv_in(1,:,i)-env.gl) * F, conv_in(2,:,i) * F, conv_in(3,:,i) * F);
  h = CL_g_select(gce(), "Polyline");
  h.foreground = max(1,min(nbcol,round(i * nbcol / nb_dates)));   
end

subplot(2,2,2);
a=gca();
a.title.text = "Halo manifold (conv out) - "+tab_Lpoint(Lpoint)+" - "+tab_env(ibodies) + " - unit = 0.001";
for i=1:nb_dates
  param3d((conv_out(1,:,i)-env.gl) * F, conv_out(2,:,i) * F, conv_out(3,:,i) * F);
  h = CL_g_select(gce(), "Polyline");
  h.foreground = max(1,min(nbcol,round(i * nbcol / nb_dates)));   
end

subplot(2,2,3);
a=gca();
a.title.text = "Halo manifold (div in) - "+tab_Lpoint(Lpoint)+" - "+tab_env(ibodies) + " - unit = 0.001";
for i=1:nb_dates
  param3d((div_in(1,:,i)-env.gl) * F, div_in(2,:,i) * F, div_in(3,:,i) * F);
  h = CL_g_select(gce(), "Polyline");
  h.foreground = max(1,min(nbcol,round(i * nbcol / nb_dates)));   
end

subplot(2,2,4);
a=gca();
a.title.text = "Halo manifold (div out) - "+tab_Lpoint(Lpoint)+" - "+tab_env(ibodies) + " - unit = 0.001";
for i=1:nb_dates
  param3d((div_out(1,:,i)-env.gl) * F, div_out(2,:,i) * F, div_out(3,:,i) * F);
  h = CL_g_select(gce(), "Polyline");
  h.foreground = max(1,min(nbcol,round(i * nbcol / nb_dates)));   
end

// adjustments
for j = 1 : 4
  subplot(2,2,j);
  a=gca();
  a.margins = [0.13, 0.13, 0.16, 0.1];  // left, right, top, bottom
  CL_g_stdaxes(a, colg=GrisC, ft=2, fl=1, fg=1);
  h=CL_g_select(a, "Polyline");
  h.thickness=2;
  a.rotation_angles = [75, 245]; 
end

f.immediate_drawing="on";
f.visible="on";