%Demonstration of how to calculate the orbit of a planet about a star, from
%observations of polar angles of the Sun and their time derivative. 
%Could log these using the angular deviations of a distant star in the night sky over a year?
%
% LAST UPDATED by Dr Andy French March 2025

function earth_orbit_calc

%Start from known orbital parameters!
T = 10.000;   %Orbital period (years)
dt = T/52;   %Time interval for measurements (years)
ecc = 0.5;   %Orbital eccentricity (Earth is 0.01, Mars is 0.09).
M = 3;   %Star mass (in solar masses)
theta0 = 0;  %Initial polar angle /rad.
a = (M*T^2)^(1/3);   %Semi-major axis of orbit /AU

%Calculate times (years) and polar angles (rad)
t = 0: dt : T;
theta = angle_vs_time( t, T, ecc, theta0 );

%Calculate rate of change of angle (rad/year)
theta_dot = (2*pi/T)* ( (1-ecc^2)^(-3/2) ) .*(1-ecc*cos(theta)).^2 ;

%Plot sqrt( theta_dot * T /(2*pi) ) vs cos(theta). y intercept should
%enable a calculation of eps
y = sqrt( theta_dot * T /(2*pi) ); x = cos(theta);
[yfit,xfit,r,m,c,dm,dc,yupper,ylower,s] = bestfitc(x,y);
ecc_from_bestfit = sqrt( 1- c^(-4/3) );
plot(xfit,yfit,'b-',x,y,'r+'); grid on; 
xlabel( 'cos\theta'); ylabel('sqrt( d\theta/dt * T/(2*pi) ');
set(gca,'fontsize',16); 
title(['T=',num2str(T),' yr, a=', num2str(a),...
    ' AU, orbital eccentricity=',num2str(ecc_from_bestfit)]);
print(gcf,'compute eccentricity.png','-dpng','-r300'); close(gcf);

%%

%Numeric method to compute polar angle theta (rad) vs orbit time.
%T is the orbital period, and time vector t has the same units as T.
%ecc is the eccentricity and theta0 is the initial polar angle.
function theta = angle_vs_time( t, T, ecc, theta0 )

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

%Define array of polar angles for orbits. N is number of orbits.
N = ceil( t(end)/T );
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 = T*( (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' );

%%

%Line of best fit function yfit = m*x +c, with product moment correlation
%coefficient r
function [yfit,xfit,r,m,c,dm,dc,yupper,ylower,s] = bestfitc(x,y)

%Find any x or y values that are NaN or Inf
ignore = isnan(abs(x)) | isnan(abs(y)) | isinf(abs(x)) | isinf(abs(y));
x(ignore) = []; y(ignore) = [];
n = length(x);

%Compute line of best fit
xbar = mean(x); ybar = mean(y); xybar = mean(x.*y);
xxbar = mean(x.^2 ); yybar = mean(y.^2 );
Vx = xxbar - xbar^2; Vy = yybar - ybar^2;
COVxy = xybar - xbar*ybar;
m = COVxy/Vx; c = ybar - m*xbar; r = COVxy/sqrt( Vx*Vy );
[x,i] = sort(x); y = y(i);
yfit = m*x + c; xfit = x;

%Compute errors in gradient m and intercept c
s = sqrt( (1/(n-2))*sum( (y - yfit).^2 ) );
dm = s/sqrt(n*Vx);
dc = ( s/sqrt(n) )*sqrt( 1 + (xbar^2)/Vx );

%Determine envelope lines
yupper = (m+dm)*x + c - dc;
ylower = (m-dm)*x + c + dc;

%End of code
