function [ P, t, r, theta, x, y ] = orbit( Mstar_Solmass, aAU, e )

% orbit : [ P, t, r, theta, x, y ] = orbit( Mstar_Solmass, aAU, e )
% Given the mass of a massive central object Mstar_Solmass (given
% in solar masses) around which a smaller mass is orbiting at a semi-major
% axis aAU (in AU) with eccentricity e, this program calculates the period
% P of the orbit (in seconds) and also generates a  timestep t (in seconds), 
% the r and theta coordinates corresponding to t, and the x and y 
% coordinates corresponding to t.  Note that r, x, and y are in SI units 
% (meters), and theta is in radians.  The origin of coordinates is the 
% focus of the ellipse.

% Include standard constants
constants

% Convert conventional units to SI units
Mstar = Mstar_Solmass*M_Sun;
a = aAU * AU;

% Period of the orbit, in seconds
P = sqrt(4*pi^2*a^3/(G*Mstar));

% Number of time steps.  This is a variable we can control.  The larger the
% number of time steps, the smaller the time interval, and the smaller the
% error in calculating numerical derivatives.
n = 100000;
% Divide the period of the orbit into equal time steps.
dt = P/(n-1);
% Make a list of time steps from 0 to the period with equal spacing given
% by dt.
t = 0 : dt : P;
% The theta step is not uniform; we will calculate it below.  Here we just
% define a list to hold the values of theta we'll calculate.
theta = zeros(1,n+1);
% Similarly for r.
r = zeros(1,n);

% Calculate the angular momentum per unit mass, L/m (Eq. 2.30).  Since this
% quantity is conserved, we only need to calculate it once.
LoM = sqrt(G*Mstar*a*(1 - e^2));

% Now we calculate for every time in the list.
for i = 1 : n
    % First, calculate the radial position corresponding to the current
    % theta.  We start (by definition) with theta = 0 at perihelion.
    r(i) = a*(1 - e^2)/(1 + e*cos(theta(i)));
    % Compute the next value for theta using the fixed time step by 
    % combining Eq. (2.31) with Eq. (2.32), which is Kepler's Second Law.
    dtheta = LoM/r(i)^2*dt;
    theta(i+1) = theta(i)+dtheta;
end

% We needed one extra theta point to get the last value of r, but here
% we'll trim it off.
theta = theta(1:n);

% Convert the r, theta coordinates to x, y so we have these as well.
[x,y] = pol2cart(theta,r);

end