//  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 gammal = CL__3b_gamma(MU,Lpoint)
// Author
// 20/09/2010 E.Canalias (CNES, DCT/SB/MO) - initial version
// 16/11/2013 AL - improvements 
//
// Distance between libration point and closest primary
// Based on GR_GAM.F from Development of a software library for libration point mission analysis
// CNES Contract R-S06/VF-0001-052
// MU: Mass parameter (in ]0,1]) = m2/m1 (smallest mass / biggest mass)
// Lpoint: (string) Libration point ("L1", "L2" or "L3")
// gammal: Distance (>0) from the libration point to the closest primary in RTBP adimensional units

// Declarations:

// Code:

// internal function
function [u] = tb_gamma_H(Lpoint, X, MU)

  select Lpoint
    case 'L1'
      mu2 = 2*MU;
      a4 = 3-MU;
      a3 = 3-mu2; 
      u1 = -MU+X*(mu2+X*(-MU+X*(a3+X*(-a4+X))));
      u2 = mu2+X*(-mu2+X*(3*a3+X*(-4*a4+5*X))); 
      u = u1/u2; 
    
    case 'L2'
      mu2 = 2*MU; 
      a4 = 3-MU; 
      a3 = 3-mu2; 
      u1 = -MU+X*(-mu2+X*(-MU+X*(a3+X*(a4+X)))); 
      u2 = -mu2+X*(-mu2+X*(3*a3+X*(4*a4+5*X))); 
      u = u1/u2; 
    
    case 'L3'
      mum = 1-MU; 
      a3 = 1+2*MU; 
      a4 = 2+MU; 
      u1 = -mum+X*(-2*mum+X*(-mum+X*(a3+X*(a4+X)))); 
      u2 = -2*mum+X*(-2*mum+X*(3*a3+X*(4*a4+5*X))); 
      u = u1/u2;  
      
 end
 
endfunction

// initialisation of "X" argument for tb_gamma_H
function [X] = tb_gamma_initX(Lpoint, MU)

  select Lpoint
    case 'L1'
      X = (MU/3)^(1/3);
    case 'L2'
      X = (MU/3)^(1/3);
    case 'L3'
      X = 1-7*MU/12;
  end
 
endfunction


maxiter = 20;
tol = %eps; // tolerance 
 
X = tb_gamma_initX(Lpoint, MU); 
dX = %inf;
iter = 1;

while (abs(dX) > tol & iter <= maxiter)
  dX = tb_gamma_H(Lpoint, X, MU);      
  X = X - dX;
  iter = iter+1;
end

if (iter > maxiter); 
  CL__message("No convergence in computation");
end 
 
gammal = X;
 
endfunction