

f=@(x,y) 2*(cos(x)+x/5)+sin(y+1);


% första uppdelningen, hörnpunkt vald
Q=-1;
for N=[10 10];
Q=Q+1; %f(x2,y2) används
% sidoytor och ytans kontur 
x=linspace(-1,1,90);y=linspace(0,2,71);z=linspace(0,4,71);
[X,Z]=meshgrid(x,z);

% sidoyta y=0
Y=zeros(size(X));
Z(find(Z>f(X,Y)))=NaN;
surf(X,Y,Z,ones(size(X)),'FaceAlpha',0.1)
shading interp
axis([-1.5 1.5 -0.5 2.5 0 4])
hold on

% sidoyta y=2
Y=2*ones(size(X));
Z(find(Z>f(X,Y)))=NaN;
surf(X,Y,Z,ones(size(X)),'FaceAlpha',0.1)
shading interp

% sidoyta x=1
[Y,Z]=meshgrid(y,z);
X=ones(size(Y));
Z(find(Z>f(X,Y)))=NaN;
surf(X,Y,Z,ones(size(X)),'FaceAlpha',0.1)
shading interp

% sidoyta x=-1
[Y,Z]=meshgrid(y,z);
X=-ones(size(Y));
Z(find(Z>f(X,Y)))=NaN;
surf(X,Y,Z,ones(size(X)),'FaceAlpha',0.1)
shading interp


x=linspace(-1,1,21);y=linspace(0,2,21);
[X,Y]=meshgrid(x,y);
Z=f(X,Y);

% konturerna till kroppen
plot3([-1 1 1 -1 -1 -1],[0 0 2 2 0 0],[0 0 0 0 0 f(-1,0)],'black')
plot3([1 1],[2 2],[0 f(1,2)],'black')
plot3([1 1],[0 0],[0 f(1,0)],'black')
plot3([-1 -1],[2 2],[0 f(-1,2)],'black')

plot3(ones(size(y)),y,f(ones(size(y)),y),'black')
plot3(-ones(size(y)),y,f(-ones(size(y)),y),'black')
plot3(x,0*ones(size(x)),f(x,0*ones(size(x))),'black')
plot3(x,2*ones(size(x)),f(x,2*ones(size(x))),'black')


% rätblocken

x=linspace(0,0.1,90);y=linspace(0.5,0.6,71);z=linspace(0,4,71);
int=0;h=2/N;
te=text(-1,2,3,['summa = ' ]);
for k=0:(N-1)
    for m=0:(N-1)
       
       x1=k*h-1;x2=(k+1)*h-1;y1=m*h;y2=(m+1)*h;
      
       x0=x2-Q*h*rand;y0=y2-Q*h*rand;
       z1=f(x0,y0);
       int=int+z1*h^2;
       % ruta i xy-planet
       plot3([x1 x1 x2 x2 x1],[y1 y2 y2  y1 y1], [0 0 0 0 0])
       
       
       if (k+m)==0
           %linje från slumppunkt/vald punkt till ytan
       plot3(x0,y0,0,'r*')
       p1 = plot3(x0,y0,z1,'r*')
       p2 = plot3([x0 x0],[y0 y0],[0 z1],'r-');
       
       pause
       delete(p1,p2)
       % rätblocket
        %linjer upp till f-nivå i hörnen
       p1=plot3([x2 x2],[y2 y2],[0 z1],'-.');
       p2=plot3([x1 x1],[y1 y1],[0 z1],'-.');
       p3=plot3([x1 x1],[y2 y2],[0 z1],'-.');
       p4=plot3([x2 x2],[y1 y1],[0 z1],'-.');
        % ram runt f-nivån 
       plot3([x1 x1 x2 x2 x1],[y1 y2 y2  y1 y1],[z1 z1 z1 z1 z1 ])
       
        % toppytan
       X=[x1 x1;x2 x2];
       Y=[y1 y2; y1 y2];
       Z=[z1 z1;z1 z1];
       surf(X,Y,Z,'FaceAlpha',1)
       shading interp
       delete(te)
       te=text(-1,2,3,['summa = ' num2str(int)]);
           pause
       else
               %slumppunkt
       plot3(x0,y0,0,'r*')
       % p1 = plot3(x0,y0,z1,'r*');
       % p2 = plot3([x0 x0],[y0 y0],[0 z1],'r-');
       % pause
       % delete(p1)
       % rätblocket
        %linjer upp till f-nivå i hörnen
       p1=plot3([x2 x2],[y2 y2],[0 z1],'-.');
       p2=plot3([x1 x1],[y1 y1],[0 z1],'-.');
       p3=plot3([x1 x1],[y2 y2],[0 z1],'-.');
       p4=plot3([x2 x2],[y1 y1],[0 z1],'-.');
        % ram runt f-nivån 
       plot3([x1 x1 x2 x2 x1],[y1 y2 y2  y1 y1],[z1 z1 z1 z1 z1 ])
       
        % toppytan
       X=[x1 x1;x2 x2];
       Y=[y1 y2; y1 y2];
       Z=[z1 z1;z1 z1];
       surf(X,Y,Z,'FaceAlpha',1)
       shading interp
       delete(te)
       te=text(-1,2,3,['summa = ' num2str(int)]);
            pause(.003)
       end
       delete(p1,p2,p3,p4)
       
    end
end

rotate3d on
%grid off
hold off
pause
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % fortsatta uppdelningar, punkter slumpas
% Q=1; 
% for N=[10 10 ];
% x=linspace(-1,1,90);y=linspace(0,2,71);z=linspace(0,4,71);
% [X,Z]=meshgrid(x,z);
% 
% % sidoyta y=0
% Y=zeros(size(X));
% Z(find(Z>f(X,Y)))=NaN;
% surf(X,Y,Z,ones(size(X)),'FaceAlpha',0.1)
% shading interp
% axis([-1.5 1.5 -0.5 2.5 0 4])
% hold on
% 
% % sidoyta y=2
% Y=2*ones(size(X));
% Z(find(Z>f(X,Y)))=NaN;
% surf(X,Y,Z,ones(size(X)),'FaceAlpha',0.1)
% shading interp
% 
% % sidoyta x=1
% [Y,Z]=meshgrid(y,z);
% X=ones(size(Y));
% Z(find(Z>f(X,Y)))=NaN;
% surf(X,Y,Z,ones(size(X)),'FaceAlpha',0.1)
% shading interp
% 
% % sidoyta x=-1
% [Y,Z]=meshgrid(y,z);
% X=-ones(size(Y));
% Z(find(Z>f(X,Y)))=NaN;
% surf(X,Y,Z,ones(size(X)),'FaceAlpha',0.1)
% shading interp
% 
% 
% x=linspace(-1,1,21);y=linspace(0,2,21);
% [X,Y]=meshgrid(x,y);
% Z=f(X,Y);
% 
% % konturerna till kroppen
% plot3([-1 1 1 -1 -1 -1],[0 0 2 2 0 0],[0 0 0 0 0 f(-1,0)],'black')
% plot3([1 1],[2 2],[0 f(1,2)],'black')
% plot3([1 1],[0 0],[0 f(1,0)],'black')
% plot3([-1 -1],[2 2],[0 f(-1,2)],'black')
% 
% plot3(ones(size(y)),y,f(ones(size(y)),y),'black')
% plot3(-ones(size(y)),y,f(-ones(size(y)),y),'black')
% plot3(x,0*ones(size(x)),f(x,0*ones(size(x))),'black')
% plot3(x,2*ones(size(x)),f(x,2*ones(size(x))),'black')
% 
% x=linspace(0,0.1,90);y=linspace(0.5,0.6,71);z=linspace(0,4,71);
% int=0;h=2/N;
% for k=0:(N-1)
%     for m=0:(N-1)
%        
%        x1=k*h-1;x2=(k+1)*h-1;y1=m*h;y2=(m+1)*h;
%        % z1=f(x2,y2);
%        x0=x2-Q*h*rand;y0=y2-Q*h*rand;
%        z1=f(x0,y0);
%        int=int+z1;
%        % ruta i xy-planet
%        plot3([x1 x1 x2 x2 x1],[y1 y2 y2  y1 y1], [0 0 0 0 0])
%        if N<=10
%            if (k+m)==0
%            %linje från slumppunkt/vald punkt till ytan
%        plot3(x0,y0,0,'r*')
%        p1 = plot3(x0,y0,z1,'r*');
%        p2 = plot3([x0 x0],[y0 y0],[0 z1],'r-');
%        pause
%        delete(p1,p2)
%        % rätblocket
%         %linjer upp till f-nivå i hörnen
%        p1=plot3([x2 x2],[y2 y2],[0 z1],'-.');
%        p2=plot3([x1 x1],[y1 y1],[0 z1],'-.');
%        p3=plot3([x1 x1],[y2 y2],[0 z1],'-.');
%        p4=plot3([x2 x2],[y1 y1],[0 z1],'-.');
%         % ram runt f-nivån 
%        plot3([x1 x1 x2 x2 x1],[y1 y2 y2  y1 y1],[z1 z1 z1 z1 z1 ])
%        
%         % toppytan
%        X=[x1 x1;x2 x2];
%        Y=[y1 y2; y1 y2];
%        Z=[z1 z1;z1 z1];
%        surf(X,Y,Z,'FaceAlpha',1)
%        shading interp
%        
%            pause
%        else
%         plot3(x0,y0,0,'r*')
%         %linjer upp till f-nivå i hörnen
%       % p1=plot3([x2 x2],[y2 y2],[0 z1],'-.');
%       % p2=plot3([x1 x1],[y1 y1],[0 z1],'-.');
%       % p3=plot3([x1 x1],[y2 y2],[0 z1],'-.');
%       % p4=plot3([x2 x2],[y1 y1],[0 z1],'-.');
%         % ram runt f-nivån 
%        plot3([x1 x1 x2 x2 x1],[y1 y2 y2  y1 y1],[z1 z1 z1 z1 z1 ])
%        
%         % toppytan
%        X=[x1 x1;x2 x2];
%        Y=[y1 y2; y1 y2];
%        Z=[z1 z1;z1 z1];
%        surf(X,Y,Z,'FaceAlpha',1)
%        shading interp
%        
%            end
%        else
%        % ram runt f-nivån 
%        plot3([x1 x1 x2 x2 x1],[y1 y2 y2  y1 y1],[z1 z1 z1 z1 z1 ])
%        %linjer upp till f-nivå i hörnen
%        %plot3([x2 x2],[y2 y2],[0 z1],'-.')
%        %plot3([x1 x1],[y1 y1],[0 z1],'-.')
%        %plot3([x1 x1],[y2 y2],[0 z1],'-.')
%        %plot3([x2 x2],[y1 y1],[0 z1],'-.')
%        
%        % plot3(x0,y0,z1,'r*')
%        % plot3([x0 x0],[y0 y0],[0 z1],'-')
%        
%        X=[x1 x1;x2 x2];
%        Y=[y1 y2; y1 y2];
%        Z=[z1 z1;z1 z1];
%        surf(X,Y,Z,'FaceAlpha',1)
%        shading interp
%        
%        % pause(.003)
%        end
%     end
% end
% int=int*h^2;
% text(-1,2,3,['summa = ' num2str(int)])
% rotate3d on
% %grid off
% hold off
% 
% pause
% end