clf; close all; clc;
%% Rock Properties and Initial Conditions
% Limestone, Sandstone, and  Slate Thermal Diffusivity UNITS cm^2/s
lmstnK = 12.18*10^-3; sstnK = 16.45*10^-3; sltK = 9.26*10^-3;         
% Temperatures
IDykeTemp = 1200;               % Tdyke(0,x) UNITS Celcius
IRockTemp = 0;                  % Trock(0,x) UNITS Celcius

%% Create 1D Initial Temperature Profile
N = 250; L = 2500; X=0:L/N:L;                                      
T = zeros(1,N+1);
T(1:N*2/5) = IRockTemp; T((N*3/5)+1:N) = IRockTemp; % Set T(0,x) Rock
T(((N*2/5)+1)+1:N*3/5) = IDykeTemp;                 % Set T(0,x) Dyke

%% Create Cooling Equation Variables
DELTAx = L/N;                           % spacial resolution UNITS cm
ETAc = 0.5;                             % critical eta
DELTAtc = ETAc*DELTAx^2/lmstnK;         % critical timestep
DELTAt = DELTAtc*0.99;                  % calculate large timstep
ETA = lmstnK*(DELTAt/DELTAx^2);         % calculate eta

%% Cooling Equation
tyear = 3600*24*365;            % no. of seconds in a year
iteryear = tyear/DELTAt;        % corresponding no. of iterations
Tnew = T;
for t = 1:iteryear
    for i = 2:N
        Tnew(i) = lmstnK*(DELTAt/DELTAx^2)*(T(i+1)-2*T(i)+T(i-1))+T(i);
    end 
    T = Tnew;
end 

maxT = max(T);

plot(X,T,'k-')
hold on
plot([0 1250],[maxT maxT],'r-')
title('1D Temperature Profile in Limestone After One Year') ; 
xlabel('x / centimeters'); ylabel('Temperature / ^{o}C');
annotation('textbox',...
    [0.15178125 0.867156208277704 0.103125 0.0300400534045394],...
    'String',{'Maximum Temperature =' num2str(maxT)},'fitboxtotext','on');
axis([0 L 0 1250])
