% all units measured in meters, kilograms, seconds, degrees Celsius
% define char of ice & water
c_i = 2050;
c_w = 4180;
a_i = 0.6;
a_w = 0.06;
dense_i = 917;
dense_w = 1025;

% solar constant in W/m2
solar = 1370;

% defining area
L=10;
N=200;
dx=L/N;
dy=L/N;
x=0:dx:L;
y=0:dy:L;

% temperature profile
T=zeros(N+1, N+1);
T(:,:)=17;
l=50;
T(N/2-l:N/2+l,N/2-l:N/2+l)=-10;

%ice profile
I=zeros(N+1,N+1);
I(:,1)=1;
I(:,N+1)=1;
I(1,:)=1;
I(N+1,:)=1;

% PDE conditions
dt=0.0001;
k=1;
eta = k*dt/dx/dx;
energy = zeros(N+1,N+1);
f=0;

for iter=1:4000
    Tnew=T;
    Inew=I;
    for i=2:N
        for j=2:N
            if T(i,j)<0
               albedo = a_i;
               mass = dense_i;
               specific_heat = c_i;
            else albedo = a_w;
                 mass = dense_w;
                 specific_heat = c_w;
            end
            energy(i,j) = solar*dt*(1-albedo);
            Tnew(i,j) = T(i,j)+(energy(i,j))/mass/specific_heat+eta*(T(i+1,j)-2*T(i,j)+T(i-1,j)+T(i,j+1)-2*T(i,j)+T(i,j-1));
             
            if Tnew(i,j)>0
               Inew(i,j)=1;
            end
             
        end
    end
    T=Tnew; 
    I=Inew;
    
    if mod(iter+39,40)==0
        f=f+1;
        imagesc(I);
        set(gca, 'YDir');
        colormap(1-gray);
        axis equal;
        axis tight;
        saveas(gcf, ['icearea_' num2str(f, '%03d'), '.png'])
        pause(0.001)
    end   
    

%     if mod(iter+39,40)==0
%         f=f+1;
%         [X,Y]=meshgrid(y,x);
%         surf(X,Y,T);
%         axis([0 L 0 L -10 20])
%         saveas(gcf, ['tempice_', num2str(f, '%03d'),'.png'])
%         pause(0.001)
%     end

  
end

