%Kepler_orbit_discovery
% Demo which reverse-engineers from exact elliptical orbits of a planet
% about the Sun (e.g. Mars), the polar angles of a 'home' planet from a
% star (e.g. Earth), and the polar angle of the planet from the home
% planet. Let's assume these are the angles Kepler will have obtained from
% Tycho Brahe, vs time.
%
% In the first instance, let's assume Earth has a perfectly circular orbit
% about the Sun, with Earth-Sun distance (1AU) known from the 'Cosmic
% ladder' of Earth radius, Earth-Moon distance, Moon radius, Solar radius
% and 1AU calculations known (or 'potentially knowable') with methods from
% antiquity (e.g. Aristarcus). We also will assume we know the period of
% the orbits of (i) the Earth and (ii) the planet.From the 'generated'
% Kepler angle data, we will solve to find the orbit of the planet, and
% compare to the elliptical solution we started from.
%
% Once the idea of elliptical orbits is understood (Kepler's first law), we can use
% Kepler's second and third laws to determine the eccentricity of the home
% planet from solar angle data, and its time derivative. The orbit of the
% home planet can then be used to determine a more exact orbit of the other
% planet.
%
% Dr Andy French. March 2025.

function Kepler_orbit_discovery

%%% INPUT DATA FROM MODERN MEASUREMENTS %%%

%Star mass (in solar masses)
M = 1.23;  

%Home and 'away' planet orbital periods (in years)
P_H = 3.21; P = 6.54;

a_H = ( M^(1/3) )*( P_H^(2/3) ); %Home planet semi-major axis (AU)
ecc_H = 0.5;  %Home planet orbital eccentricity

a = ( M^(1/3) )*( P^(2/3) ); %Planet semi-major axis (AU)
ecc = 0.8; %Planet orbital eccentricity

%Initial polar orbital angle (rad) of planet and home planet from initial Earth-Sun vector
theta0_H = 0; theta0 = pi/6;

%

%%% SIMULATION INPUTS %%%

tmax = 40*P; %Data collection time (in years)
dt = 1/52;  %Data collection time interval (years)
fsize = 18;  %Fontsize for graphs

%

%Determine orbits of Home planet (H) and planet (P). Scale is in AU.
t = 0: dt : tmax; N = length(t);

%Assume circular Home planet orbit (only correct is ecc_H=0)
theta_H_c = 2*pi*t + theta0_H; x_H_c = a_H*cos(theta_H_c); y_H_c = a_H*sin(theta_H_c);
theta_H_c = rem(theta_H_c,2*pi);  %Map theta to range [0,2*pi] radians

%Assume elliptical Home planet orbit
theta_H = angle_vs_time( t, P_H, ecc_H, theta0_H );
theta_H = rem(theta_H,2*pi);  %Map theta to range [0,2*pi] radians
r_H = a_H*(1-ecc_H^2)./( 1-ecc_H*cos(theta_H) );
x_H = r_H.*cos(theta_H); y_H = r_H.*sin(theta_H);

%Assume elliptical planet orbit
theta_P = angle_vs_time( t, P, ecc, theta0 );
r_P = a*(1-ecc^2)./( 1-ecc*cos(theta_P) );
x_P = r_P.*cos(theta_P); y_P = r_P.*sin(theta_P);

%

%% CIRCULAR HOME PLANET MODEL ORBIT CALCULATION FOR OTHER PLANET %%

%Determine polar angles (rad) of planet from earth. This is what could have
%been measured before the time of Kepler.
phi_c = atan2( y_P - y_H_c, x_P - x_H_c );

%Step through time, angle data and determine planet (x,y) coordinates at
%intervals of one planet period.
dn = round( N*P/tmax ); %Data array elements corresponding to one planet period
x_c = []; y_c= []; tP_c = []; %Initialize planet position (AU) and time (years) data.
for n=1:N
    if (n+dn)<=N
        theta1 = theta_H_c(n); theta2 = theta_H_c(n+dn);  %Home planet angles
        phi1 = phi_c(n); phi2 =  phi_c(n+dn); %Planet angles
        
        %Solve for planet (x,y) coordinates
        k = tan(phi2) - tan(phi1) ;
        if k > 0.1  %Calculation can become unstable for small k  
            xx = sin(theta1) - sin(theta2) - tan(phi1)*cos(theta1) + tan(phi2)*cos(theta2);
            xx = xx/k;
            yy = xx*tan(phi1) - tan(phi1)*cos(theta1)+sin(theta1);
            x_c = [x_c,xx]; y_c = [y_c,yy];
            tP_c = [tP_c,t(n)];
        end
    end
end

%

%% DETERMINE ECCENTRICITY OF HOME PLANET

%Note we know the orbital period, so we can determine (from Kepler III),
%the semi-major axis a_H from this (as long as we know the mass of the star
%in solar masses).

%Calculate rate of change of angle (rad/year). This is the exact rate of
%change. Note this could be determined numerically from the table of
%theta_H vs time. i.e. theta_H(n) = ( theta_H(n) - theta_H(n-1) )/dt
theta_dot_H = (2*pi/P_H)* ( (1-ecc_H^2)^(-3/2) ) .*(1-ecc_H*cos(theta_H)).^2 ;

%Plot sqrt( theta_dot * P /(2*pi) ) vs cos(theta). y intercept should
%enable a calculation of eccentricity
y = sqrt( theta_dot_H * P_H /(2*pi) ); x = cos(theta_H);
[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.','markersize',10); grid on; 
xlabel( 'cos(\theta)'); ylabel('sqrt( d\theta/dt * T/(2*pi) ) ');
set(gca,'fontsize',18); 
title(['P=',num2str(P_H),' yr, a=', num2str(a),...
    ' AU, orbital eccentricity=',num2str(ecc_from_bestfit)]);
print(gcf,'Home planet eccentricity.png','-dpng','-r300'); close(gcf);

%

%% ELLIPTICAL HOME PLANET MODEL ORBIT CALCULATION FOR OTHER PLANET %%

%Determine polar angles (rad) of planet from earth. This is what could have
%been measured before the time of Kepler.
phi = atan2( y_P - y_H, x_P - x_H );

%Step through time, angle data and determine planet (x,y) coordinates at
%intervals of one planet period.
dn = round( N*P/tmax ); %Data array elements corresponding to one planet period
x = []; y= []; tP = []; %Initialize planet position (AU) and time (years) data.
for n=1:N
    if (n+dn)<=N

        %Home planet angles
        theta1 = theta_H(n); theta2 = theta_H(n+dn);
        
        %Home planet range
        r1 = a_H*(1-ecc_H^2)/(1 - ecc_H*cos(theta1) );
        r2 = a_H*(1-ecc_H^2)/(1 - ecc_H*cos(theta2) );
        
        %Home planet X,Y coordinates
        X1 = r1*cos(theta1); Y1 = r1*sin(theta1);
        X2 = r2*cos(theta2); Y2 = r2*sin(theta2);
        
        %Planet angles
        phi1 = phi(n); phi2 =  phi(n+dn); 
        
        %Solve for planet (x,y) coordinates
        k = tan(phi2) - tan(phi1) ;
        if k > 0.1  %Calculation can become unstable for small k  
            xx = Y1 - Y2 + X2*tan(phi2) - X1*tan(phi1) ;
            xx = xx/k;
            yy = Y1 + (xx-X1)*tan(phi1);
            x = [x,xx]; y = [y,yy];
            tP = [tP,t(n)];
        end
    end
end

%

%% PLOT ORBITS %%

%Create figure with white axes and black background
fig = figure('name','solar system','color',[0 0 0],'InvertHardCopy', 'off');
axes('nextplot','add','fontsize',fsize,...
    'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1],'zcolor',[1 1 1]);
axis equal; grid on; box on;

%Plot Sun
plot(0,0,'oy','markersize',10,'markerfacecolor','y');
pS = plot(0,0,'y*','markersize',30);

%Plot Home planet and planet orbits
pH = plot(x_H,y_H,'b-','linewidth',2); pP = plot(x_P,y_P,'r-','linewidth',2);
%pHc= plot(x_H_c,y_H_c,'c-','linewidth',1);
xlabel('x /AU'); ylabel('y /AU');
legend([pS,pH,pP],{'Star','Home','Planet'},...
    'fontsize',fsize,'textcolor',[1,1,1],'edgecolor',[1,1,1]);
ylim([-a,a]);
title(['P_H=',num2str(P_H),' Yr, a_H=',num2str(a_H,3),' AU, ecc_H=',num2str(ecc_H),...
    ', P=',num2str(P),' Yr, a=',num2str(a,3),' AU, ecc=',num2str(ecc)],...
    'fontsize',16,'color',[1,1,1]);

%Plot planet positions based upon this computation
plot(x,y,'r*');

%Print plot
print( gcf,'kepler orbit discovery.png','-dpng','-r300'); close(gcf);

%Plot angles vs time
figure; plot(t,theta_H,t,phi); set(gca,'fontsize', 18); xlabel('time /years'); ylabel('Angle /rad');
legend({'\theta_H','\phi'}); title('Home and planet polar angles'); grid on;
print( gcf,'kepler orbit discovery angles.png','-dpng','-r300'); close(gcf);

%

%% Export time and angle data to an Excel sheet %%
X{1,1} = ['Kepler angle data, assuming Home planet has an orbital period of ',num2str(P_H),' years.'];
X{3,1} = ['Orbital period of planet = ',num2str(P),' years. This is an interval of ',num2str(dn),' measurements.'];
X{5,1} = 'Time /years';
X{5,2} = 'Home planet polar angle /rad';
X{5,3} = 'Home planet theta_dot (rad/year)';
X{5,4} = 'Planet polar angle /rad from Home';
for n=1:length(t)
    X{5+n,1} = t(n);
    X{5+n,2} = theta_H(n);
    X{5+n,3} = theta_dot_H(n);
    X{5+n,4} = phi(n);
end
xlswrite('Kepler angle and time.xlsx',X);

%%

%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