%Flooding the Planet

clear;

H = load('topography_180x360_grid.txt'); %Loads elevation numbers
size(H);
ny = size(H,1);
nx = size(H,2);
R = 6371000; %Radius of Earth in meters

f = 0;
L = 6000; % set higher than highest point in data set, in meters
ds = 10; %change in sea level
n = L/ds;

for i = 1:n
    area = 0;
    volO = 0; %ocean volume
    areaO = 0; %ocean area
    
%      if mod(n+9,1) == 0  
%          f = f+1;
%          saveas (gcf, ['project_', sprintf('%03d',f), '.png'])
%          pause(0.01)
%      end
    
    for iy = 1:ny
        lat = 90 - 180/ny*(iy-0.5); %latitude
    for ix = 1:nx
        H(iy,ix) = H(iy,ix) - ds; %changes the sea level
        
        long = 360/nx*(ix-1); %longitude
        dx = 2*pi*R/nx;
        dy = pi*R/ny;
        patch = dx*dy*cosd(lat);
        area = area + patch;
        
        if(H(iy,ix) < 0); %determines what points are below sea level
            areaO = areaO + patch;
            volO = volO + patch*abs(H(iy,ix));
        end
    end
    end
    
%     fprintf('%d area \n',areaO) %will print new ocean area
     fprintf('%d # of iterations \n',i)

    subplot(2,1,1);
	image(H)
    axis equal

    subplot(2,1,2);
    plot(i,areaO,'-.or')
    xlabel('iteration number')
    ylabel('Ocean area in sq. meters')
    title('New ocean area')
    hold on
    axis([0 n 3.6*10^14 5.5*10^14])
    drawnow
end