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

// -----------------------------------------------------------
//> Lissajous manifold around a given libration point :
//> Coordinates are adimensional 
// -----------------------------------------------------------

// Defaut values
ibodies = 1;
Lpoint = 1;
Ax = 0.0005;
phix = 0;
Az = 0.001;
phiz = 0;
t_max = 3.14;
t_max_m = 1.5; 

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 X axis", Ax, valid='$x > 0'),...
   CL_defParam("Phase in X", phix, units=['rad', 'deg'], valid='$x >= 0'),...
   CL_defParam("Amplitude along Z axis", Az, valid='$x > 0'),...
   CL_defParam("Phase in Z", phiz, units=['rad', 'deg'], valid='$x >= 0'),...
   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,Ax,phix,Az,phiz,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);


epsilon = 1.e-10;
[lissajous_orb,omega,nu] = CL_3b_lissajous(env,Ax,phix,Az,phiz,epsilon,t_orb);

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

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

f=scf();
f.figure_size=[750,550];
f.visible="on";
f.immediate_drawing="off";
nbcol = min(max(nb_dates,3),100); 
f.color_map = rainbowcolormap(nbcol) * 0.6;
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 = "Lissajous 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 = "Lissajous 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 = "Lissajous 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 = "Lissajous 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";