%Ole Romer's calculation of the speed of light by recording the variations
%in the time of occultations of Io by Jupiter, as viewed from Earth.
%
% LAST UPDATED by Dr Andy French. September 2025

function ole_romer_io_speed_of_light

%Speed of light (m/s)
c = 2.998e8;

%Astronomical Unit (m) (approx mean Earth orbital radius about Sun)
AU = 1.496e11;

%Orbital period of Io /years about Jupiter
T_Io = 1.769/365.25;

%Period of Jupiter Orbit /years about Sun
T_J = 11.86;

%Radius of Jupiter Orbit (m) about Sun
r_J = 5.20*AU;

%Radius of Io Orbit (m) about Jupiter
r_Io = 421000*1e3;

%Time interval /years
dt = 1/(1000*365.25);

%Fontsize for graphs
fsize = 16;

%%

%Starting with Earth, Jupiter and Io in alignment, run time for one orbit
%of Jupiter. Calculate, w.r.t. 'absolute time', the times when Earth-Io and
%Earth-Jupiter vectors are within less than dtheta from each other.
dtheta = 2*pi/(360*3600);

t = 0 : dt : 1*T_J;  %Time vector in years
x_E = AU*cos(2*pi*t); y_E = AU*sin(2*pi*t);  %Earth x,y coordinates
x_J = r_J*cos(2*pi*t/T_J); y_J = r_J*sin(2*pi*t/T_J);  %Jupiter x,y corodinates
x_Io = x_J + r_Io*cos(2*pi*t/T_Io); y_Io = y_J + r_Io*sin(2*pi*t/T_Io);  %Io x,y coordinates

%Calculate angles between Earth Io and Earth Jupiter vectors and hence
%determine times of alignment. These are the times when Io is half way
%through an occultation. Ignore relative movement of Earth and Jupiter
%during one orbit of Io.
theta_E_Io = atan2( y_Io - y_E, x_Io - x_E );
theta_E_J = atan2( y_J - y_E, x_J - x_E );
a = find( abs( theta_E_Io - theta_E_J) < dtheta );
t_a = t(a);

%Determine Earth-Io distances at these times of alignment (m)
r_E_Io = sqrt( (x_Io(a)-x_E(a)).^2 + (y_Io(a)-y_E(a)).^2  );

%Determine time delay (s) between Io and Earth
t_delay = r_E_Io/c;

%Determine time (as recieved at Earth) of Jupiter-Io mid occultation
%(years)
t_a_E = t_a + t_delay/ (365.25*24*3600 );

%Plot t_a_E vs Earth-Io distance (in AU)
plot(r_E_Io/AU, t_a_E ); grid on; set(gca,'fontsize',fsize);
xlabel('Earth to Io distance /AU'); ylabel('Times / years');
title('Times of mid occultation of Io by Jupiter as viewed from Earth');
print(gcf,'Occultation time vs distance.png','-r300','-dpng'); close(gcf);

%Plot t_a_E vs Earth-Io distance (in AU)
plot(t_a, t_delay ); grid on; set(gca,'fontsize',fsize);
xlabel('Model times /years'); ylabel('Time delay /s');
title('Times of mid occultation of Io by Jupiter as viewed from Earth');
print(gcf,'Occultation time model vs measured.png','-r300','-dpng'); close(gcf);

%Plot t_delay vs Earth-Io distance (in AU)
plot(r_E_Io/AU, t_delay ); grid on; set(gca,'fontsize',fsize);
xlabel('Earth to Io distance /AU'); ylabel('Time delay /s');
title('Time difference (s) between predicted and measured');
print(gcf,'Occultation time delay vs distance.png','-r300','-dpng'); close(gcf);

%Plot Earth and Jupiter orbits
plot( x_E/AU,y_E/AU,'b-',x_J/AU,y_J/AU,'r-'); hold on; axis equal; grid on;
n = floor( 0.2*length(t) ); xlim([-6,6]); ylim([-6,6]);
plot( x_E(n)/AU,y_E(n)/AU,'b*',x_J(n)/AU,y_J(n)/AU,'r*');  set(gca,'fontsize',fsize);
xlabel('x / AU'); ylabel('y / AU'); title('Earth and Jupiter orbits');
print(gcf,'Earth and Jupiter orbits.png','-r300','-dpng'); close(gcf);

%Plot Earth and Jupiter spirograph
plot( x_E/AU,y_E/AU,'b-',x_J/AU,y_J/AU,'r-'); hold on; axis equal; grid on;
xlim([-6,6]); ylim([-6,6]); set(gca,'fontsize',fsize); nmax = 1000;
for n=1:nmax
    nn = ceil( n*length(t)/nmax );
    plot( [x_E(nn)/AU,x_J(nn)/AU],[y_E(nn)/AU,y_J(nn)/AU],'k-','linewidth',0.05);  
end
plot( x_E/AU,y_E/AU,'b-',x_J/AU,y_J/AU,'r-'); 
xlabel('x / AU'); ylabel('y / AU'); title('Earth and Jupiter orbits');
print(gcf,'Earth and Jupiter orbit spirograph.png','-r300','-dpng'); close(gcf);

%End of code