%SHELLMETHOD                       <<last updated  8/13/01>>
%
%This is an instructional demo for calculus.
%It illustrates graphically the shell method for computing volumes
%of solids of revolution for regions of area with the following constraints:
%       >The region is in the first quadrant.
%       >The region is between a curve y = f(x) & 
%            the x-axis over an interval [a,b].
%       >The rotation is about the y-axis.
%       >THREE shells are generated as specified by the user.
%
%Notes: >The region can have a hole in the center.
%       >It is not designed for regions between two curves y = f(x) & y = g(x).
%       >It works best for monotonic curves over [a,b].
%       >Some shells can 'hide' other shells.
%
%DETAILS: There are SIX built-in demos. These are displayed in a menu.
%         There is a sequence of pictures generated in subplots.
%              1. graph of f(x) over [a,b] with region filled in color
%              2. graph of f(x) over [a,b] with three vertical strips shown
%                 that will generate the shells
%              3. three graphs of the individual shells generated by 
%                 rotating the 3 strips
%              4. a final picture showing all 3 shells together
%                 (in some cases this last picture will be rotated 
%                  to clearly show the 3 shells; this feature is demo
%                  dependent & must be planned for by the user)
%
%EXTENSIONS: Functions can be added to the menu with minor modifications.
%            Below is a sample of the code that must be supplied 
%            with annotations.
%     if ch==4,f='cos(pi/2*x)';                    <=function formula
%         strip1=[1/2 1/2 5/8 5/8 1/2];x1=strip1;  <=x-coordinates of strips
%         strip2=[5/8 5/8 3/4 3/4 5/8];x2=strip2;
%         strip3=[3/4 3/4 1 1 3/4];x3=strip3;
%         a=.5;b=1;grdom=[0 1];                    <=interval [a,b] &
%         %setup function info                       grdom, the x-domain used
%         vf=vectorize(f);                           in the sketch
%         x=a:.05:b;xval=x;
%         fval=eval(vf);                           <=computing f over [a,b]
%         height=max(fval);                        <=max function value
%         %region boundary
%         vertpiecex=[xval(1) xval(1)];vertpiecey=[0 fval(1)]; <=coordinates
%                             to draw a vertical line to close the region
%         fboundaryx=[xval xval(1)];fboundaryy=[fval 0];  <=coordinates used
%                             in the fill command to shade the region
%         goodch='Y';    <=signal that we have a good input choice code
%      end
%          
%
%AUXILIARY ROUTINES  addaxes by D. Hill
%
%     By: David R. Hill, MATH DEPT, Temple University,
%         Philadelphia, Pa. 19122   Email: dhill001@temple.edu

%strings
s0=' ';
contmess='Press ENTER to see next figure.';
endmess='Routine is over; press ENTER.';
header='            Solids of Revolution by the SHELL METHOD.';
menu=['            OPTIONS                               ';
      '                                                  ';
      '1. Demo with f = x^2     over [0,1].              ';
      '                                                  ';
      '2. Demo with f = 1 - x^2 over [0,1].              ';
      '                                                  ';
      '3. Demo with f = 2 - x over [0,1]                 ';
      '                                                  ';
      '4. Demo with f = cos(pi/2*x) over [1/2,1]         ';
      '                                                  ';
      '5. Demo with f = sqrt(x) over [0,1].              ';
      '                                                  ';
      '6. Demo with f = -2x^2 + 8x - 6 over [1,3].       ';
      '                                                  ';
      '0. QUIT!                                          ';
      '                                                  '];
choice='Enter the number of your choice:  ==> ';
again='Improper choice; TRY AGAIN.';
enter='Press ENTER to continue.';
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%COLORMAP DEFINITIONS for solid colors
myred=zeros(64,3);
myred(:,1)=1*ones(64,1);
mygreen=zeros(64,3);
mygreen(:,2)=1*ones(64,1);
myblue=zeros(64,3);
myblue(:,3)=1*ones(64,1);
%++++++++++++++++++++++++++++
%
%START  <><><><><><><><><><><><><><><><><><><><><><><><><><><>
%
%INPUT MENU & CHOICES
goodch='N';
while goodch=='N'
      clc
      disp(header),disp(s0),disp(s0)
      disp(menu)
      ch=input(choice);
      if ch==0,return,end
      if ch==1,f='x^2';
         strip1=[1/4 1/4 1/2 1/2 1/4];x1=strip1;
         strip2=[1/2 1/2 3/4 3/4 1/2];x2=strip2;
         strip3=[3/4 3/4 1 1 3/4];x3=strip3;
         a=0;b=1;grdom=[0 1];
         %setup function info
         vf=vectorize(f);
         x=a:.05:b;xval=x;
         fval=eval(vf);
         height=max(fval);
         %region boundary
         vertpiecex=[1 1];vertpiecey=[0 1];
         fboundaryx=[xval 1 xval(1)];fboundaryy=[fval 0 0];
         goodch='Y';
      end
      if ch==2,f='1-x^2';
         strip1=[1/4 1/4 1/2 1/2 1/4];x1=strip1;
         strip2=[1/2 1/2 3/4 3/4 1/2];x2=strip2;
         strip3=[3/4 3/4 1 1 3/4];x3=strip3;
         a=0;b=1;grdom=[0 1];
         %setup function info
         vf=vectorize(f);
         x=a:.05:b;xval=x;
         fval=eval(vf);
         height=max(fval);
         %region boundary
         vertpiecex=[xval(1) xval(1)];vertpiecey=[0 fval(1)];
         fboundaryx=[xval xval(1)];fboundaryy=[fval 0];
         goodch='Y';
      end
      if ch==3,f='2-x';
         strip1=[1/4 1/4 1/2 1/2 1/4];x1=strip1;
         strip2=[1/2 1/2 3/4 3/4 1/2];x2=strip2;
         strip3=[3/4 3/4 1 1 3/4];x3=strip3;
         a=0;b=1;grdom=[0 1];
         %setup function info
         vf=vectorize(f);
         x=a:.05:b;xval=x;
         fval=eval(vf);
         height=max(fval);
         %region boundary
         n=length(xval);
         vertpiecex=[xval(n) xval(n)];vertpiecey=[0 fval(n)];
         fboundaryx=[xval xval(n) xval(1)];fboundaryy=[fval 0 0];
         goodch='Y';
      end
      if ch==4,f='cos(pi/2*x)';
         strip1=[1/2 1/2 5/8 5/8 1/2];x1=strip1;
         strip2=[5/8 5/8 3/4 3/4 5/8];x2=strip2;
         strip3=[3/4 3/4 1 1 3/4];x3=strip3;
         a=.5;b=1;grdom=[0 1];
         %setup function info
         vf=vectorize(f);
         x=a:.05:b;xval=x;
         fval=eval(vf);
         height=max(fval);
         %region boundary
         vertpiecex=[xval(1) xval(1)];vertpiecey=[0 fval(1)];
         fboundaryx=[xval xval(1)];fboundaryy=[fval 0];
         goodch='Y';
      end
      if ch==5,f='sqrt(x)';
         strip1=[1/4 1/4 1/2 1/2 1/4];x1=strip1;
         strip2=[1/2 1/2 3/4 3/4 1/2];x2=strip2;
         strip3=[3/4 3/4 1 1 3/4];x3=strip3;
         a=0;b=1;grdom=[0 1];
         %setup function info
         vf=vectorize(f);
         x=a:.05:b;xval=x;
         fval=eval(vf);
         height=max(fval);
         %region boundary
         n=length(xval);
         vertpiecex=[xval(n) xval(n)];vertpiecey=[0 fval(n)];
         fboundaryx=[xval xval(n) xval(1)];fboundaryy=[fval 0 0];
         goodch='Y';
      end
      if ch==6,f='-2*x^2+8*x-6';
         strip1=[1.5 1.5 2 2 1.5];x1=strip1;
         strip2=[2 2 2.5 2.5 2];x2=strip2;
         strip3=[2.5 2.5 3 3 2.5];x3=strip3;
         a=1;b=3;grdom=[0 3];
         %setup function info
         vf=vectorize(f);
         x=a:.05:b;xval=x;
         fval=eval(vf);
         height=max(fval);
         %region boundary
         vertpiecex=[0];vertpiecey=[0];
         fboundaryx=[xval];fboundaryy=[fval];
         goodch='Y';
      end
      if goodch=='N'  %BAD input number
         disp(s0),disp(again),disp(enter),pause
      end
end %WHILE
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%
%SET UP STRIP y-data
x=x1(1);
y1=[0 eval(vf) eval(vf) 0 0];
x=x2(1);
y2=[0 eval(vf) eval(vf) 0 0];
x=x3(1);
y3=[0 eval(vf) eval(vf) 0 0];
%+++++++++++++++++++++++++++++++++
%setting up circle domain values
t=0:.5:2*pi+.1;
%++++++++++++++++++++++++++++++++
%GENERATING the graphics figure
figh=figure('color','w','units','normal','position',[0 0 1 1]);
%VANITY
vanity=uicontrol(figh, 'units', 'normalized', 'style', 'text',...
   'string', 'by D.R.Hill', 'position', [.80 .02 .29 .04],...
   'fontsize', 10,'fontweight','bold','fontangle','oblique','backgroundcolor','w');
%SET UP MESSAGE AREA
messsub=subplot('position',[.3 .05 .5 .05]);axis('off')
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%generate graph of region in 2D
subplot(3,3,1)
plot(xval,fval,'-k','linewidth',2)
ax2=[-grdom(2),grdom(2),0,height]; %assumed grdom(2)>0
axis(ax2)
addaxes
hold on
%plot vertical piece of region
plot(vertpiecex,vertpiecey,'-k','linewidth',2)
%FILL the region
fill(fboundaryx,fboundaryy,'yellow')
title('The region to rotate')
xlabel('about the y-axis.')
%++++
%DISPLAY formula
subplot(334)
axis('off')
text(0,.5,['f(x) = ' f],'fontweight','bold','color','red','fontsize',14)
%++++
subplot(3,3,2)
plot(xval,fval,'-k','linewidth',2)
ax2=[-grdom(2),grdom(2),0,height]; %assumed grdom(2)>0
axis(ax2)
addaxes
hold on
plot(x1,y1,'-m','linewidth',2)
plot(x2,y2,'-m','linewidth',2)
plot(x3,y3,'-m','linewidth',2)
title('Strips to rotate.')
%+++++
subplot(messsub)
mess1=text(0,0,contmess,'fontweight','bold','fontsize',14,'color','magenta');
pause
%+++++
subplot(332)
fill(x1,y1,'red')
plot(xval,fval,'-k','linewidth',2)
subplot(333)
surfl(cos(t)'*x1,sin(t)'*x1,ones(size(t))'*y1);
b=grdom(2);
ax=[-b b -b b 0 height];%assumes b positive & height positive
axis(1.5*ax);
colormap(myred);
title('Shell from strip #1.','fontweight','bold','color','red')
pause
%+++++
subplot(332)
fill(x2,y2,'green')
plot(xval,fval,'-k','linewidth',2)
subplot(336)
surfl(cos(t)'*x2,sin(t)'*x2,ones(size(t))'*y2);
axis(1.5*ax);
colormap(mygreen)
title('Shell from strip #2.','fontweight','bold','color','green')
pause
%+++++
subplot(332)
fill(x3,y3,'blue')
plot(xval,fval,'-k','linewidth',2)
subplot(339)
surfl(cos(t)'*x3,sin(t)'*x3,ones(size(t))'*y3);
axis(1.5*ax);
colormap(myblue)
title('Shell from strip #3.','fontweight','bold','color','blue')
pause
%+++++
subplot(337)
surfl(cos(t)'*x1,sin(t)'*x1,ones(size(t))'*y1);
axis(1.5*ax);
colormap(prism)
hold on
pause(1)
surfl(cos(t)'*x2,sin(t)'*x2,ones(size(t))'*y2);
pause(1)
surfl(cos(t)'*x3,sin(t)'*x3,ones(size(t))'*y3);
title('The three shells.','fontweight','bold')
if ch==1|ch==5|ch==6
   pause(2)
   view(-41,64);
   title('Rotated to see the shells.','fontweight','bold')
end
delete(mess1);
pause(2)
subplot(messsub)
mess1=text(0,0,endmess,'fontweight','bold','fontsize',14,'color','red');