%%%%%
% File Info
% Description: File used to create a basic simulation of an impact tsunami
% for educational use. Based off the paper "Tsunami" by SN Ward in 
% Encyclopedia of Solid Earth Geophysics, Springer Press 2010.
%
% Author: 
% J Naylor
% School of Aerospace, Mechanical and Mechatronic Engineering, Faculty of
% Engineering, University of Sydney, August 2019
%%%%%

%% Setup
clear all
close all
runfirst

%% Variable Definitions

% Wavelength Vector
lambda = 10e3:100:35e3;

% Wavenumber vector
k = 2*pi./lambda;

% Fluid depth
h = 4e3;

% Fraction of energy that goes into tsunami wave
epsilon = 0.15;

% Radius vector for solution
r1 = -60e3:100:60e3;

% Specific gravity of asteroid
SG = 4.5/0.997;

% Velocity
vi = 20e3;

% Impact radius
Ri = 100;

% Gravity
g = 9.81;

% Depth of cavity
Dc = (8*epsilon*SG*vi^2/(9*g))^0.25*Ri^0.75;

% Diameter of cavity
dc = 3*Dc;

% Radius of cavity
Rc = dc/2;

% Adjusted radius
Rd = sqrt(2)*Rc;


%% Plotting

% Increase fig size
fig = figure('Position', [10 10 900 800]);

% Set colour for ease of viewing
colormap winter
cmap = colormap;

% Restrict colourmap axis
caxis([3.5e3 4.8e3])

% Force white background
set(gcf, 'color', 'white');

% Iterate through time
for t = 1:60*10
    
    % Reset surface vector
    u_surf = [];
    
    % Solve for each radius
    for r = r1
        
        % Initial disturbance to zero
        usurf = 0;
        
        % Iterate through wave numbers
        for k1 = k
            
            % Solve for disturbance
            u_surf1 = 1/(2*pi)*k1*cos(sqrt(g*k1*tanh(k1*h))*t)*besselj(0,k1*r)*2*pi*Dc*Rd*besselj(3,k1*Rd)/k1*(k(2)-k(1));
            
            % Add to disturbance at that point
            usurf = usurf+u_surf1;
        end
        
        % Store solution in vector
        u_surf = [u_surf usurf];
    end
    
    % Setup plot
    [R PHI] = meshgrid(r1,linspace(0,2*pi,length(r1)));
    Z = ones(length(r1),length(r1));
    
    % Create Z matrix
    for i = 1:length(r1)
        Z(i,:) = (h+u_surf(i))*Z(i,:);
    end
    
    % Plot cross section
    ax1 = subplot(2,1,1);
    set(ax1,'Position',[0.13 0.6 0.775 0.3])
    plot(r1./1000,(h+u_surf)./1000)
    ylim([2.5 5])
    xlabel('$r$ (km)')
    ylabel('$z$ (km)')
    title(['Cross-section of Impact Wave @ t = ' num2str(t) 's'])
    drawnow
    
    % Plot 3D Surface
    ax2 = subplot(2,1,2);
    set(ax2,'Position',[0.13 .05 0.775 0.5])
    mesh(R.*cos(PHI)./1e3, R.*sin(PHI)./1e3, (Z./1e3)');
    zlim([2.5 5])
    xlim([-45 45])
    ylim([-45 45])
    caxis([3.6 4.4])
    xlabel('$x$ (km)')
    ylabel('$y$ (km)')
    zlabel('$z$ (km)')
    drawnow
    
    % Obtain image to generate GIF
    frame = getframe(fig);
    im{t} = frame2im(frame);
end

%% Save File
% Save as GIF
filename2 = 'waterEjecta2.gif';

% Iterate through number of frames
for idx2 = 1:600
    [A2,map2] = rgb2ind(im{idx2},256);
    
    % Write frames to file
    if idx2 == 1
        imwrite(A2,map2,filename2,'gif','LoopCount',Inf,'DelayTime',1e-36);
    else
        imwrite(A2,map2,filename2,'gif','WriteMode','append','DelayTime',1e-36);
    end
end