% The oscillation of a mass being dropped through a radial shaft through a
% solid planet of uniform density has the same period as an orbit about the
% planet at the radius of the planet. This fact is used to generate a
% simulation of a planet-sized gravity clock!
function gravity_clock
r = 1;  %Radius of masses (relative units, for screen)
R = 5;  %Radius of Planet (relative units, for screen)
rho = 5.972e24/( (4/3)*pi*(6371*1000)^3 );  %Density of planet /kgm^-3
G = 6.67e-11;  %Universal Gravitational Constant /Nm^2 kg^-2
dt = 10; %Timestep /s
fsize = 18; %Fontsize

%Calculate period of oscillations /s
P = sqrt(3*pi/(G*rho)); Pmin= P/60;

%Set up figure
fig = figure('color',[1 1 1],'name','Gravity Clock','units','normalized',...
    'position',[0.2 0.2 0.5 0.5],'renderer','painters');

%Draw initial setup
x1 = 0;  x2 = 0; y1 = R;  y2 = R;
theta = linspace(0,2*pi,1000);
xs = R*cos(theta); ys = R*sin(theta);
patch(xs,ys,'g','edgecolor','none');
hold on; axis equal;
p1 = plot(x1,y1,'r*','markersize',fsize*2);
p2 = plot(x2,y2,'ko','markersize',fsize*2); axis off;

%Set title
title( ['\rho = ',num2str(rho,4),' kgm^{-3}, Period = ',...
    num2str(Pmin,3),' minutes.'],'fontsize',fsize )

%Step simulation and update
stop = 0; t = 0;
while stop ==0
    %Update time
    t = t + dt;
    
    %Update mass being vertically dropped through radial shaft
    y1 = R*cos( 2*pi*t/P );
    
    %Update mass orbiting at the surface
    y2 = R*cos( 2*pi*t/P ); x2 = R*sin( 2*pi*t/P );
    
    %Update plot
    try
        set(p1,'ydata',y1);
        set(p2,'xdata',x2,'ydata',y2); drawnow;
    catch
        return
    end
end

%End of code