//  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 [X0,omegacorr] = CL__3b_condInitHalo(Az, direction, env)
// Author:
// B. Meyssignac ( CNES, DCT/SB/MO)
//
// Genere les conditions initiales de halo en dynamique approchee autour du point de
// Lagrange choisi a partir de l amplitudes et du sens donnes en entree
// On utilise ici la methode de Richardson (D.L. Richardson 1979) pour trouver
// des conditions initiales proche d une halo dans le plan XZ
// Ca donne des conditions initiales du type (x0;0;z0;0;vy0;0) 
// (les halo sont symetriques par rapport au plan XZ)
// 
// entree : 
// Az: amplitude de l'orbite de halo desires: .
//    Ax > Axmin car sinon  les effets non lineaires sont trop faibles
//    pour donner deux frequences propres (plan et hors plan) egales
// direction: direction of motion around the +X axis
//    +1: direct (prograde) around +X axis
//    -1: indirect (retrograde) around +X axis
//
// sortie : 
// vecteur colonne des conditions initiales correspondantes dans le repere synodique normalise
// et om pulsation de la halo corrigee
//
// constantes pour la methode de Richardson
// et normalisation en distance par rapport au point de lagrange considere

// Declarations:


// Code:
if (direction <> 1 & direction <> -1)
  CL__error("Invalid value for direction"); 
end

Lpoint = env.Lpoint;
gammal = env.gammal;
MU = env.MU;
omega = env.omega_init;
k = env.k;
gl = env.gl;
c2 = env.c2

if (Lpoint == "L1")
  c3 = 1/(gammal^3)*(MU-(1-MU)*(gammal^4)/((1-gammal)^4));
  c4 = 1/(gammal^3)*(MU+(1-MU)*(gammal^5)/((1-gammal)^5));

elseif  (Lpoint == "L2")
  c3 = 1/(gammal^3)*(-MU-(1-MU)*(gammal^4)/((1+gammal)^4));
  c4 = 1/(gammal^3)*(MU+(1-MU)*(gammal^5)/((1+gammal)^5));

else // if (Lpoint == "L3")
  c3 = 1/(gammal^3)*(1-MU+MU*(gammal^4)/((1+gammal)^4));
  c4 = 1/(gammal^3)*(1-MU+MU*(gammal^5)/((1+gammal)^5));
end

Az = Az/gammal;

delta = omega^2-c2;

d1 = 3*omega^2*(k*(6*omega^2-1)-2*omega)/k;
d2 = 8*omega^2*(k*(11*omega^2-1)-2*omega)/k;
d21 = -c3/(2*omega^2);

b21 = -3*c3*omega/(2*d1)*(3*k*omega-4);
b22 = 3*c3*omega/d1;

a21 = 3*c3*(k^2-2)/(4*(1+2*c2));
a22 = 3*c3/(4*(1+2*c2));
a23 = -(3*c3*omega*(3*k^3*omega-6*k*(k-omega)+4))/(4*k*d1);
a24 = -(3*c3*omega*(2+3*k*omega))/(4*k*d1);
a31 = -9*omega/(4*d2)*(4*c3*(k*a23-b21)+k*c4*(4+k^2))+(9*omega^2+1-c2)/(2*d2)*(3*c3*(2*a23-k*b21)+c4*(2+3*k^2));
a32 = -(9*omega/4*(4*c3*(k*a24-b22)+k*c4)+3/2*(9*omega^2+1-c2)*(c3*(k*b22+d21-2*a24)-c4))/d2;

a1 = -3/2*c3*(2*a21+a23+5*d21)-3/8*c4*(12-k^2);
a2 = 3/2*c3*(a24-2*a22)+9/8*c4;

b31 = 3/(8*d2)*(8*omega*(3*c3*(k*b21-2*a23)-c4*(2+3*k^2))+(9*omega^2+1+2*c2)*(4*c3*(k*a23-b21)+k*c4*(4+k^2)));
b32 = 1/d2*(9*omega*(c3*(k*b22+d21-2*a24)-c4)+3/8*(9*omega^2+1+2*c2)*(4*c3*(k*a24-b22)+k*c4));

d31 = 3/(64*omega^2)*(4*c3*a24+c4);
d32 = 3/(64*omega^2)*(4*c3*(a23-d21)+c4*(4+k^2));

s1 = 1/(2*omega*(omega*(1+k^2)-2*k))*(3/2*c3*(2*a21*(k^2-2)-a23*(k^2+2)-2*k*b21)-3/8*c4*(3*k^4-8*k^2+8));
s2 = 1/(2*omega*(omega*(1+k^2)-2*k))*(3/2*c3*(2*a22*(k^2-2)+a24*(k^2+2)+2*k*b22+5*d21)+3/8*c4*(12-k^2));

l1 = a1+2*omega^2*s1;
l2 = a2+2*omega^2*s2;


// en dessous d'une amplitude minimum les effets des non linearites sont trop faibles
// pour produire des frequences propres egales entre le mouvement plan et le mouvement hors plan
Ax = sqrt((-l2*Az^2-delta)/l1);

if (Ax < (sqrt(abs(delta/l1))))
  CL__error("Too small halo amplitude");
end

// pulsation corrigee
om = 1+s1*Ax^2+s2*Az^2;
omegacorr = om*omega;

x0 = a21*Ax^2+a22*Az^2-Ax+a23*Ax^2-a24*Az^2+a31*Ax^3-a32*Ax*Az^2;
y0 = 0;
z0 = Az-2*d21*Ax*Az+d32*Az*Ax^2-d31*Az^3;

vx0 = 0;
vy0 = omegacorr*k*Ax+2*omegacorr*(b21*Ax^2-b22*Az^2)+3*omegacorr*(b31*Ax^3-b32*Ax*Az^2);
vz0 = 0;

// Adjust the sign of z0 so that the motion is prograde or retrograde 
// (as desired) in the (Y,Z) plane around +X axis. 
// Prograde motion (direction = +1) <=> z0 * vy0 <= 0

if (direction * z0 * vy0 > 0); 
  z0 = -z0;
end

X0 = [x0; y0; z0; vx0; vy0; vz0];

// renormalisation par rapport a AU et
// par rapport au repere centre sur le centre de gravite du systeme a trois corps
if (Lpoint == "L1")
  X0 = (1-MU-gl)*X0 + [gl;0;0;0;0;0];
elseif (Lpoint == "L2")
  X0 = (gl-1+MU)*X0 + [gl;0;0;0;0;0];
else // if (Lpoint == "L3")
  X0 = (-gl-MU)*X0 + [gl;0;0;0;0;0];
end


endfunction