> #sfb sep 25, 2001
> with(plots);#Just need the display function to combine plots
> circle:=t=0..2*Pi;#To help us describe the rotation, used in plots
> #
> #The GIVENS (x0 is any point between a and b but it needs to be given a value at this point)
> #x^2,1..3;1/x,1..10,1+cos,0..4pi
> f:=x^2;a:=1;b:=3;x0:=(a+2*b)/3;
> f:=1/x;a:=10;b:=100;x0:=(a+2*b)/3;
> f:=1+cos(x);a:=0;b:=4*Pi;x0:=(a+2*b)/3;
> fa:=subs(x=a,f);fb:=subs(x=b,f);fx0:=subs(x=x0,f);
> xrange:=x=a..b;xgraph:=x=0..b;#Used in plots
> ygraph:=y=0..max(fa,fb);#assumes f is monotone
> F:=plot(f,xgraph,ygraph,title="y=f(x)"):F;#initial plot 
> W:=plot3d([x,f,0],xrange,ygraph,axes=BOXED,thickness=3,color=plum):
> W:#interesting once
> #
> #Our first look at rotation, first about the X axis and then the Y axis
> #Look at orientation "-90,0" is the standard view & spin the zero
> X:=plot3d([x,f*cos(t),f*sin(t)],xrange,circle): display(X,W,title="about X axes");
> Y:=plot3d([cos(t)*x, f, x*sin(t)],xrange,circle): display(Y,W,title="about Y axes");
> display(X,Y,W,title="have it both ways");
> #
> # Show the area so we can show the volume
> Ff:=plot(f,xrange,filled=true,title="y=f(x)"):Ff;
> Wf:=plot3d([x,s*f,0],xrange,s=0..1,axes=BOXED,thickness=3,color=plum):
> Wf:#interesting once
> Xa:=plot3d([a,fa*cos(t)*s,fa*sin(t)*s],s=0..1,circle): Xb:=plot3d([b,fb*cos(t)*s,fb*sin(t)*s],s=0..1,circle):
> display(X,Xa,Xb,Wf,title="about X axes");
> Yb:=plot3d([cos(t)*b, s*fb, b*sin(t)],s=0..1,circle):  Ya:=plot3d([cos(t)*a, s*fa, a*sin(t)],s=0..1,circle): display(Y,Ya,Yb,Wf,title="about Y axes");
> #
> # Computing the volume--derive the integral
> Fdx:=plot([x0,s*fx0,s=0..1],xrange):display(F,Fdx);
> dxx:=plot3d([x0,s*fx0*cos(t),s*fx0*sin(t)],circle,s=0..1,thickness=3,color=aquamarine):
> dxx:#interesting once
> display(X,dxx,W,title="a little dx about the X axes");
> dxy:=plot3d([cos(t)*x0,s*fx0,x0*sin(t)],s=0..1,circle,thickness=3,color=aquamarine): 
> dxy:#interesting once
> display(Y,Yb,Ya,dxy,W,title="a little dx about the Y axes");
> #
> #Computing the volume
> VolX:=int(Pi*f^2,x=a..b);evalf(VolX);
> VolY:=int(2*Pi*x*f,x=a..b);evalf(VolY);
> #
> #
> #Suppose the area inside x^2+(y-R)^2=r^2 is rotated about the x axes
> funs:=solve(x^2+(y-R)^2=r^2,y);
> 
> f:=funs[1];g:=funs[2];
> r:=1;R:=2;a:=-r;b:=r;x0:=(a+2*b)/3;xrange:=x=a..b;fx0:=subs(x=x0,f); gx0:=subs(x=x0,g);
> Ff:=plot3d([x,f*cos(t),f*sin(t)],xrange,circle,style=WIREFRAME):Ff;
> Fg:=plot3d([x,g*cos(t),g*sin(t)],xrange,circle,style=WIREFRAME):Fg;
> Fx0:=plot3d([x0,s*cos(t),s*sin(t)],circle,s=gx0..fx0,color=plum):Fx0;
> gx0..fx0;
> display(Ff,Fg,Fx0,title="donut",scaling=constrained);
> volume:=int(Pi*(f^2-g^2),x=a..b);
> expand(f^2-g^2);
>