clf; close all; clc;
%% Rock Properties and Initial T
% Limestone, Sandstone and Slate: Thermal Diffusivity UNITS m^2/s
lmstnK = 12.18*10^-3; sstnK = 16.45*10^-3; sltK = 9.26*10^-3;                        
% Initial Temperatures UNITS oC
IDykeTemp = 1200; IRockTemp = 0;              
% Metamorphism Temperature Thresholds UNITS oC
lmstnTmmax = 500; sstnTmmax = 400; sltTmmax = 300;

%% Create Initial 2D Temperature Profile
N = 50; L = 10; X = 0:L/N:L; Y = X; T = zeros(N+1,N+1);
T(1,1:N+1)=IRockTemp; T(1:N+1,1)=IRockTemp; T(1:N+1,N+1)=IRockTemp; 
T(N+1,1:N+1)=IRockTemp; T(1:51,(N/2+1)-N/10:(N/2+1)+N/10)=IDykeTemp;

%% Cooling Equation Variables
DELTAx = L/N;                           % Spacial Resolution UNITS m
ETAc = 0.25;                            % Critical ETA
DELTAtc = ETAc*DELTAx^2/lmstnK;           % Critical Timestep
DELTAt = DELTAtc*0.99;                  % Timestep UNITS s
ETA = lmstnK*(DELTAt/DELTAx^2);           % Thermal Diffusivity UNITS m^2/s

%% Metamorphism and Dyke Graphic
C = zeros(N+1,N+1); 
D = zeros(N+1,N+1); D(2:N,(N/2+1)-N/10:(N/2+1)+N/10) = 1;
F = 0;
%% 2D Cooling Equation and Plot
Tnew = T;
for t = 1:200
    for i = 2:N
        for j = 2:N
            Tnew(i,j) = (lmstnK*DELTAt/DELTAx^2)*(T(i+1,j)...
            +T(i-1,j)+T(i,j+1)+T(i,j-1)-4*T(i,j))+T(i,j);
            % Boundary Conditions
            Tnew(1,2:N) = T(2,2:N); Tnew(2:N,1) = T(2:N,2);
            Tnew(2:N,N+1) = T(2:N,N); Tnew(1,1) = 0.5*(T(2,1)+T(1,2));
            Tnew(1,N+1)=0.5*(T(2,N+1)+T(1,N));
            Tnew(N+1,1) = 0.5*(T(N,1)+T(N+1,2));
            Tnew(N+1,N+1) = 0.5*(T(N,N+1)+T(N+1,N)); Tnew(N+1,2:N) = 0;
            % Auerole
            if T(i,j) >= lmstnTmmax
                C(i,j) = C(i,j)+1/5;
                if C(i,j) >= 1
                    C(i,j) = 1;
                end 
            end 
        end 
    end
T = Tnew;
%if mod(t,2) == 0;
    % Plot Figure
    colormap jet
    pcolor(X,Y,T)
    shading interp
    hold on
    % Plot Heat Distribution
    contour(X,Y,T,[100 200 300 400 500 600 700 800 900 1000 1100 1200],'k');
    % Plot Dyke and Aureole
    [q,r] = contour(X,Y,D,1,'r'); set(r,'LineWidth',1); 
    [u,v] = contour(X,Y,C,1,'w'); set(v,'LineWidth',2);
    j = legend([r v],'Dyke','Metamorphic Zone'); 
    j.Box = 'Off'; j.TextColor = 'White';
    title('The Development of a Metamorphic Zone in Limestone');
    xlabel('x / meters'); ylabel('y / meters');
    axis([0 L 0 L 0 1200])
    drawnow
%     F = F + 1;
%     saveas(gcf, ['SlateMeta', num2str(F, '%03d'), '.png'])
%     pause(0.01)
%end 
end 