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