 clf
f=@(x,y) 2*(cos(x)+x/5)+sin(y+1);
g1=@(x) sin(x+1)/2;
g2=@(x) sin(2*x+1)/2+1.5;
x=linspace(-1,1,60);y=linspace(0,2,61);z=linspace(0,4,71);
[X,Z]=meshgrid(x,z);

Y=g1(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

Y=g2(X);
Z(find(Z>f(X,Y)))=NaN;
surf(X,Y,Z,ones(size(X)),'FaceAlpha',0.1)
shading interp

plot3(x,g1(x),f(x,g1(x)),'black')
plot3(x,g2(x),f(x,g2(x)),'black')

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

[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

[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

f=@(x,y) (2*(cos(x)+x/5)+sin(y+1)).*(y>sin(x+1)/2).*(y<sin(2*x+1)/2+1.5);
x=linspace(-1,1,120);y=linspace(0,2,121);
[X,Y]=meshgrid(x,y);
Z=f(X,Y);
Z(find(Z==0))=NaN;
ytan=surf(X,Y,Z);
shading interp
rotate3D on

plot3(x,g1(x),zeros(size(x)),'black')
plot3(x,g2(x),zeros(size(x)),'black')
plot3(x,g1(x),f(x,g1(x)),'black')
plot3(x,g2(x),f(x,g2(x)),'black')
pause
delete(ytan)
xs=linspace(-1,1,20);
for n=1:19
    xa=xs(n);
    xb=xs(n+1);
    [Y,Z]=meshgrid(y,z);
    X=xa*ones(size(Y));
    Xb=xb*ones(size(Y));
    Z(find(Z>f(Xb,Y)))=NaN;
    sur1=surf(X,Y,Z,3*ones(size(X)),'FaceAlpha',0.1);
    shading interp

    X=xb*ones(size(Y));
    Z(find(Z>f(X,Y)))=NaN;
    sur2=surf(X,Y,Z,3*ones(size(X)),'FaceAlpha',0.1);
    shading interp
    alpha(0.2)
    
    xr=linspace(xa,xb,30);y=linspace(g1(xa),g2(xa),31);
    [X,Y]=meshgrid(xr,y);
    Z=f(xb*ones(size(X)),Y);
    Z(find(Z==0))=NaN;
    surf(X,Y,Z,1*ones(size(X)),'FaceAlpha',0.4)
    shading interp
    
    ft=@(x,y) (2*(cos(x)+x/5)+sin(y+1));
    kur1=plot3([xa xa xa xa],[g1(xb) g1(xb) g2(xb) g2(xb)],[ft(xb,g1(xb)) 0 0 ft(xb,g2(xb))]);
    y=linspace(g1(xb),g2(xb),71);
    plot3(xa*ones(size(y)),y,f(xb*ones(size(y)),y))
    kur2=plot3([xb xb xb xb],[g1(xb) g1(xb) g2(xb) g2(xb)],[ft(xb,g1(xb)) 0 0 ft(xb,g2(xb))]);
    y=linspace(g1(xb),g2(xb),71);
    plot3(xb*ones(size(y)),y,f(xb*ones(size(y)),y))
    pause(0.1)
%     plot3([0 0 0 0],[g1(0) g1(0) g2(0) g2(0)],[f(0.1,0) 0 0 f(0.1,2)])
%     y=linspace(g1(0.1),g2(0.1),71);
%     plot3(0*ones(size(y)),y,f(0.1*ones(size(y)),y))
%     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')
%     pause
%     y=linspace(g1(-1),g2(-1),71);
%     plot3(-ones(size(y)),y,f(-ones(size(y)),y))
%     y=linspace(g1(1),g2(1),71);
%     plot3(ones(size(y)),y,f(ones(size(y)),y))
%     plot3([-1 -1 -1 -1],[g1(-1) g1(-1) g2(-1) g2(-1)],[f(-1,g1(-1)) 0 0 f(-1,g2(-1))],'black')
%     plot3([1 1 1 1],[g1(1) g1(1) g2(1) g2(1)],[f(1,g1(1)) 0 0 f(1,g2(1))],'black')
    
     delete(sur1);
     delete(sur2);
%     delete(kur1);
%     delete(kur2);
end

% for k=1:9
%     x=k*0.1-1;
%     plot3([x x],[g1(x) g2(x)],[0 0],'cyan')
% end
% for k=12:19
%     x=k*0.1-1;
%     plot3([x x],[g1(x) g2(x)],[0 0],'cyan')
% end

hold off