//  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'.


// ------------------------------------------------------------------
// Solution of kepler's equation (M to E) - Hyperbola only
// (EQ)  M = E*sinh(E) - E   
// NB: 
// - not supposed to be called directly 
// - no tests on arguments (eccentricity) 
// - ecc and M must have the same dimension
// ------------------------------------------------------------------

function [E] = CL__kp_M2Ehyp(ecc,M)

  // Initial value for E: 
  // asinh(M/e) or (6M/e)^(1/3) (if M > 0)  
  // obtained by : 
  // - asinh(M/e)   : neglecting E in (EQ) 
  // - (6M/e)^(1/3) : developing exp(E) = 1+E+E^2/2+E^3/6  in (EQ) and neglicting first order
  // Note that we have : asinh(M/e) <=   F   <=   ( 6M/e )^(1/3)    

  E = asinh(M ./ ecc); // default value

  I = find( (ecc-1).^2 + (0.5*abs(M)-1).^2 < 1.01 );  // empirical - based on tests
  E(I) = sign(M(I)) .* (6 * abs(M(I)) ./ ecc(I)).^(1/3); 

  // -------------------------------------------
  // resolution iterative par methode de newton
  //    f(x) = 0
  //    xi+1 = xi - (f(xi) / f'(xi))
  // -------------------------------------------

  nb_max_iter = 20;         // max number of iterations
  errmax = 100 * %eps;      // accuracy 

  err1 = ones(E) * %inf;    // initial value for accuracy variables
  err2 = ones(E) * %inf; 
  K = 1:length(E);          // indices to be dealt with
  iter = 0;                 // iter = current number of iterations

  while (K <> [] & iter <= nb_max_iter) 

    E1 = E;  // previous value 
    E(K) = E(K) - (ecc(K) .* sinh(E(K)) - E(K) - M(K)) ./ (ecc(K) .* cosh(E(K)) - 1); 
 
    err1(K) = abs(E1(K)-E(K)) ./ max(abs(E(K)),1);    
    err2(K) = abs(ecc(K) .* sinh(E(K)) - E(K) - M(K)) ./ max(abs(M(K)),1); 

    K = find( err1 > errmax | err2 > errmax );
    iter = iter + 1; 

  end

  if (K <> []) 
    CL__error('Maximum number of iterations reached'); 
  end

endfunction