%%% Solving the 1D Heat Equation in Real Time

clear all
close all

N = 50;
L = 20; % width of lava dike and adjacent host rock in meters

x = 0:L/(5*N):L;
dx = L/(5*N);

%{
d = 4; % width of lava dike in meters
dx = d/N; % descritization of lava dike
h = 2*d; % width of host rock adjacent to lava dike in meters
dh = h/(2*N); % descritization of adjacent host rock
%}

T(2*N) = 0;
T(3*N) = 0;

for i = 1:2*N;
    T(i) = 0;
end

for i = 3*N:(5*N)+1;
    T(i) = 0;
end

for i = (2*N)+1:(3*N)-1;
    T(i) = 1175;
end

plot(x,T,'ko-');

% temperature distribution
for r = 1:5000;
    Tnew = T;
    dt = 120; % change in time
    delta = (dt/(dx*dx));
    kappa = 0.025; % heat conduction coefficient, j/cm*seconds*kelvin
    cp = 1.05; % heat capacity, joule/gram/kelvin
    rho = 2.7; % density, gram/cubic centimeters
    k = (kappa/(cp*rho)); % coefficient in squared meters/second
    c = k*delta;
    for i = 2:N;
        Tnew(i) = c*(T(i+1)-2*T(i)+T(i-1))+T(i);
    end
    T = Tnew;
    plot(x,T,'ko-');
end