
clear 
clf

N=50;
L=50;
x=0:L/N:L;
T=zeros(N+1,N+1);
H=zeros(N+1,N+1);


Etamedium=0.01;
Etairon=.20408;
Etastone=.0043367;

maxi=N;
mini=1;
maxj=N;
minj=1;

f=0;
j=1;
i=1;
k=0;
hi=1;
hj=1;

T(1:15,1:N+1)=2200;
T(1:N+1,1:15)=2200;
T(35:N+1,1:N+1)=2200;
T(1:N+1,35:N+1)=2200;
T(16:34,16:34)=100;

medium=1;
stone=2;
iron=3;

dx=L/N;
dt=1*10^-4;

it_max=10000;


 for hi=1:N; 
  for hj=1:N;
   if hi>0
   if hi<N+1
   if hi>0
   if hi<N+1
        H(hi,hj) = medium;
   end
   end
   end
   end
   if hi<35
   if hi>15 
   if hj>15
   if hj<35
     r = rand()*2;
    if r>1
     H(hi,hj) = stone;
    end
    if r<1
     H(hi,hj) = iron;
    end
   end
   end
   end
   end
  end
 end


xx=0:L/N:L;
yy=0:L/N:L;
[X,Y]=meshgrid(xx,yy);

Tnew=T;

for k=1:it_max; 
T(1:15,1:N)=2200;
T(1:N,1:15)=2200;
T(35:N,1:N)=2200;
T(1:N,35:N)=2200;
      for i=2:N+1; 
       for j=2:N+1; 
            if H(i,j) == medium
              Tnew(i,j)=Etamedium*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1)-4*T(i,j))+T(i,j);
            end
            if H(i,j) == stone
              Tnew(i,j)=Etastone*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1)-4*T(i,j))+T(i,j);
                   if T(i,j)>1500
                    H(i,j)=medium;
                   end
            end
            if H(i,j) == iron
              Tnew(i,j)=Etairon*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1)-4*T(i,j))+T(i,j);
                   if T(i,j)>1800
                    H(i,j)=medium;
                   end
            end
       end
      end
 
      
      T=Tnew;
      
     
     
%      for j=minj
%          for i=1:N+1
%          Tnew(i,j)=T(i,j+1);
%          end
%      end
%      for i=mini
%          for j=1:N+1
%          Tnew(i,j)=T(i+1,j);
%          end
%      end
%      for j=maxj+1
%          for i=1:N+1
%          Tnew(i,j)=T(i,j-1);
%          end
%      end
%      for i=maxi+1
%          for j=1:N+1
%          Tnew(i,j)=T(i-1,j);
%          end
%      end
          
      
      
      %M=mod(k,1000);
      %if M==0
      
      if mod(k+99,100) == 0
          f=f+1;
          subplot(1,4,1)
          colormap(jet);
          surf(X,Y,T);
          axis([0 L 0 L -300 2200]); 
          %saveas(gcf, ['ablation',sprintf('%03d',f),'.png'])
          pause(1.00);
          
          subplot(1,4,2)
          colormap(jet);
          surf(X,Y,T);
          axis([0 L 0 L]); 
          %saveas(gcf, ['ablation',sprintf('%03d',f),'.png'])
          pause(1.00);
          
          subplot(1,4,3)
          colormap(jet);
          surf(X,Y,H);
          axis([0 L 0 L 0 3]);
          %saveas(gcf, ['materials',sprintf('%03d',f),'.png'])
          pause(1.00);
          
          subplot(1,4,4)
          colormap(jet);
          surf(X,Y,H);
          axis([0 L 0 L]); 
          saveas(gcf, ['ablation',sprintf('%03d',f),'.png'])
          pause(1.00);
          

      end
      
      
end

