//  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 orbit 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 = 6.28;

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

[ibodies,Lpoint,Ax,phix,Az,phiz,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);


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


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

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

f.color_map = jetcolormap(10) * 0.8;

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=1);
a=gca();
a.title.text = "Lissajous 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=1);
a=gca();
a.title.text = "Lissajous orbits 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=1);
a=gca();
a.title.text = "Lissajous orbits 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 = 1;
a=gca();
a.title.text = "Lissajous orbits 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";