//  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 [env] = CL_3b_environment(bodies, Lpoint, mrat, dist, omega)
// Structure creation for the 3-body problem functions
//
// Calling Sequence
// env = CL_3b_environment(bodies, Lpoint [, mrat, dist, omega])
//
// Description
// <itemizedlist><listitem>
// <p>Creates a structure to be used by the other "3b" functions.</p> 
// <p>This structure contains the necessary information about a specific libration point 
// in a specific system.</p>
// <p></p>
// <p>The body system can be: </p>
// <p>- "S-EM": Sun + barycenter of Earth-Moon </p>
// <p>- "E-M": Earth + Moon </p>
// <p>- "S-J": Sun + Jupiter </p>
// <p>- "-": user-defined  </p>
// <p></p>
// <p>If the body system is "-", the 3 additional arguments (<b>mrat</b>: mass ratio, <b>dist</b>: 
// distance between the 2 main bodies, <b>omega</b>: angular velocity of the line between the 2 main bodies)
// must be given, otherwise they should be either absent or empty ([]). </p>
// <p></p>
// <p>The Lagrange point can be either "L1", "L2" or "L3" ("L4" and "L5" are not accepted). </p>
// <p><imageobject><imagedata fileref="3b_frame.png"/></imageobject></p>
// <p></p>
// </listitem>
// <listitem>
// <p> Description of the output structure:</p>
// <p>- bodies: (string) System of primaries ("S-EM", "E-M", "S-J" or "-")</p>
// <p>- Lpoint: (string) Libration point ("L1", "L2" or "L3")</p>
// <p>- D: Distance between the primaries (m) </p>
// <p>- OMEGA: Inertial rotation rate of the primaries around their center of mass (rad/s) </p>
// <p>- MU: Mass ratio (in ]0, 1]): mass of smallest primary / mass of biggest primary (=mu2/mu1)</p>
// <p>- gammal: Adimensional distance (positive) between the libration point and the closest primary: d(P2,L1), d(P2,L2) or d(P1,L3) </p>
// <p>- gl: Adimensional abscissa (positive or negative) from the center of mass (G) to the libration point: G-L1 or G-L2 or G-L3</p>
// <p></p>
// <p> Other quantities for internal use:</p>
// <p>- omega_init: Theoretical pulsation of in-plane motion for orbits around the libration point</p>
// <p>- pas: internal generation step (adimensional)</p>
// <p>- k </p>
// <p>- c2 </p>
// <p>- nu </p>
// <p>- invM </p>
// <p>- ESCAPEDIR </p>
// <p>- ESCAPEDIRNORM </p>
// <p><imageobject><imagedata fileref="3b_L2_param.png"/></imageobject></p>
// <p></p></listitem>
// <listitem>
// <p><b>Note</b>:</p>
// <p>Conversion between "dimensional" and "adimensional" quantities:</p>
// <p>- for lengths: dx = dX / D </p>
// <p>- for durations: dt = dT * OMEGA </p>
// <p>Where: </p>
// <p>- dx and dX: respectively adimensional and dimensional lengths, </p>
// <p>- dt and dT: respectively adimensional and dimensional time durations, </p>
// <p>- D and OMEGA: as described above. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// bodies: (string) Body system: "S-EM", "E-M", "S-J" or "-"
// Lpoint: (string) Type of libration point: "L1", "L2" or "L3"
// mrat: (optional) Mass ratio = mass of smaller body divided by mass of bigger one if bodies == "-". (1x1)
// dist: (optional) Distance between the 2 main bodies [m] if bodies == "-". (1x1) 
// omega: (optional) Angular velocity of the line between the 2 main bodies in inertial frame [rad/s] if bodies == "-". (1x1) 
// env: Resulting structure. 
//
// Authors
// CNES - DCT/SB
//
// Examples
// env = CL_3b_environment("S-EM","L1")
//

// Error checking
if (~exists("mrat", "local")); mrat = []; end
if (~exists("dist", "local")); dist = []; end
if (~exists("omega", "local")); omega = []; end

if (size(bodies, "*") <> 1)
  CL__error("Invalid size for bodies"); 
end

if (size(Lpoint, "*") <> 1)
  CL__error("Invalid size for Lpoint"); 
end

if (find(Lpoint == ["L1", "L2", "L3"]) == [])
  CL__error("Invalid value for Lpoint");
end

if (find(bodies == ["S-EM", "E-M", "S-J", "-"]) == [])
  CL__error("Invalid value for bodies");
end

if ( (bodies == "-" & (mrat == [] | dist == [] | omega == [])) | ..
     (bodies <> "-" & (mrat <> [] | dist <> [] | omega <> [])) )
  CL__error("Invalid value for mrat, dist or omega");
end

if (bodies == "-" & (mrat <= 0 | mrat > 1 | dist <= 0 | omega <= 0)) 
  CL__error("Invalid value for mrat, dist or omega");
end

// Load constants for three-body problem (if bodies <> "-")
I = find(bodies == ["S-EM", "E-M", "S-J"]);
if (I<>[])
  tabnames = ["Sun_EarthMoon", "Earth_Moon", "Sun_Jupiter"];
  name = tabnames(I);
  str = "threebody." + name + [".mrat", ".dist", ".omega"];
  mrat = CL_dataGet(str(1)); 
  dist = CL_dataGet(str(2)); 
  omega = CL_dataGet(str(3)); 
end

// old notations
MU = mrat; 
D = dist; 
OMEGA = omega; 

// distance to closest primary
gammal = CL__3b_gamma(MU,Lpoint);

select Lpoint
  case "L1"
    c2 = 1/(gammal^3)*(MU+(1-MU)*(gammal^3)/((1-gammal)^3));
    gl = 1 -gammal - MU;
    
  case "L2"
    c2 = 1/(gammal^3)*(MU+(1-MU)*(gammal^3)/((1+gammal)^3));
    gl = gammal + 1 - MU;
    
  else // case "L3"
    c2 = 1/(gammal^3)*(1-MU+MU*(gammal^3)/((1+gammal)^3));
    gl= -gammal -MU;
    
end

// internal generation step
pas = 0.01; 

lambda = sqrt(1/2*(-2+c2+sqrt(9*c2^2-8*c2)));
omega_init = sqrt(1/2*(2-c2+sqrt(9*c2^2-8*c2)));
nu = sqrt(c2);
c = (lambda^2-1-2*c2)/(2*lambda);
k = (omega_init^2+1+2*c2)/(2*omega_init);
b1 = c*lambda + k*omega_init;
b2 = c*omega_init - k*lambda;

// Linearization of motion
// coordinates:  x-gl, y, z, vx, vy, vz (rows)
// in basis of eigen vectors (columns):
//  A*e^(lambda*t), B*e^-(lambda*t), C*e^(i*omega_init*t), D*e^(-i*omega_init*t), E*e^(i*nu*t), F*e^(-i*nu*t)

M = [1          1          1               1              0        0;
     c         -c          %i*k           -%i*k           0        0;
     0          0          0               0              1        1;
     lambda    -lambda     %i*omega_init  -%i*omega_init  0        0;
     lambda*c   lambda*c  -omega_init*k   -omega_init*k   0        0;
     0          0          0               0              %i*nu   -%i*nu];

// matrix of eigen vectors (rows) = inverse of M: 
// A*e^(lambda*t), B*e^-(lambda*t), C*e^(i*omega_init*t), D*e^(-i*omega_init*t), E*e^(i*nu*t), F*e^(-i*nu*t)
// in basis:  x-gl, y, z, vx, vy, vz (columns)

invM = 1/2* [k*omega_init/b1   omega_init/b2     0  -k/b2      1/b1   0;
             k*omega_init/b1  -omega_init/b2     0   k/b2      1/b1   0;
             c*lambda/b1       %i*lambda/b2      0  -%i*c/b2  -1/b1   0;
             c*lambda/b1      -%i*lambda/b2      0   %i*c/b2  -1/b1   0;
             0                 0                 1   0         0     -%i/nu;
             0                 0                 1   0         0      %i/nu];

// escape_direction vector in plane (vx, vy)
ESCAPEDIR = real(invM(1,4:5).');
ESCAPEDIRNORM = ESCAPEDIR / norm(ESCAPEDIR);

// non escape_direction vector in plane (vx, vy)
NESCAPEDIR = [1/b1;k/b2];
NESCAPEDIRNORM = 1/norm(NESCAPEDIR)*NESCAPEDIR;

// definition of structure 
env = struct('bodies',bodies, 'Lpoint',Lpoint, ...
             'MU',MU, 'D',D, 'OMEGA',OMEGA, 'gammal',gammal,  ...
             'gl',gl, 'omega_init',omega_init, 'pas', pas, 'k',k,   ...
             'c2',c2, 'nu',nu, 'invM',invM, 'ESCAPEDIR',ESCAPEDIR, ...
             'ESCAPEDIRNORM',ESCAPEDIRNORM);
             
endfunction