source("http://macosa.dima.unige.it/r.R")    # If I have not already loaded the library
---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
         
Volumes of solids of revolution
 
Rotation of a plane curve around a straight line that lies on the same plane
      [disc integration which integrates along the axis of revolution]

An example (ellipsoid):
      
BF=3; HF=2.5
f = function(x) sqrt(1-(x/2)^2)*1.5
PLANE(-2,2, -1.5, 1.5)
g = function(x) -f(x)
graph(f, -2,2, "brown"); graph1(g, -2,2, "brown")
segm(0.6, f(0.6), 0.6,0, "seagreen")
line(0.6,-f(0.6), 0.6,0, "seagreen")
H = function(x) f(x)^2*pi
integral(H, -2,2)
# 18.84956
integral(H,-2,2)/pi
# 6               The volume is 6*π

Shell integration, as if the solid was composed of many small cylindrical shells

An example: volume of a kind of "glass":
                  
BF=2.5; HF=2.5
PLANE(-2.5,2.5, 0,5)
f = function(x) x^2; graph2(f, -2,2, "brown")
F <- function(x,y) x^2+y^2-1; multrarot(F, 2,1/3, 0,4, 0, "brown")
line(0,0, -2,4, "brown"); line(0,0, 2,4, "brown")
# cylindrical shells of radius x, height 2x-x^2, thickness dx
H = function(x) 2*pi*x*(2*x-x^2)
integral(H, 0,2)
# 8.37758
fraction(integral(H, 0,2)/pi)
# 8/3               The volume is 8*π/3

The volume of a barrel (with parabolic section)
 

D = 7.9; d = 6.65; h = 9.5   # dm
# y=a*x^2+c fromx=-h/2, y=-d/2 to x=h/2,y=d/2 passing through x=0,y=D/2.
# x=0, y=D/2 → c=D/2;  x=h/2, y=d/2 → a*(h/2)^2+D/2 = d/2 → a=(d/2-D/2)/(h/2)^2
# y = (d/2-D/2)/(h/2)^2*x^2 + D/2
f = function(x) (d/2-D/2)/(h/2)^2*x^2+D/2
g = function(x) f(x)^2*pi; integral(g, -h/2,h/2)
# 418.8702    (419 liters)
# It can be find that:
V = function(h,D,d) pi*h*(3*d^2 + 4*d*D + 8*D^2)/60
# verification:
V(9.5, 7.9, 6.65)
# 418.8702       OK
#
# The figure:
BF=4; HF=4; BOXW(0,0,0,0) # I use BOXW to reduce white space
x = c(-7.9/2,7.9/2); y = c(-7.9/2,7.9/2) ; w = c(0,9.5) # the box
z = array(c(w[1],w[1],w[1],w[1]), dim=c(2,2))
th = 120; ph = 40
F = persp(x,y,z, theta=th,phi=ph, scale=FALSE, xlim=x,ylim=y,zlim=w,d=3,border="brown")
# border="white", box=FALSE  and  co="white", if I don't want the box
co="red"
axes = function(F) { 
  lines(trans3d(c(0,0),c(0,0),c(0,w[2]),pmat=F),col=co)
  lines(trans3d(c(0,x[2]),c(0,0),c(0,0),pmat=F),col=co)
  lines(trans3d(c(0,0),c(0,y[2]),c(0,0),pmat=F),col=co)
  lines(trans3d(c(0,0),c(0,0),c(w[1],0),pmat=F),col=co,lty=3)
  lines(trans3d(c(x[1],0),c(0,0),c(0,0),pmat=F),col=co,lty=3)
  lines(trans3d(c(0,0),c(y[1],0),c(0,0),pmat=F),col=co,lty=3) }
axes(F)
t=seq(0,9.5,0.1)
lines(trans3d(0,f(t-9.5/2),t,pmat=F),col="black")
for(a in seq(0,360,360/20)) lines(trans3d(cos(a*pi/180)*f(t-9.5/2),sin(a*pi/180)*f(t-9.5/2),t,pmat=F),col="black")
a = (0:50)*2*pi/50
lines(trans3d(cos(a)*d/2,sin(a)*d/2,0,pmat=F),col="blue")
lines(trans3d(cos(a)*d/2,sin(a)*d/2,9.5,pmat=F),col="blue")
lines(trans3d(cos(a)*D/2,sin(a)*D/2,9.5/2,pmat=F),col="blue")


Other examples of use