# 
#  Center of mass of a solid of revolution
# 
# If  y *`=` *f(x) >= 0 for  a <= x *`<=` b, then let S be the solid of
# revolution obtained by rotating
# the region under the graph of f around the x axis. We know how to
# express the volume of S as an integral:
# 
# Just integrate from a to b the  crossectional area Pi*f(x)^2 of the
# solid S. 
> Volume = Int(Pi*f(x)^2,x=a..b);
# Now how to find the center of mass of the solid, assuming it's made of
# a homogeneous material?
# 
# Well, it's  clear that the center of mass will be somewhere along the
# x-axis between a and b.   
# 
# Let CM be the center of mass.  Partition  [a,b] into n subintervals
# [x[i],x[i+1] ] and using planes 
# perpendicular to [a,b] approximate the solid S with the n disks where
# the ith one has volume   Delta V[i] = Pi*f(x[i])^2*Delta x[i]
# 
# Now the signed moment of the ith disk about the point CM is   M[i] =
# (x-CM)* Delta V[i]  and the sum of these moments will be approximately
# 0, since CM is the center of mass.  If we let Delta x[i] go to zero 
# this  approximate equation becomes an equation for the center of mass:
> CMequation :=  int((x-CM)*Pi*f(x)^2,x=a..b) =0;
# Useing properties of integrals, we can solve this equation for CM.
> sol := solve(CMequation,{CM} )   ;
# We can define a word cenmass  which takes a function  f, an interval
# [a,b], and locates the center of the solid of revolution.
> cenmass := proc(f,a,b)  
>  int(x*f(x)^2,x=a..b)/int(f(x)^2,x=a..b) end;
> cenmass(cos,0,1);
# Now we can define a word to draw the solid and locate the center of
# mass.
> drawit := proc(f,a,b) 
>  local cm, solid;
>   cm :=
> plots[pointplot3d]([evalf(cenmass(f,a,b)),0,0],color=red,symbol=box,th
> ickness=3):
>   solid :=
> plots[tubeplot]([x,0,0],x=a..b,radius=f(x),numpoints=20,style=wirefram
> e); 
> plots[display]([cm,solid],scaling=constrained); end; 
>  drawit(cos+2,0,7);
>  
> plots[pointplot3d]([evalf(cenmass(cos,0,1)),0,0],color=red,symbol=box,
> thickness=3);
> 
# We can animate the motion of the center of mass as the solid changes.
# 
> plots[display]([seq(drawit(2+cos,0,t),t=1..10)],insequence=true);