%%%% Numerical Model of Freezing Pool
% Assumptions:
% Pool is fresh liquid water
% Liquid starts at uniform temperature
% Constant in X and Y directions (only varies in vertical direction)
% No heat enters system

% 0 Thermal energy defined as when water loses all latent heat and becomes
% ice

close
clear

% Initial conditions:
T0 = 0; % Bottom water temperature [°C]
TW = 0; % Initial temperature of water [°C]
T1 = -5; % Temperature of surface (air) [°C]
d = .1; % Depth of pool [m]
tmax = 1*3600; % Length of simulation [s]

% Material properties
pW = 1000;  % Density of water [kg/m^3]
pI = 916.7; % Density of ice [kg/m^3]
L = 334000; % Latent heat of fusion of water [J/kg]
cW = 4186; % Heat capacity of water [J/(kg?°C)]
cI = 2100; % Heat capacity of ice [J/(kg?°C)]
kW = 0.56; % Thermal conductivity of water [J/(m?s?°C)]
kI = 2; % Thermal conductivity of ice [J/(m?s?°C)]

% Boundary conditions:
% Bottom temperature constant
% Air temp constant

dz = 0.001; % Spatial discretization for liquid
dZ = dz*(pW/pI)^(1/3); % Spatial discretization for solid, conserving mass
dt = 0.002; % Temporal discretization
N = round(d/dz);
niter = round(tmax/dt); 

T = TW * ones(N+1,1); % Initial water temperature
Q0 = (L*pW + pW*cW*TW)*dz^3; % Assume cells are cubes with side length dz
Q = Q0 * ones(N+1,1); % Initial heat
yW = 1; % Uppermost liquid layer
yI = 0; % Lowermost solid layer
delt = NaN(1,N);

% Make colormap
cmap = [(1:-0.01:0)',(1:-0.01:0)',ones(101,1);0 0 0.5];
colormap(cmap)

count = 1;
f = 1;

% y positive down
for t = 1:niter
    % Heat flow 
    if yI ~= 0 % Case: Ice has formed, calculate temperature change 
        % separately for ice, water, and boundary layer
        T(yW) = max([(Q(yW)-L*pW*dz^3)/(pW*cW*dz^3),0]);
        DTI = dt*kI/(dZ^2*cI*pI) * ([T1;T(1:yI-1)] - 2*T(1:yI) + ...
            [T(2:yI);TW]);
        DTW = dt*kW/(dz^2*cW*pW) * (T(yW:end-1) - 2*T(yW+1:end) + ...
            [T(yW+2:end);T0]);
        Q(1:yI) = Q(1:yI) + pI*dZ^3*cI*DTI;
        Q(yW+1:end) = Q(yW+1:end) + pW*dz^3*cW*DTW;
        T(1:yI) = T(1:yI) + DTI;
        T(yW+1:end) = T(yW+1:end) + DTW;
    else % Case: Ice has not formed, calculate temperature change for water
        % only
        T(1:end) = dt*kW/(dz^2*cW*pW) * ([T1;T(1:end-1)] - 2*T(1:end) ...
            + [T(2:end);T0]) + T(1:end);
    end
    
    % heat conducted from uppermost liquid layer
    if yI == 0
        Q(yW) = Q(yW) -kW*dz^2 * (T(yW)-T1)/dz;
    elseif any(yW)
        Q(yW) = Q(yW) -kW*dz^2 * (T(yW)-T(yI))/dz; % Fourier heat equation
    end
    
    if any(Q>0)
        yW = find(Q>0,1,'first');
        if yW-1 ~= yI % Did boundary move downward?
            delt(count) = t*dt;
            count = count + 1;
        end
        yI = yW-1;
    else
        yI = N;
    end
    if mod(t,10000) == 0
        S1 = subplot(3,2,[1,3]);
        plot([T1,T0],-[yI, yI]*dz,'r-',T(:),0:-dz:-d,'LineWidth',2)
        axis([min([T0 TW T1]) max([T0 TW T1]) -d 0])
        S1.YTickLabel = linspace(d,0,11);
        title('Temperature')
        xlabel('Temperature [°C]')
        ylabel('Depth [m]')
        S2 = subplot(3,2,[2,4]);
        imagesc(T)
        title(['Time Elapsed: ', num2str(t*dt),' seconds'])
        S2.XTickLabelMode = 'manual';
        S2.YTickLabelMode = 'manual';
        cbar = colorbar;
        cbar.Label.String = 'Temperature [°C]';
        S3 = subplot(3,2,5:6);
        plot(delt,(1:N)*dz,'r-','LineWidth',2)
        axis([0 1500 0 d])
        title('Depth of Ice-Water Boundary')
        xlabel('Time [s]')
        ylabel('Boundary Depth [m]')
%         saveas(gcf,['stefan',num2str(f,'%03d'),'.png'])
        f = f+1;
        drawnow
    end
end