function [E] = NewtonRaphson(e, M)
%NewtonRaphson Iteratively computes E, given M and e.
%   Iteratively computes elliptical anomaly E.
%   Requires eccentricity e, and initial guess M = 2*pi*(t-tp)/T.

f       = @(E) E-e*sin(E)-M;
fprim   = @(E) 1-e*cos(E);

% Newton-Raphson Method 
E = M;
Eprev = 0;
iterations = 0;
tol = 1e-6;
while (abs(E - Eprev) > tol)
    Eprev = E;
    E = E - f(E)/fprim(E);
    
    iterations = iterations + 1;
end
fprintf("E = %f rad / %f deg\n(Converged to within %1.0e in %d iterations.)\n", E, rad2deg(E), tol, iterations); 

end

