%% Eccentric Anomaly Numerical Calculation
% Author: Jack Naylor
% SID: 470354333

% Description:
% This function will numerically solve the Kepler equation:
%
%                   M = E-e*sin(E)
%
% using a 3rd order accurate numerical scheme centred around the angle E.

% Inputs:
% e = eccentricity
% M = mean anomaly (deg)
% tol = tolerance

% Outputs:
% E = eccentric anomaly (deg)

function E = eccentricAnomaly(e, M, tol)
%% Setup
% Limit the mean anomaly input to be between 0 and 2*pi
M = mod(deg2rad(M),2*pi);

% Use a 3rd order Taylor expansion expression 
E0 = M+(-0.5*e^3+e+(e^2+3/2*cos(M)*e^3)*cos(M))*sin(M);

% Set the defined difference in the true value of E(n) and E(n-1) to be
% larger than the tolerence to enter the iteration
dE = tol+1;

%% Iterative solver
% Run until the difference is lower than the specified tolerance
while dE>tol
    
    % Calculate the new E using the initial guess and a defined error
    % function
    E = E0 - epsE(e,M,E0);
    
    % Calculate the difference between this value and the last
    dE = abs(E-E0);
    
    % Set the current E value to the "initial" one in the scheme
    E0 = E;
    
end

%% Output
% Output eccentric anomaly in degrees
E = rad2deg(E);

end

%% Epsilon Function
% Calculate the epsilon value in a separate sub-routine for neatness

function out = epsE(e,M,x)
%% Setup
% Define components of expression
exp1 = cos(x);
exp2 = -1+e*exp1;
exp3 = sin(x);
exp4 = e*exp3;
exp5 = -x+exp4+M;
exp6 = exp5/(0.5*exp5*exp4/exp2+exp2);

%% Output
% Calculate expression and output
out = exp5/((0.5*exp3-1/6*exp1*exp6)*e*exp6+exp2);

end