%% Final Project: Applying Monte Carlo in 3D to Heat Equation


% EPS109 Computer Simulations in Earth and Planetary Science, Fall 2012
% Professor Burkhard Militzer

% Antoine Bichara, Freshman

clc;
clear all;
clf;
close;


x = linspace(0,100,99); 
y = linspace(0,100,99); 


dx = 1; 
dy = 1;
f = 0;

ix = length(x); 
iy = length(y);

T = zeros(ix,iy);

% Initial Structure, could be modified to mirror a more defined landscape

T(15:25,5:45) = 3.5;
T(75:95,5:95) = 4.8;
T(25:75,50:70) = 4;
T(25:75,85:95) = 5;
T(5:10,30:45) = 4.5;
T(40:60,10:70) = 2;
T(60:65,20:40) = 3.7;
T(85:93,20:33) = 5;
T(5:15,80:95) = 4.3;

dt=0.7; % anywhere between 0.5 and 1 works well

kconstant = 0.25;
alpha = kconstant*dt/dx^2;

set(gca,'ydir','normal')
C = [];
riem2 = [];
error = [];

nn = 18000; % Number of points

for a=1:1000; % number of iterations
    
    %Tnew = T;
    
    for i = 2:98;
        
        for j = 2:98;
            
            T(i,j) = (alpha)*((T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1) - ...
                4*T(i,j))) + T(i,j);
           
            
        end
        
       
       Tnew = T + alpha*(T(i+1)-T(i));
       
     
        
    end
    
    % subplot(1,2,1)
%   hold on
%   c = 0;
%   rr = [];
%   rrr = rand(99,99)*5;
%   
%           for i = floor(rand(1,99)*99)+1;
%               for j = floor(rand(1,99)*99)+1;
%                   if rrr(i,j) > T(i,j)
%                      scatter3(x(i),y(i),rrr(i,j),'.') 
%                       c = c+1;
%                    
%                   else
%                       c = c;
%                   end
%               end
%           end
%           

h = subplot(2,3,1);
cla(h);
  
subplot(2,3,1)

points1 = rand(nn,2)*99;
points2 = rand(nn,2)*99;
points3 = rand(1,nn)*5;

w = rand(99,99)*5;

p1 = floor(points1)+1;
p2 = floor(points2)+1;

%if points1(floor(rand()*6000)) > T(

%scatter3(points(:,1),points(:,2),points1(:,1),'.'); 

c = 0;
qy = [];
qqy = [];
qn = [];
qqn = [];
wwy = [];
wwn = [];

for i = 1:nn
    if w(p1(i),p2(i))>T(p1(i),p2(i))
        
        qy(end+1) = p1(i);
        
        qqy(end+1) = p2(i);
        wwy(end+1) = w(p1(i),p2(i));
        
    else
        qn(end+1) = p1(i);
        qqn(end+1) = p2(i);
        
        
        wwn(end+1) = w(p1(i),p2(i));
        c = c+1;
    end
end
C(end+1) = c;

%scatter3(p(:,1),p(:,2),points1(:,1),'.');
scatter3(qy,qqy,points3(1:length(qy)),'.');


aaa = rand(6000,1)*5;
bbb = rand(6000,1)*5;

  [X,Y] = meshgrid(x,y);
  
  hold on
  
  subplot(2,3,1)
  surf(X,Y,T)
  view(3);
  axis([0 100 0 100 0 5]);
  drawnow
  title('Monte Carlo in 3D')
  
  subplot(2,3,3)
  pcolor(X,Y,T)
  title('Contour Plot')
  
  subplot(2,3,2)
  surf(X,Y,T)
   title('Melting Landscape')
  view(3);
  axis([0 100 0 100 0 5]);
  drawnow

  
  subplot(2,3,5)
  
  plot(C/nn*49005)
  set(gca,'Ytick',0:4000:24000,'YTickLabel',{'0','4000'...
      ,'8000','12000','16000','20000','24000'})
  axis([0 400 0 24000]);
 xlabel('# of Iterations')
 ylabel('Volume (m^3)')
 grid on
  
 title('Volume in m^3 as Function of Iterations')
 
 
 subplot(2,3,6)
 plot(100-(C/nn*49005/20123*100))
 set(gca,'Ytick',0:10:100,'YTickLabel',{'0','10','20','30','40','50','60','70','80','90','100'})
  axis([0 400 0 100]);
  grid on
  
 title('% of Volume Depleted as Function of Iterations')
 xlabel('# of Iterations')
 ylabel('% Depleted')
 
 
 subplot(2,3,4)
 format long
 riem1 = sum(T);
 
 riem2(end+1) = sum(riem1);
 format short
 error(end+1) = 100 * abs(riem2-(C(1:length(riem2))/nn*49005))/(C(1:length(riem2))/nn*49005);
 plot(error)
 
 
 set(gca,'Ytick',0:1:10,'YTickLabel',{'0','1','2','3','4','5','6','7','8','9','10'})
  axis([0 400 0 10]);
  grid on
  
 title('Relative Percentage Error With Riemann Sum')
 xlabel('# of Iterations')
 ylabel('% Error')
% set(gca,'Ytick',0:4000:24000,'YTickLabel',{'0','4000'...
%       ,'8000','12000','16000','20000','24000'})
%   axis([0 300 0 24000]);
%  grid on
%   
%  title('Riemann Sum over Domain of Temperature')
%  xlabel('# of Iterations')
%  ylabel('Volume m^3')
 
% subplot(2,3,6)

% 
% for n = n_iter
%     tic; fib_iter(n); ti = toc;
%     t_iter = [t_iter ti];
%     
%  
% end
% 
%    
% loglog(n_rec,t_rec,n_iter,t_iter)
% 
% xlabel('log N'), ylabel('log Time')
% legend('recursion','iteration')
% 
%  end
%  



end



