clear all
close all

T0 = 1200; %initial temperature of the lava
T_air = 25;
T_ground = T_air;
T_threshold = 400;

% a1 = 10;        %initial concentration
% T1 = 500;        %Initial Temperature of concentration a1
% T2 = 600;        %second temperature 
% A = (a1^(T1(T1-T2))); %constant value 
% alpha = log(a1)/(T1-T2); %the other constant value

%a = A*exp(alpha*T1); %concentration gas inside the rock as a function of tempereture 
% A = 10^6;
% alpha = 0.023;

% D1/a^2   %frequency factor
% D1       %diffusion coefficient
% a=       %lenght of crystal 
% Ea=      %activation energy 1MeV ~ 1.61 * 10^-13 Joules 
% R =      %gas constant 
% T =      %what temperature? 
% 
% D(T)/a^2 = (D1/a^a)*exp(-Ea/RT);

p = 2.65; %density of rock
k = (0.031)/(1.1*2.65); %diffusivity constant (heat conductance /heat capacity of the lava times the density of the host rock)
c = 1.1;   %The heat capacity of the lava in Joule/gramm/Kelvin

N=100;  %number of discritized intervals 
L=10;   %lenght my x axis in meters 
step = L/N;  %0.100 step
tmax = 5000;  %time but can change

M = 1; %cooling length scale

x = -M*L:step:L; %goes from -10 metets to +10 meters in steps of 0.100

T = zeros(size(x)); %matrix of the size of x
Tnew = zeros(size(T)); % memory allocation only once

dx = x(2)-x(1);  %leght between steps (this code 0.1000)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INITIAL CONDITION
real_time = zeros(1+tmax,1);

T(1:M*N+1) = T_ground;
T(M*N+2:M*N+1+N) = T0;

eta = 0.5;

dt = eta*(dx^2)/k; 

%(dt*k/dx^2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TIME ITERATIONS

D_pert = zeros(tmax,1);

a = zeros(2:M*N+1,tmax);

f = 0;
for t = 1:tmax %time count index
  
  for i = 2:M*N+N % PDE valid only inside the domain 
      Tnew(i) =  T(i) + eta*(T(i+1)-2*T(i)+T(i-1));
  end
  
  
  %BOUNDARY CONDITIONS (DIRICHLET TYPE) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % at x = +L   upper surface (air)
  Tnew(M*N+1+N) = T_air;
    
  % at x = -M*L lower surface (bedrock, basement)
  Tnew(1) = T_ground;
  
%   for i = 1:M*N+1;
%          a(i,t) = A*exp(alpha*T(i)); %concentration gas inside the rock as a function of tempereture 
%   end 
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   plot(a(:,t)T 'b')
% %   axis (-M*L L);
%   title('concentration',14);
%   xlabel('ELEVATION [m]', 'Fontsize',12); 
%   ylabel('TEMPERATURE [C]', 'Fontsize',12);
%   grid on;
  
%   for i = 1:M*N+1;
%       a(i,t) = (D1/a^2)*exp(-Ea/RT);
%       plot 
  %  title('gas diffusion as temp changes',14);
  %   xlabel('ELEVATION [m]', 'Fontsize',12); 
  %   ylabel('TEMPERATURE [C]', 'Fontsize',12);
  %   grid on;
      
      
      
      
      
  plot(x,T,'b',x,T,'b.');
  axis([-M*L L 0 1.1*T0]); 
  title('COOLING LAVA FLOW','Fontsize',14);
  xlabel('ELEVATION [m]', 'Fontsize',12); 
  ylabel('TEMPERATURE [C]', 'Fontsize',12);
  grid on;
  
  %this loop will stop the iteration at time 1, then press H to continue
  %iterations 
  
  if t==1
    fprintf('  real_time: %f  mean Temp: %f \n', real_time(t), mean(T));  
    pause 
  end
  
  pause(0.0005);
 
    
  real_time(t+1) = real_time(t) + dt;
  T = Tnew; %%%%% UPDATE THE SOLUTION 
  
  fprintf('  real_time: %f  mean Temp: %f \n', real_time(t), mean(T));
  
  %%%%%CAPTURING THE DEPTH OF THRESHOLD TEMPERATURE
  
  for i = 1:M*N+1+N
    
      if T(i) >= T_threshold
          
        D_pert(t) = x(i);
        
        break;
        
      end
      
  end
  
  if mod(t+49,50)==0
            f=f+1;
         
%       surf(x,y,T)
     axis([-M*L L 0 1.1*T0])
      saveas(gcf, ['andresprojmov_',sprintf('%03d',f),'.png'])
    pause(0.001)
   end
  
 end


figure
plot(real_time(2:tmax+1),D_pert,'k');
title('DEPTH OF TEMPREATURE PERTURBATION','Fontsize',14);
xlabel('TIME', 'Fontsize',12); 
ylabel('ELEVATION', 'Fontsize',12);
grid on;