%orbit_angle_vs_time
% Plots orbit angle vs time for an elliptical orbit and compares to a
% cicular equivalent with the same period.
%
% LAST UPDATED by Andy French August 2022

function orbit_angle_vs_time

%Object name
object_name = 'pluto';

%Number of orbits
N = 3;

%Period /years
P = 248.348;

%Orbital eccentricity 
ecc = 0.25;

%Initial angle /rad
theta0 = 0;

%Fontsize for graph
fsize = 18;

%

%Compute orbit times and angles for circular case
t = linspace( 0,N*P,1000 );
theta_c = 2*pi*t/P;

%Numeric method to compute polar angle vs orbit time
theta = angle_vs_time( t, P, ecc, theta0 );

%Plot theta vs t and compare to circular version
figure('name','orbit angle vs time','color',[1 1 1]);
plot( t, theta_c, t, theta,'linewidth',2 ); grid on;
set(gca,'fontsize',fsize,'gridlinestyle','-');
ylabel('orbit polar angle /rad'); xlabel('time /years');
title(['Orbit angle vs time for ',object_name])
legend({'Circular \epsilon=0',['\epsilon=',num2str(ecc)]},'location','northwest')
print( gcf, 'orbit angle vs time.png','-dpng','-r300' );
close(gcf);

%%

%Numeric method to compute polar angle vs orbit time
function theta = angle_vs_time( t, P, ecc, theta0 )

%Angle step for Simpson's rule
dtheta = 1/1000;

%

%Number of orbits
N = ceil( t(end)/P );

%Define array of polar angles for orbits
theta = theta0 : dtheta : ( 2*pi*N + theta0 );

%Evaluate integrand of time integral
f = (1 - ecc*cos(theta) ).^(-2);

%Define Simpson rule coefficients c = [ 1, 4, 2, 4, 2, 4, ....1 ]
L = length(theta);
isodd = rem( 1:(L-2),2 ); isodd( isodd==1 ) = 4; isodd( isodd==0 ) = 2;
c = [1, isodd, 1];

%Calculate array of times
tt = P*( (1-ecc^2)^(3/2) )*(1/(2*pi))*dtheta*(1/3).*cumsum( c.*f );

%Interpolate the polar angles for the eccentric orbit at the circular orbit
%times
theta = interp1( tt, theta, t, 'spline' );

%End of code