%tidesim
% Simulation of tidal bulges due to gravitational attraction of a spherical
% planet and a nearby moon. The planet and moon are in a mutually
% elliptical orbit. Concentric rings of essentially massless particles are
% set in circular motion about the planet. These represent liquid oceans on
% the surface. Assume a uniform covering of liquid initially i.e. no land
% masses! The massless particles feel the gravitational pull of both planet
% and moon, but have no mutual gravitational interaction.
%
% The elliptical orbits of the planet and moon are defined exactly,
% although a numeric integral determines time. A Verlet finite timestep
% solver is used to solve for the motion of the tidal masses.
%
% Key actions: 
%
% h    toggles hold simulation
% p    prints a PNG file of the current screen
% q    quit simulation
%
% LAST UPDATED by Andy French May 2017

function tidesim
global d

%%% Astronomical data %%%

%Universal Gravitational constant / m^3kg^-1s^-2
G = 6.67384e-11;

%Mass of Earth /kg
M_Earth = 5.977e24;

%Radius of Earth /m
R_Earth = 6371e3;

%Mass of Moon /kg
m_Moon = 7.35e22;

%Radius of Moon /m
R_Moon = 1737e3;

%Semi-major axis of Earth-Moon orbit /m
a_Earth_Moon = 0.3844e9;

%Eccentricity of Earth-Moon orbit /m
ecc_Earth_Moon = 0.0549;

%Rotation period of Earth /s
Trot_Earth = 24*3600;

%

%%% Astronomical parameters for this model %%%

M = M_Earth;
R = R_Earth;
m = M_Earth/2;
%m = m_Moon;
Rm = R_Moon;
a = R_Earth*2.5;
%a = a_Earth_Moon;
ecc = ecc_Earth_Moon;
%ecc = 0.3;


%Calculate orbital period /s using Kepler's Third Law
T = sqrt( ( (4*pi^2)/(G*(m+M)) )*a^3 );

%Calculate orbital period in Earth days
Tdays = T/(24*3600);

%Planet rotation period /s
Trot = T/5;

%

%%% Model parameters %%%

%Initial polar angle of Planet and Moon /radians
theta0 = 0;

%Define number of orbital values
N = 1000;

%Define number of rows of masslets
rows = 100;

%Define row spacing /m
row_spacing = R/100;

%Define arc spacing of masslets /m
arc_spacing = 2*pi*R/100;

%Fontsize for graphs
fsize = 18;

%Masslet marker size
msize = 1;

%Axis length scale /m
scale = a;

%Axis limits / scale units
xlimits = [-1 - ecc,1];
ylimits = [-1,1];

%Frame rate
FPS = 50;

%

%%% Define planet and moon orbit %%%%

%Determine timestep /s
dt = T/(N-1);

%Compute Planet, Moon orbit
[planet,moon,tt] = orbit( M,m,a,ecc,G,theta0,N );

%Interpolate at equi-spaced times
t = linspace(0,Tdays,N);
planet.x = interp1( tt, planet.x, t );
planet.vx = interp1( tt, planet.vx, t );
planet.y = interp1( tt, planet.y, t );
planet.vy = interp1( tt, planet.vy, t );
planet.t = t;
moon.x = interp1( tt, moon.x, t );
moon.vx = interp1( tt, moon.vx, t );
moon.y = interp1( tt, moon.y, t );
moon.vy = interp1( tt, moon.vy, t );
moon.t = t;
clear tt t

%Plot orbits on an x,y graph
fig = figure('color',[1 1 1],'name','Planet, Moon & tide simulation','KeyPressFcn',...
    @keypressfunc,...
    'renderer','opengl','units','normalized','position',[0.2,0.2,0.5,0.5] );
plot(planet.x/scale,planet.y/scale,'r--',moon.x/scale,moon.y/scale,'g--');
legend({'Planet','Moon'},'fontsize',fsize);
set(gca,'fontsize',fsize)
axis equal
hold on;
xlim(xlimits);
ylim(ylimits);
grid on;
xlabel('x /scale','fontsize',fsize)
ylabel('y / scale','fontsize',fsize)
d.titlestr = {'Planet, Moon & tide simulation',...
    ['scale = ',num2str(scale,3),' m, t = ',num2str(0),' days. T = ',num2str(Tdays,3),' days.'],...
    ['M = ',num2str(M,3), 'kg, m = ',num2str(m,3),' kg. T_{rot} = ',num2str(Trot/(24*3600),3),' days.']};
d.title = title( d.titlestr,'fontsize',fsize);

%

%%% Initial conditions for dynamic simulation %%%

%Determine initial positions and velocities of masslets
x = [];
y = [];
vx = [];
vy = [];
for n=1:rows
    %How many masslets in this row?
    current_radius = R + row_spacing*(n-1);
    nm = floor( 2*pi*current_radius/arc_spacing );
    
    %Determine initial position and velocity (circular orbits rotating with planet)
    angles = linspace( 0,2*pi, nm);
    for k=1 : nm
        x = [x,current_radius*cos(angles(k)) + planet.x(1)];
        y = [y,current_radius*sin(angles(k)) + planet.y(1)];
        vx = [vx,-current_radius*sin(angles(k))*2*pi/Trot + planet.vx(1)];
        vy = [vy,current_radius*cos(angles(k))*2*pi/Trot + planet.vy(1)] ;
    end
end

%Overlay planet and moon
theta = linspace(0,2*pi,500);
xp = R*cos(theta) + planet.x(1);
yp = R*sin(theta) + planet.y(1);
xm = Rm*cos(theta) + moon.x(1);
ym = Rm*sin(theta) + moon.y(1);
p_planet = patch( xp/scale, yp/scale, 'r','edgecolor','r' );
p_moon = patch( xm/scale, ym/scale, 'g','edgecolor','g' );

%Overlay line on planet to indicate rotation
rot_line = plot( [planet.x(1),planet.x(1)+R]/scale,...
    [planet.y(1),planet.y(1)]/scale,'k' );

%Add masslets to orbit plot
masslets = plot(x/scale,y/scale,'b.','markersize',msize);

%

%%% Run simulation %%%
d.go = 1;
d.pause = 0;
n = 1;
d.t = 0;
while d.go == 1
    if d.pause==0
        
        %Extract planet and moon positions
        Xp = planet.x(n);
        Yp = planet.y(n);
        Xm = moon.x(n);
        Ym = moon.y(n);
        
        %Determine x,y separations from planet to masslets
        dxp = x - Xp;
        dyp = y - Yp;
        
        %Determine x,y separations from moon to masslets
        dxm = -x + Xm;
        dym = -y + Ym;
        
        %Calculate x and y accelerations of masslets. This is the extra
        %gravitational aceleration towards the moon + the centripetal
        %acceleration towards the centre of the planet, to model the assumption
        %that the masslets are dragged around by the rotating planet. The
        %gravity of the planet is assumed to be balanced by the normal reaction
        %force.
        ax = G*m*dxm./( dxm.^2 + dym.^2 ).^(3/2) - ((2*pi/Trot)^2)*dxp;
        ay = G*m*dym./( dxm.^2 + dym.^2 ).^(3/2) - ((2*pi/Trot)^2)*dyp;
        
        %Calculate new masslet positions (Verlet method)
        x = x + vx*dt + 0.5*ax*dt^2;
        y = y + vy*dt + 0.5*ay*dt^2;
        
        %Determine x,y separations from planet to masslets
        dxp = x - Xp;
        dyp = y - Yp;
        
        %Determine x,y separations from moon to masslets
        dxm = -x + Xm;
        dym = -y + Ym;
        
        %Calculate updated x and y accelerations of masslets
        axx = G*m*dxm./( dxm.^2 + dym.^2 ).^(3/2) - ((2*pi/Trot)^2)*dxp;
        ayy = G*m*dym./( dxm.^2 + dym.^2 ).^(3/2) - ((2*pi/Trot)^2)*dyp;
        
        %Calculate new velocities of masslets (Verlet method)
        vx = vx + 0.5*dt*( ax + axx );
        vy = vy + 0.5*dt*( ay + ayy );
        
        %Determine whether any masslets have been 'eaten' by the moon or the earth
        %These will be removed from the simulation.
        eaten = ( ( dxm.^2 + dym.^2 ) < Rm^2 ) | ( ( dxp.^2 + dyp.^2 ) < R^2 );
        x(eaten) = [];
        y(eaten) = [];
        vx(eaten) = [];
        vy(eaten) = [];
        
        %Update time /days
        d.t = d.t + dt/(24*3600);
        
        %Update orbital step counter. Reset if beyond N (i.e. one orbit)
        n = n+1;
        if n>N
            n = 1;
        end
        
        %Update masslet, planet and moon positions
        xp = R*cos(theta) + planet.x(n);
        yp = R*sin(theta) + planet.y(n);
        xm = Rm*cos(theta) + moon.x(n);
        ym = Rm*sin(theta) + moon.y(n);
        set( p_planet,'xdata',xp/scale,'ydata',yp/scale );
        set( p_moon,'xdata',xm/scale,'ydata',ym/scale );
        set(masslets,'xdata',x/scale,'ydata',y/scale);
        
        %Overlay line on planet to indicate rotation
        set( rot_line, 'xdata', [planet.x(n),planet.x(n)+R*cos( 2*pi*d.t*24*3600/Trot)]/scale,...
            'ydata',[planet.y(n),planet.y(n)+R*sin( 2*pi*d.t*24*3600/Trot)]/scale );
        
        %Update title
        d.titlestr = {'Planet, Moon & tide simulation',...
            ['scale = ',num2str(scale,3),' m, t = ',num2str(d.t,3),' days. T = ',num2str(Tdays,3),' days.'],...
            ['M = ',num2str(M,3), 'kg, m = ',num2str(m,3),' kg. T_{rot} = ',num2str(Trot/(24*3600),3),' days.']};
        set(d.title,'string', d.titlestr)
        
        %Flush graphics buffer and pause before next frame
        drawnow
    end
    pause(1/FPS)
end
close(fig);

%%

%Binary system eliptical orbit function. Outputs are:
% planet.x
% planet.y
% planet.vx
% planet.vy
%
% and the same for Moon. Coordinate origin is barycentre of Planet, Moon
% system. t is time /days and N is the number of data points per orbit.
function [planet,moon,t] = orbit( M,m,a,ecc,G,theta0,N )

%Polar angle /radians
theta = linspace(theta0,theta0+2*pi,N);

%Time /days
t = sqrt( ( ( a*(1-ecc^2) )^3 )/( G*(m+M) ) ) * thetafunc( theta0, theta, ecc );
t = t/(24*3600);

%Define separation r /m
r = a*(1-ecc^2)./( 1-ecc*cos(theta) );

%x and y coordinates /m
planet.x = (m/(M+m))*r.*cos(theta);
moon.x = -(M/(M+m))*r.*cos(theta);
planet.y = (m/(M+m))*r.*sin(theta);
moon.y = -(M/(M+m))*r.*sin(theta);

%Velocity /ms^-1
Vx = (1-ecc*cos(theta))*sqrt( G*(m+M)/((1-ecc^2)*a) ).*...
    ( -sin(theta) - ( ecc*sin(theta)./(1 - ecc*cos(theta) )).*cos(theta) );
Vy = (1-ecc*cos(theta))*sqrt( G*(m+M)/((1-ecc^2)*a) ).*...
    ( cos(theta) - ( ecc*sin(theta)./(1 - ecc*cos(theta) )).*sin(theta) );
planet.vx = (m/(M+m))*Vx;
planet.vy = (m/(M+m))*Vy;
moon.vx = -(M/(M+m))*Vx;
moon.vy = -(M/(M+m))*Vy;

%%

%Eliptic integral function
function f = thetafunc( theta0, theta, ecc )
y = 1./( 1 - ecc*cos(theta) ).^2;
f = integrator( theta, y, theta0, theta );

%%

%integrator
% Integration function using cubic splines.
%
% Syntax: I = integrator( xx, yy, x0, x )
%
% xx, yy     Samples from function y(x) which is being integrated. These
%            must be the same dimensions and have at least three elements each.
% x0         Starting point (i.e. lower limit) of the integral
% x          Upper limit of the integral. This can be an array of any
%            dimensions.
%
% I          Output of integral corresponindg to array x
%
% The function is based on the code written on pages 292 and 294
% of "Mastering Matlab 6" by Duane Hanselman & Bruce Littlefield.

function I = integrator( xx, yy, x0, x )

%Check inputs are the right size
if ( numel(xx) == numel(yy) ) && ( numel(xx) >= 3 )
    
    %Sort xx in ascending order
    [xx,i] = sort(xx);
    yy = yy(i);
    
    %Compute piecewise polynomial based upon cubic splines.
    pp=spline(xx,yy);
    
    %Take apart piecewise polynomial
    [br,co,npy,nco] = unmkpp(pp);
    
    %Scale factors for integration
    sf=nco:-1:1;
    
    %Integral coefficients
    ico=[co./sf(ones(npy,1),:) zeros(npy,1)];
    
    %Integral spline has higher order
    nco=nco+1;
    
    %Set integration constant of all splines to be zero
    ico(1,nco)=0;
    
    %Find constant terms in polynomials
    for k=2:npy
        ico(k,nco)=polyval(ico(k-1,:),br(k)-br(k-1));
    end
    
    %Build pp form for integral
    ppi=mkpp(br,ico);
    
    %Evaluate definite integral between limits of x0 to x
    I = ppval(ppi,x) - ppval(ppi,x0);
    
else
    I = NaN(size(x));
end

%%

%Key press handler
function keypressfunc( fig,evnt )

%Allow read-write usage of global variable d
global d

%Get the string of the current key press
k = get( fig,'currentkey');

%Actions, depending on the key pressed...
if strcmp(k,'p') == 1
    %Print PNG image
    print( gcf,['tide t=',strrep(num2str(d.t),'.','p'),'.png'] ,'-dpng','-r300' );
elseif strcmp(k,'q') == 1
    %Stop simulation
    d.go = 0;
elseif strcmp(k,'h') == 1
    if d.pause == 0
        d.pause = 1;
    else
        d.pause = 0;
    end
end

%End of code
