%Kepler_orbit_discovery_v1
% 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 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.
%
% Dr Andy French. February 2025.

function Kepler_orbit_discovery_v1

%%% INPUT DATA FROM MODERN MEASUREMENTS %%%

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

%Initial polar orbital angle (rad) of planet from initial Earth-Sun vector
theta0 = 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

%

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

%Assume circular Earth orbit
theta = 2*pi*t; x_E = cos(theta); y_E = sin(theta);
theta = rem(theta,2*pi);  %Map theta to range [0,2*pi] radians

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

%

%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]);

%

%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_E, x_P - x_E );

%Step through time, angle data and determine planet (x,y) coordinates at
%intervals of one planet period.
dn = round( N*T/tmax ); %Data array elements corresponding to one planet period
x = []; y= []; tP = 0; %Initialize planet position (AU) and time (years) data.
for n=1:N
    if (n+dn)<=N
        theta1 = theta(n); theta2 = theta(n+dn);  %Earth angles
        phi1 = phi(n); phi2 =  phi(n+dn); %Planet angles
        
        %Solve for planet (x,y) coordinates
        M = [-tan(phi1), 1; -tan(phi2), 1 ] ;
        if ( tan(phi2) - tan(phi1) ) > 1e-4   %Check matrix is not singular
            x_y_P = ( M^(-1) )*[ sin(theta1) - tan(phi1)*cos(theta1);...
                sin(theta2) - tan(phi2)*cos(theta2) ];
            x = [x,x_y_P(1)]; y = [y,x_y_P(2)];
            tP = [tP,t(n)];
        end
    end
end

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

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

figure;
plot(t,theta,t,phi);

%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