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

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 (adimensional)", t_max, valid='$x > 0')...
);

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


// -----------------------------------------------------------
// Computation
// -----------------------------------------------------------

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

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

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


// -----------------------------------------------------------
// computation / plot
// -----------------------------------------------------------

f=scf();
f.immediate_drawing="off";
f.figure_size=[750,550];
f.visible="on";

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

tab_system_names = ['(Sun-EarthMoon)' , '(Earth-Moon)'];

// unit factor
F = 1000; 

subplot(2,2,1);
plot2d((orb(1,:)-env.gl)*F, orb(2,:)*F, style=2);
a=gca();
a.title.text = "Halo orbit at "+tab_Lpoint(Lpoint)+" "+tab_system_names(ibodies) + " - unit = 0.001";
a.x_label.text = "X";
a.y_label.text = "Y";

subplot(2,2,2);
plot2d((orb(1,:)-env.gl)*F, orb(3,:)*F, style=2);
a=gca();
a.title.text = "Halo orbit at "+tab_Lpoint(Lpoint)+" "+tab_system_names(ibodies) + " - unit = 0.001";
a.x_label.text = "X";
a.y_label.text = "Z";

subplot(2,2,3);
plot2d(orb(2,:)*F, orb(3,:)*F, style=2);
a=gca();
a.title.text = "Halo orbit at "+tab_Lpoint(Lpoint)+" "+tab_system_names(ibodies) + " - unit = 0.001";
a.x_label.text = "Y";
a.y_label.text = "Z";

subplot(2,2,4);
param3d((orb(1,:)-env.gl)*F, orb(2,:)*F, orb(3,:)*F);
h = CL_g_select(gce(), "Polyline");
h.foreground = 2;
a=gca();
a.title.text = "Halo orbit at "+tab_Lpoint(Lpoint)+" "+tab_system_names(ibodies) + " - unit = 0.001";


// adjustments
for j = 1 : 4
  subplot(2,2,j);
  a=gca();
  a.margins = [0.15, 0.15, 0.2, 0.2];  
  CL_g_stdaxes(a, colg=GrisF, ft=2, fl=1, fg=1);
  h=CL_g_select(a, "Polyline");
  h.thickness=2;
end

subplot(2,2,4);
a=gca();
a.rotation_angles = [62, 273]; 

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