# 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%.

Download KeplerEquationError.m

% 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'); ylabel('radial error'); legend(legendText,'Location','Northwest') % Save plot as pdf set(gcf, 'PaperPosition', [0 2 8 8]); print ('-dpdf','kepler.pdf'); print ('-dpng','kepler.png','-r300');