%Kepler_orbit_discovery_v2
% Demo which reverse-engineers from exact elliptical orbits of a planet
% about the Sun, the polar angles of the Earth from the Sun, and the polar
% angle of the planet from the Earth. These are the angles Kepler will have
% obtained from Tycho Brahe, vs time. 
%
% Let's assume when the polar angle of the Earth-Sun is zero, the Earth
% is 1AU from the Sun. This will scale our orbits.
%
% 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.
%
% To make the simulation more interesting, the eccentricity of BOTH the
% 'Earth' and the planet can be set (in the range 0 < ecc < 1). In other
% words, one doesn't need to assume the Earth's orbit is entirely circular.
%
% Dr Andy French. February 2025.

function Kepler_orbit_discovery_v2

%%% INPUT DATA FROM MODERN MEASUREMENTS %%%

a = 1.523;  %Planet semi-major axis (AU)
ecc = 0.09; %Planet orbital eccentricity
ecc_E = 0.4; %Earth orbital eccentricity
a_E = 1; % Earth semi-major axis
T = a^1.5;  %Planet period (in years), using Kepler III
T_E = a_E^1.5;  %Earth period (in years), using Kepler III

%Initial polar orbital angle (rad) of planet from initial Earth-Sun vector
theta_P0 = pi/6;

%

%%% SIMULATION INPUTS %%%

tmax = 10*T; %Data collection time (in years)
dt = T/50;  %Data collection time interval (years)
fsize = 18;  %Fontsize for graphs

%

%%% EXACT ORBIT CALCULATION %%%   (Well, assuming orbits in the same plane)

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

%Assume elliptical Earth orbit
theta_E = angle_vs_time( t, 1, ecc_E, 0 );
r_E = a_E*(1-ecc_E^2)./( 1-ecc_E*cos(theta_E) );
x_E = r_E.*cos(theta_E); y_E = r_E.*sin(theta_E);

%Assume elliptical planet orbit
theta_P = angle_vs_time( t, T, ecc, theta_P0 );
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);

%Determine polar angles (rad) of planet from earth. This is what could have
%been measured in the time of Brahe, Copernicus and Kepler.
phi = atan2( y_P - y_E, x_P - x_E );

%Determine polar angles (rad) of Earth from Sun. 
theta = atan2( y_E, x_E );

%

%Plot angles vs time for Earth and Planet
figure; plot(t,theta,t,phi); set(gca,'fontsize',fsize); xlabel('Time /years'); ylabel('Angle /rad');
title('Earth,Sun and Planet,Earth polar angles vs time');
legend({'Earth','Planet'},'fontsize',fsize); grid on;
print( gcf,'kepler orbit angles vs time.png','-dpng','-r300'); close(gcf);

%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 Earth and planet orbits
pE = plot(x_E,y_E,'b-','linewidth',2); pP = plot(x_P,y_P,'r-','linewidth',2);
xlabel('x /AU'); ylabel('y /AU');
legend([pS,pE,pP],{'Sun','Earth','Planet'},'fontsize',fsize,'textcolor',[1,1,1],'edgecolor',[1,1,1]);

%

%Step through time, angle data and determine Earth (X,Y) coordinates at
%intervals of one planet period. The corresponding array interval is dn.
dn = round( N*T/tmax ); X = 1 + ecc_E; Y = 0; tt = 0;
X(2) = tan( phi(1+dn) )/( tan( phi(1+dn) ) - tan( theta(1+dn) ) );
Y(2) = X(2)*tan( theta(1+dn) ); tt(2) = t(1+dn);
figure(fig); plot(X,Y,'b*');
for n=2: (N-dn-1)
    t_phi0 = tan( phi(1) );  t_phi_i = tan( phi(n) );  t_phi_i_plus1 = tan( phi( n+dn ) );
    t_th_i = tan( theta(n) );  t_th_i_plus1 = tan( theta( n+dn ) );
    X(n+1) = ( X(n)* ( ( t_phi_i -t_th_i - t_phi0 )/ ( t_phi_i - t_phi0 ) )*...
        ( t_phi_i_plus1 - t_phi0 ) + t_phi0 )/( t_phi_i_plus1 - t_th_i_plus1 );
    Y(n+1) = X(n+1)*tan( theta(n+1) );
    tt(n+1) = t(n+dn);
    plot(X(n+1),Y(n+1),'b*');
end

%Plot planet positions based upon this computation
plot(X,Y,'b*');

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

%Export time and angle data to an Excel sheet
X{1,1} = 'Kepler angle data, assuming Earth is a circular orbit of radius 1AU, taking a year';
X{3,1} = ['Orbital period of planet = ',num2str(T),' years. This is an interval of ',num2str(dn),' measurements.'];
X{5,1} = 'Time /years';
X{5,2} = 'Earth polar angle /rad';
X{5,3} = 'Planet polar angle /rad from Earth';
for n=1:length(t)
    X{5+n,1} = t(n);
    X{5+n,2} = theta(n);
    X{5+n,3} = 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' );

%End of code