//  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 [X_t] = CL__3b_shooting(X0,t,MU,epsilon)
// This function uses a multiple shooting method to determine the orbit from reduce initial conditions
//
// Inputs: 
// X0: Halo Initial conditions - position+velocity (6x1)
// t: discretisation times of output orbit.
// MU: mass ratio (see env structure)
// epsilon: precision of the results.
// Output: 
// X_t: stabilized Trajectory 
//

// Author:
// A. BLAZQUEZ (CNES DCT/SB/MO)


// Declarations:


// Code:

max_iter = 50;
ecart_rel_max = 0.5

// initialisation
[Xa,norm_F] = CL__3b_shootingEach(X0,t,MU);

ecart_abs = max(norm_F);

if (ecart_abs < epsilon)
  // CL__message("Stabilization reached - no iterations \n");
  X_t = [X0; t]; 
  return;  // <== RETURN
end


ecart_abs_prec = ecart_abs;
ecart_rel = 0.1;
ecart_min = 15;
iter = 1;

while (ecart_abs > epsilon & ecart_rel <= ecart_rel_max & iter <= max_iter)
  // CL__message("Iteration %d: Error = %.2e\n",iter, ecart_abs); 
  if (ecart_abs < ecart_min)
    ecart_min = ecart_abs; 
    iter_min = iter; 
    X_min = Xa; 
  end

  [Xb, norm_F] = CL__3b_shootingEach(Xa, t, MU);
    
  // absolute and relative errors
  ecart_abs = max(norm_F); 
  ecart_rel = (ecart_abs - ecart_abs_prec) / ecart_abs_prec; 
  Xa = Xb;
  ecart_abs_prec = ecart_abs; 
  iter = iter + 1;

end

Xfin = Xa;

// CL__message("Iteration %d, Error: %.2e\n",iter, ecart_abs); 

if (ecart_rel > ecart_rel_max)
  // CL__message("ecart_rel = %5.3f > %5.3f\n", ecart_rel, ecart_rel_max); 
  // CL__message("Best Iteration: %d Error: %.2e\n",iter_min, ecart_min); 
  CL__warning("Approximate Stabilization"); 
  Xfin = X_min; 
end

if (iter > max_iter)
  CL__error("Stabilization not reached\n"); 
end

X_t = [Xfin; t]; 

endfunction