#  Center of mass of a wire
# Suppose  we  have a wire  l   feet long  whose  density is rho(x) 
# pounds per foot at the point x feet from the left hand end of the
# wire.  What is the total mass of the wire and where is its center of
# mass, i.e., the point cm about which the total moment of the wire is
# 0?
# Mass  Chop the wire into  n small pieces each  Delta*x[i] feet long
# and pick an arbitrary point c[i] in each piece.  An approximation to
# the mass of the ith piece of wire is rho(c[i])*Delta*x[i], so an
# approximation to the total mass   is Sum(rho(c[i])*Delta*x[i],i=1..n)
# .   This approximate mass is a Riemann sum approximating the integral 
#  Int(rho(x),x=0..l), and so the mass of the wire is defined as the
# value of this integral.
# Center of mass:    Chopping as above,  the approximate moment of the
# ith piece about the center of mass cm is 
# (c[i]-cm)*rho(c[i])*Delta*x[i] and so the total approximate moment is
# Sum((c[i]-cm)*rho(c[i])*Delta*x[i],i=1..n).  This is seen to be a
# Riemann sum approximating the integral  Int((x-cm)*rho(x),x=0..l).  
# But the center of mass is defined as the point about which the total
# moment is zero so the integral satisfies the equation
# Int((x-cm)*rho(x),x=0..l)=0.   Using properties of integrals, we can
# solve this equation for cm, to get the ratio of integrals 
# cm=Int(x*rho(x),x=0..l)/Int(rho(x),x=0..l) .   Note the top integral
# represents the total moment of the wire about its left end (x=0) and
# the bottom integral is the total mass of the wire.
> 
#   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} )   ;
# Notice that the center of mass of the solid of revolution is the same
# as the center of mass of a wire whose density at x is the area of the
# cross-section.
# 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 := proc(f, a, b)
    int(x*f(x)^2, x = a .. b)/int(f(x)^2, x = a .. b)
end

# For example, the center of  the solid obtained by rotating the region
# R under the graph of  y=cos(x) for x between 0 and pi/2 is
> cenmass(cos,0,Pi/2);

                                    2
                             1/16 Pi  - 1/4
                           4 --------------
                                   Pi

# 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 := proc(f, a, b)
local cm, solid;
    cm := plots[pointplot3d]([evalf(cenmass(f, a, b)), 0, 0],
        color = red, symbol = box, thickness = 3);
    solid := plots[tubeplot]([x, 0, 0], x = a .. b, radius = f(x),
        numpoints = 20, style = wireframe);
    plots[display]([cm, solid], scaling = constrained)
end

# Test this out.
>  drawit(cos+2,0,7);

# 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);