%tidal_g_field
% Model of the gravitational field inside and beyond a planet orbited by a
% moon. The planet also orbits a star, whose gravtational influence is also
% computed. The planet is assumed to be of uniform density and perfectly
% spherical. The differing field strength due to the interaction of the
% star and moon will contribute to tidal effects such as heating and
% bulges, in particular when liquids are present on the surface of the
% planet.
%
% LAST UPDATED by Andy French May 2017

function tidal_g_field
global d

%%% Astronomical parameters %%%

%Universal gravitational constant / Nm^2kg^-2
G = 6.67384e-11;

%Mass of sun /kg
M_s = 1.99e30;

%Earth sun mean distance /m
r_es = 1.49597871e11;

%Earth mass /kg
M_e = 5.972e24;

%Earth radius /m
R_e = 6.371e6;

%Earth moon mean distance /m
r_em = 0.3844e9;

%Moon mass /kg
M_m = 7.35e22;

%Moon Radius /m
R_m = 1737e3;

%Surface gravity on Earth /ms^-2
g_e = G*M_e/(R_e^2);

%

%%% Simulation specific parameters %%%

%Star, planet and moon masses /kg
M_star = 1e3*M_s;
M_planet = M_e;
M_moon = 1e5*M_m;

%Star to planet and planet to moon orbital distances /m
r_sp = r_es;
r_pm = r_em;

%Planet radius /m
R = R_e;

%Initial elevation angle of moon from planet-star line of sight /radians
d.theta0 = 0;

%Scaling factor /m
scale = R_e;

%Axis limits (in terms of scale)
xlimits = [-3,3];
ylimits = [-3,3];

%Number of field evaluation points in x and y directions
N = 200;

%Number of evaluation points for field direction arrows
NN = 40;

%Fontsize for graphs
fsize = 14;

%Moon elevation angle steps /radians for animation
dtheta = 2*pi/360;

%Frames per second for animation
FPS = 30;

%

%Calculate g field for map underlay and quiver plot
[g,gx,gy,x,y] = gcalc( G, M_star, M_planet, M_moon, d.theta0, R, r_sp, r_pm, scale*xlimits, scale*ylimits, N );
[gg,gxx,gyy,xx,yy] = gcalc( G, M_star, M_planet, M_moon, d.theta0, R, r_sp, r_pm, scale*xlimits, scale*ylimits, NN );

%Plot magnitude of g field
fig = figure('name','tidal_g_field','color',[1 1 1],'units','normalized',...
    'position',[0.2,0.2,0.5,0.5],'KeyPressFcn',...
    @keypressfunc,'renderer','opengl');
p = pcolor( x/scale,y/scale, g ); colormap('jet');
view(2);
shading interp;
interp_colormap(N);
colorbar;
caxis([0,max(max(g))]);
axis equal;
xlim(xlimits);
ylim(ylimits);
set(gca,'fontsize',fsize);
xlabel('x','fontsize',fsize);
ylabel('y','fontsize',fsize);
box on;
hold on;
title({'g field /Nkg^{-1} of planet, star, moon system',...
    ['M_s = ',num2str(M_star,3),'kg, M_p = ',...
    num2str(M_planet,3),'kg, M_m = ',num2str(M_moon,3),'kg'],...
    ['scale = ',num2str(scale,3),'m, r_{sp} = ',num2str(r_sp,3),'m, r_{pm} = ',num2str(r_pm,3),'m']},'fontsize',fsize)

%Overlay planet boundary
t = linspace(0,2*pi,N);
xp = R*cos(t);
yp = R*sin(t);
plot(xp/scale,yp/scale,'w-')

%Overlay g field arrows
q = quiver( xx/scale,yy/scale,gxx,gyy,'w' );

%%% Run simulation %%%
d.go = 1;
d.pause = 0;
d.t = 0;
while d.go == 1
    if d.pause==0
        
        %Update Moon elevation angle /radians
        d.theta0 = d.theta0 + dtheta;
        if d.theta0 > 2*pi
            d.theta0 = 0;
        end
        
        %Calculate g field for map underlay and quiver plot. Note we have to do this twice since the
		%plot density for the map (N) is much migher than the arrow density for the vector field.
        [g,gx,gy,x,y] = gcalc( G, M_star, M_planet, M_moon, d.theta0, R, r_sp, r_pm, scale*xlimits, scale*ylimits, N );
        [gg,gxx,gyy,xx,yy] = gcalc( G, M_star, M_planet, M_moon, d.theta0, R, r_sp, r_pm, scale*xlimits, scale*ylimits, NN );
        
        %Update plot
        set(p,'xdata',x/scale,'ydata',y/scale,'Cdata',g);
        set(q,'xdata',xx/scale,'ydata',yy/scale,'Udata',gxx,'Vdata',gyy);
        
        %Flush graphics buffer
        drawnow
    end
    pause(1/FPS)
end
close(fig);


%%

%g field calculator
function [g,gx,gy,x,y] = gcalc( G, M_star, M_planet, M_moon, theta0, R, r_sp, r_pm, xlimits, ylimits, N )

%Generate mesh of x,y coordinates centred on the centre of the planet
x = linspace( xlimits(1),xlimits(2),N );
y = linspace( ylimits(1),ylimits(2),N );
[x,y] = meshgrid(x,y);

%Determine mass of planet to use for each point
M = M_planet*ones(N,N);
inner = (x.^2 + y.^2) < R^2;
density = M_planet/( (4/3)*pi*R^3 );
M(inner) = density*(4/3)*pi*( ( x(inner).^2 + y(inner).^2 ).^(3/2) );

%Compute g field strength at x,y points due to planet, moon and star
gx_s = G*M_star*( -x - r_sp )./(( (-r_sp - x).^2 + y.^2 ).^(3/2));
gx_m = G*M_moon*( -x + r_pm*cos(theta0) )./(( ( -x + r_pm*cos(theta0) ).^2 + ( -y + r_pm*sin(theta0) ).^2 ).^(3/2));
gx_p = -G*M.*x./(x.^2 + y.^2).^(3/2);
gy_s = G*M_star*( -y )./( (-r_sp - x).^2 + y.^2 ).^(3/2);
gy_m = G*M_moon*( -y + r_pm*sin(theta0) )./(( ( -x + r_pm*cos(theta0) ).^2 + ( -y + r_pm*sin(theta0) ).^2 ).^(3/2));
gy_p = -G*M.*y./((x.^2 + y.^2).^(3/2));

%Determine magnitude of g field /ms^-2
gx = gx_s + gx_m +gx_p;
gy = gy_s + gy_m +gy_p;
g = sqrt( gx.^2 + gy.^2 );

%%

%interp_colormap
% Function which interpolates current colourmap to yield better graduated
% shading. N is number of possible colours.
function interp_colormap( N )

%Get current colourmap
map = colormap;

%Initialise new colormap
new_map = ones(N,3);

%Get size of current colormap and initalise red,green,blue vectors
dim = size(map);
R = ones(1,dim(1));
G = ones(1,dim(1));
B = ones(1,dim(1));
RR = ones(1,N);
GG = ones(1,N);
BB = ones(1,N);

%Populate these with current colormap
R(:) = map(:,1);
G(:) = map(:,2);
B(:) = map(:,3);

%Interpolate to yield new colour map
x = linspace( 1, dim(1), N );
RR = interp1( 1:dim(1), R, x );
GG = interp1( 1:dim(1), G, x );
BB = interp1( 1:dim(1), B, x );
new_map(:,1) = RR(:);
new_map(:,2) = GG(:);
new_map(:,3) = BB(:);

%Set colormap to be new map
colormap( new_map );

%%

%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,['g field theta=',strrep(num2str(d.theta0*180/pi),'.','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