%%% Imitation of Massive Planet-Asteroid Crater Transformation
%%% Author : Matthew Mansat
%%% EPS 109 Fall 2015 Final Project

clear all;
close all;

% Initial conditions
f = 0;
explode = 0;
N = 200; %km wide crater
L = 200; %
T = zeros(N,N);

x = 0:L/N:L;
dx = L/N;
dy = dx;
k = 1.3; %constant k
dt = 0.2; %time step dt
d = 2;
Nhalf = (N*0.5);
    
T(:,1) = 0;
T(1,:) = 0;
T(N,:) = 0;
T(:,N) = 0;    
    
for i =Nhalf-5:Nhalf+5 %Constructs 'meteor' N distance above ground
	for j = Nhalf-5:Nhalf+5               
		T1(i,j) = d * (i - 1);
		T2(i,j) = d * (N + 1 - i);
		T3(i,j) = d * (j - 1);
		T4(i,j) = d * (N + 1 - j);
    T(i,j) = min([T1(i,j) , T2(i,j) , T3(i,j) , T4(i,j)]);
    end
end



for iter = 1:5000
    Tnew = T;
      
    % Initial point of impact
    if explode < 1
        for i = 2:N-1
            for j = 2:N-1
                Tnew(i,j) = (k * (dt/dx^2) * (( T(i+1,j) - (2 * T(i,j)) + T(i-1,j) ) ) ) + (k * (dt/dy^2) * (( T(i,j+1) - (2 * T(i,j)) + T(i,j-1) ) ) ) + T(i,j);
            end
        end
    end
    
     
    T = Tnew;
    xx = 1:N;
    yy = 1:N;
    [X,Y] = meshgrid(yy,xx);
    
        %surf(X,Y,T);
        %shading flat;
        %if T(Nhalf, Nhalf) > 0 && explode < 1
            %colormap summer;
        %end
        %if T(Nhalf,Nhalf) < 0
            %colormap autumn;
        %end
        %view([+34.5 14]);
        %axis tight
            %xlabel('x-axis (km)');
            %ylabel('y-axis (km)');
            %zlabel('height (km)');
            %title(sprintf('t = %6.2f seconds', dt*iter));
        %axis([0 N 0 N -N 2*N]);
    
    whitebg('black');
    
    
    
    if mod(dt*iter,0.1) == 0
        f = f + 1;
        surf(X,Y,T);
        
        shading flat;
        if T(Nhalf, Nhalf) > 0 && explode < 1
            colormap summer;
        end
        if T(Nhalf,Nhalf) < 0
            colormap autumn;
        end
        view([+34.5 14]);
        axis tight
            xlabel('x-axis (km)');
            ylabel('y-axis (km)');
            zlabel('height (km)');
            title(sprintf('t = %6.2f seconds', dt*iter));
        axis([0 N 0 N -N 2*N]);
        saveas(gcf, ['impact_',sprintf('%03d',f),'.png']);
        pause(0.01)
    end
    
    if T(Nhalf,Nhalf) < -200
        explode = explode + 1; %event has exploded
        saveas(gcf, ['impact_',sprintf('%03d',f),'.png']);
        break;
    end
 
    pause(0.00005);
end