# Addendum to Kepler's Equation: Error Analysis

In the code KeplerEquation.m we tracked an orbiting body through one period of an elliptical orbit. It used a series expansion involving Bessel functions to solve Kepler's equation. In the code keplerEquationError.m we evaluate the accuracy of this method when a finite number of terms are included in the series expansion.

The following plot shows the error as a function of the angle theta for three different cases: when the infinite series is truncated at the 3rd, 5th and 7th terms for an ellipse with an eccentricity of 0.4. With only 5 terms, the error as dropped to less than 0.1%.

% KeplerEquationError.m
%
% This code estimates errors in the orbital positions calculated using
% the Bessel function expansion of Kepler's equation (see KeplerEquation.m)
% It plots the errors as a function of the angle theta for three different
% cases:  when the infinite series is truncated at the 3rd, 5th and 7th terms
%
% MATLAB-Monkey.com  10/1/2013

%%%%%%%%%%
%%%%%%%%%%  Parameters & Initialization
%%%%%%%%%%

clear           % clear variables
clc             % clear command window

e = 0.4;        % eccentricity
a = 1;          % semi-major axis
N = 100;        % number of points defining orbit

nTerms = [3 5 7]; % number of terms kept in eccentric anomaly series

M = linspace(0,2*pi,N);   % mean anomaly parameterizes time
% M varies from 0 to 2*pi over one orbit

alpha = zeros(1,N);       % preallocate space for eccentric anomaly array

%%%%%%%%%%
%%%%%%%%%%  Calculations & Plotting
%%%%%%%%%%

%%%%%%%%%%  Estimate the "true" solution by using 30 terms in the infinite series

% Calculate eccentric anomaly at each point in orbit
for j = 1:N
% initialize eccentric anomaly to mean anomaly
alpha(j) = M(j);

% include first nTerms in infinite series
for n = 1:30
alpha(j) = alpha(j) + 2 / n * besselj(n,n*e) .* sin(n*M(j));
end
end

% calcualte polar coordiantes (theta, r) from eccentric anomaly
theta30 = 2 * atan(sqrt((1+e)/(1-e)) * tan(alpha/2));
theta30(theta30<0) = theta30(theta30<0) + 2*pi;
r30 = a * (1-e^2) ./ (1 + e*cos(theta30));

%%%%%%%%%%  Orbit calculation with restricted number of terms

% loop through elements of matrix nTerms
for k = 1:length(nTerms)

% Calculate eccentric anomaly at each point in orbit
for j = 1:N
% initialize eccentric anomaly to mean anomaly
alpha(j) = M(j);

% include first nTerms in infinite series
for n = 1:nTerms(k)
alpha(j) = alpha(j) + 2 / n * besselj(n,n*e) .* sin(n*M(j));
end
end

% calcualte polar coordiantes (theta, r) from eccentric anomaly
theta = 2 * atan(sqrt((1+e)/(1-e)) * tan(alpha/2));
theta(theta<0) = theta(theta<0) + 2*pi;
r = a * (1-e^2) ./ (1 + e*cos(theta));

% find error
thetaError = theta - theta30;
rError = r - r30;

% plot errors as a function of theta
subplot(2,1,1)
plot(theta30,thetaError)
hold all

subplot(2,1,2)
plot(theta30,rError)
hold all

% create entry for legend indicating max number of terms
legendText{k} = sprintf('%i terms',nTerms(k));
end

% label plots and create legend
subplot(2,1,1)
title(sprintf('Truncation Error for Kepler''s equation (e = %.2f)',e));
xlim([0 2*pi]);
xlabel('\theta');
ylabel('\theta error');
legend(legendText,'Location','Northwest')

subplot(2,1,2)
xlim([0 2*pi]);
xlabel('\theta');