% KeplerEquation.m
%
% Follow an orbiting body through one period of an elliptical orbit
% Use Bessel's solution to Kepler's equation
%
% MATLAB-Monkey.com  9/28/2013


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

clear           % clear variables
clc             % clear command window

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

nTerms = 10;    % number of terms to keep in infinite series defining
                % eccentric anomaly

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


% 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
        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));
r = a * (1-e^2) ./ (1 + e*cos(theta));

% polar plot of orbit
polar(theta,r,'bo');

% Label plot with eccentricity
title(['Elliptical Orbit with e = ' sprintf('%.2f',e)]);

% Save plot as pdf
set(gcf, 'PaperPosition', [0 2 8 8]);
print ('-dpdf','kepler.pdf');
print ('-dpng','kepler.png','-r300');

