//  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 [M] = CL_kp_E2M(ecc,E)
// Eccentric anomaly to mean anomaly 
//
// Calling Sequence
// M = CL_kp_E2M(ecc,E)
//
// Description
// <itemizedlist><listitem>
// <p>Computes mean anomaly from eccentric anomaly for elliptical 
// and hyperbolic orbits:</p>
// <p> - Elliptical orbits ( ecc < 1 ) : <b>M = E - ecc * sin(E)</b> </p>
// <p> - Hyperbolic orbits ( ecc > 1 ) : <b>M = ecc * sinh(E) - E</b> </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// E: Eccentric anomaly [rad] (PxN or Px1)
// ecc: Eccentricity (PxN or Px1)
// M: Mean anomaly [rad] (PxN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Orbital Mechanics for Engineering Students, H D Curtis, Chapter 3, equation 3.11
// 2) Mecanique spatiale, CNES - Cepadues 1995, Tome I, section 4.6.6, page 181
//
// See also
// CL_kp_M2E
//
// Examples
// // Elliptic and hyperbolic orbit: 
// E = [1.7, 5.7];
// ecc = [0.3, 1.5];
// M = CL_kp_E2M(ecc,E);
//
// // Consistency test:
// E - CL_kp_M2E(ecc,M) // => 0

  
// Declarations:

// Code:

I = find(ecc < 0); 
if (I <> []) 
  CL__error("Invalid eccentricity (<0)");
end

// check / resize
nr = size(ecc,1); 
[ecc,E] = CL__checkInputs(ecc,nr,E,nr); 

M = %nan * ones(E);

// ------------------------
// Parabolic orbit: 
// result is %nan
// ------------------------

// ------------------------
// Ellipse
// ------------------------
// Standard case
Cond = (ecc < 0.999 | E > 0.1); 

I = find(Cond & ecc < 1); 
M(I) = E(I) - ecc(I).*sin(E(I));

// Special case: numerical accuracy would not be good with above formula.
// Use a specific function to compute E-e*sin(E)
// NB: These tresholds have been found so that a precision of 1e-13 is guaranteed.
I = find((~Cond) & ecc < 1); 
M(I) = CL__E_esinE(ecc(I),E(I));

// ------------------------
// Hyperbola
// ------------------------
I = find(ecc > 1);  
M(I) = ecc(I).*sinh(E(I)) - E(I);


endfunction