clc;
clear all;
clf;
close;


xxx = [];
errorav = [];
 ccc = 0;

for i = 1:200
    ccc = ccc+1
x = linspace(0,100,99); 
y = linspace(0,100,99); 

 nn = i*300;
 
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=1; % anywhere between 0.5 and 1 works well

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

set(gca,'ydir','normal')



   
C = [];
riem2 = [];

error = [];

%nn = 14000; % Number of points

    

for a=1:1:300; % 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
  


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);
  
  

 
 riem1 = sum(T);
 riem2(end+1) = sum(riem1);
 error(end+1) = abs(riem2-(C(1:length(riem2))/nn*49005))/(C(1:length(riem2))/nn*49005)*100;
 
%  
%  set(gca,'Ytick',0:1:10,'YTickLabel',{'0','1','2','3','4','5','6','7','8','9','10'})
%   axis([0 300 0 10]);
%   grid on
%   
%  title('Relative Percentage Error With Riemann Sum')
%  xlabel('# of Iterations')
%  ylabel('% Error')



end
Error = mean(error(90:100));
errorav(end+1) = Error;

xxx(end+1) = nn;
plot(xxx,errorav)
 xlabel('Number of Points in 490050 m^3')
 ylabel('Relative Percentage Error')

%set(gca,'Xtick',0:7000:35000,'XTickLabel',{'0','7000','14000','21000','28000','35000'})
% set(gca,'Ytick',0:10:100,'YTickLabel',{'0','10','20','30','40','50','60','70','80','90','100'})
% axis([0 35000 0 100]);
title('Relative Percentage Error vs Number of Points')
   grid on


end
