# Volumes by slicing


{{< include ../_common_code.qmd >}}

This section uses these add-on packages:


```{julia}
using CalculusWithJulia
using Plots
plotly()
using QuadGK
using Unitful, UnitfulUS
using Roots
using SymPy
```

```{julia}
#| echo: false
#| results: "hidden"
import LinearAlgebra: norm, cross
using SplitApplyCombine
nothing
```


```{julia}
#| echo: false
# commands used for plotting from https://github.com/SigurdAngenent/WisconsinCalculus/blob/master/figures/221/09surf_of_rotation2.py


#linear projection of R^3 onto R^2
function _proj(X, v)
    # a is ⟂ to v and b is v × a
	vx, vy, vz = v
	a = [-vy, vx, 0]
	b = cross([vx,vy,vz], a)
	a′, b′ = a/norm(a), b/norm(b)
	return 	(a′ ⋅ X, b′ ⋅ X)
end

# project a curve in R3 onto R2
pline(viewp, ps...) = [_proj(p, viewp) for p in ps]

# determinant of Jacobian; area multiplier
# det(J); used to identify folds
function jac(X, u, v)
	return det(ForwardDiff.jacobian(xs -> collect(X(xs...)), [u,v]))
end

function _fold(F, t, θmin, θmax)
	λ = θ -> jac(F, t, θ) # F is projected surface, psurf
	iszero(λ(θmin)) && return θmin
	iszero(λ(θmax)) && return θmax

	return solve(ZeroProblem(λ, (θmin, θmax)))
end

nothing
```

---


```{julia}
#| hold: true
#| echo: false
imgfile = "figures/michelin-man.jpg"
caption = """

Hey Michelin Man, how much does that costume weigh?

"""
# ImageFile(:integrals, imgfile, caption)
nothing
```

![Hey Michelin Man, how much does that costume weigh?](./figures/michelin-man.jpg)

An ad for a summer job says work as the Michelin Man! Sounds promising, but how much will that costume weigh? A very hot summer may make walking around in a heavy costume quite uncomfortable.


A back-of-the envelope calculation would start by


  * Mentally separating out each "tire" and lining them up one by one.
  * Counting the number of "tires" (or rings), say $n$.
  * Estimating the radius for each tire, say $r_i$ for $1 \leq i \leq n$.
  * Estimating the height for each tire, say $h_i$ for $1 \leq i \leq n$


Then the volume would be found by adding:


$$
V = \pi \cdot r_1^2 \cdot h_1 + \pi \cdot r_2^2 \cdot h_2 + \cdots + \pi \cdot r_n^2 \cdot h_n.
$$

The weight would come by multiplying the volume by some appropriate density.


Looking at the sum though, we see the makings of an approximate integral. If the heights were to get infinitely small, we might expect this to approach something like $V=\int_a^b \pi r(h)^2  dh$.


In fact, we have in general:


> **Volume of a figure with a known cross section**: The volume of a solid with known cross-sectional area $A_{xc}(x)$ from $x=a$ to $x=b$ is given by
>
> $V = \int_a^b A_{xc}(x) dx.$
>
> This assumes $A_{xc}(x)$ is integrable.



This formula is derived by approximating the volume by "slabs" with volume $A_{xc}(x) \Delta x$ and using the Riemann integral's definition to pass to the limit. The discs of the Michelin man are an example, where the cross-sectional area is just that of a circle, or $\pi r^2$.


## Solids of revolution


We begin with some examples of a special class of solids - solids of revolution. These have an axis of symmetry from which the slabs are then just circular disks.


Consider the volume contained in this glass, it will depend on the radius at different values of $x$:


```{julia}
#| hold: true
#| echo: false
imgfile = "figures/integration-glass.jpg"
caption = L"""

A wine glass oriented so that it is seen as generated by revolving a
curve about the $x$ axis. The radius of revolution varies as a function of $x$
between about $0$ and $6.2$cm.

"""
#ImageFile(:integrals, imgfile, caption)
nothing
```

![A wine glass oriented so that it is seen as generated by revolving a
curve about the $x$ axis. The radius of revolution varies as a function of $x$
between about $0$ and $6.2$cm.
](./figures/integration-glass.jpg)

If $r(x)$ is the radius as a function of $x$, then the cross sectional area is $\pi r(x)^2$ so the volume is given by:


$$
V = \int_a^b \pi r(x)^2 dx.
$$

:::{.callout-note}
## Note
The formula is for a rotation around the $x$-axis, but can easily be generalized to rotating around any line (say the $y$-axis or $y=x$, ...) just by adjusting what $r(x)$ is taken to be.

:::

::: {#fig-solid-of-revolution}

```{julia}
#| echo: false
plt = let
    gr()
	# Follow lead of # https://github.com/SigurdAngenent/WisconsinCalculus/blob/master/figures/221/09surf_of_rotation2.py
	# plot surface of revolution around x axis between [0, 3]
	# best if r(t) descreases

	rad(x) = 2/(1 + exp(x))
    trange = (0,3)
    θrange = (0, 2pi)
    viewp = [2,-2, 1]

    ##
	proj(X) = _proj(X, viewp)

	# surface of revolution
    surf(t, z) = [t, rad(t)*cos(z), rad(t)*sin(z)]

	# project the surface at (t, a=theta)
	psurf(t,z) = proj(surf(t,z))

	# create shape holding project disc
    drawdiscF(t) = Shape(invert([psurf(t, 2*i*pi/100) for i in 1:101])...)

	α = 1.0 # opacity

	line_style = (; line=(:black, 1))
	plot(; empty_style..., aspect_ratio=:equal)



	# by layering, we get x-axis as desired
	plot!(pline(viewp, [-1,0,0], [0,0,0]); line_style...)
	plot!(drawdiscF(0);   fill =(:lightgray, α))

	plot!(pline(viewp, [0,0,0], [1,0,0]); line_style...)
	plot!(drawdiscF(1);   fill =(:black, α)) #  black to lightgray gives thickness
	plot!(drawdiscF(1.1); fill=(:lightgray, α))

	plot!(pline(viewp, [1.1,0,0], [2,0,0]); line_style...)
	plot!(drawdiscF(2);   fill=(:lightgray, α))

	plot!(pline(viewp, [2,0,0], [3,0,0]); line_style...)
	plot!(drawdiscF(3);   fill=(:lightgray, α))

	plot!(pline(viewp, [3,0,0], [4,0,0]); line_style..., arrow=true, side=:head)
	plot!(pline(viewp, [0,0,0], [0,0,1.25]); line_style..., arrow=true, side=:head)


	tt = range(trange..., 30)

	curve = psurf.(tt, pi/2)
	plot!(curve; line=(:black, 2))

	f1 = [(t, _fold(psurf, t, 0, pi)) for t in tt]
	curve = [psurf(f[1], f[2]) for f in f1]
	plot!(curve; line=(:black,1))

	f2 = [(t, _fold(psurf, t, pi, 2*pi)) for t in tt]
    curve = [psurf(f[1], f[2]) for f in f2]
	plot!(curve; line=(:black,1))



    ## find bottom edge (t,θ) again
	tt = range(0, 3, 120)
	f1 = [(t, _fold(psurf, t, pi, 2*pi)) for t in range(trange..., 100)]

	# shade bottom by adding bigger density of lines near bottom
	for (i,f) ∈ enumerate(f1)
		λ = iseven(i) ? 6 : 4 # adjust density by have some lines only extend to 6
		isnan(f[1]) || isnan(f[2]) && continue
		curve = [psurf(f[1], θ) for θ in range(f[2] - 0.2*(λ - f[1]), f[2], 20)]
		plot!(curve; line=(:black, 1))
	end

	current()

end
plt
```

```{julia}
#| echo: false
plotly()
nothing
```


Illustration of a figure being rotated around the $x$-axis. The discs have approximate volume given by the area of the base times the height or $\pi r(x)^2 \Delta x$. (Figure ported from @Angenent.)

:::


For a numeric example, we consider the original Red [Solo](http://en.wikipedia.org/wiki/Red_Solo_Cup) Cup. The dimensions of the cup were basically: a top diameter of $d_1 = 3~ \frac{3}{4}$ inches, a bottom diameter of $d_0 = 2~ \frac{1}{2}$ inches and a height of $h = 4~ \frac{3}{4}$ inches.


The central axis is straight down. If we rotate the cup so this is the $x$-axis, then we can get


$$
r(x) = \frac{d_0}{2} + \frac{d_1/2 - d_0/2}{h}x = \frac{5}{4} + \frac{5}{38}x
$$

The volume in cubic inches will be:


$$
V = \int_0^h \pi r(x)^2 dx
$$

This is


```{julia}
d0, d1, h = 2.5, 3.75, 4.75
rad(x) = d0/2 + (d1/2 - d0/2)/h * x
vol, _ = quadgk(x -> pi * rad(x)^2, 0, h)
```

So $36.9 \text{in}^3$. How many ounces is that? It is useful to know that 1 [gallon](http://en.wikipedia.org/wiki/Gallon) of water is defined as $231$ cubic inches, contains $128$ ounces, and weighs $8.34$ pounds.


So our cup holds this many ounces:


```{julia}
ozs = vol / 231 * 128
```

Full it is about $20$ ounces, though this doesn't really account for the volume taken up by the bottom of the cup, etc.


If you are poor with units, `Julia` can provide some help through the `Unitful` package. Here the additional `UnitfulUS` package must also be included, as was done above, to access fluid ounces:


```{julia}
vol * u"inch"^3 |> us"floz"
```

Before Solo "squared" the cup, the Solo cup had markings that - [some thought](http://www.snopes.com/food/prepare/solocups.asp) - indicated certain volume amounts.


```{julia}
#| hold: true
#| echo: false
imgfile = "figures/red-solo-cup.jpg"
caption = "Markings on the red Solo cup indicated various volumes"
# ImageFile(:integrals, imgfile, caption)
nothing
```

![Markings on the red Solo cup indicated various volumes.](./figures/red-solo-cup.jpg)

What is the height for $5$ ounces (for a glass of wine)? $12$ ounces (for a beer unit)?


Here the volume is fixed, but the height is not. For $v$ ounces, we need to convert to cubic inches. The conversion is $1$ ounce is $231/128 \text{in}^3$.


So we need to solve $v \cdot (231/128) = \int_0^h\pi r(x)^2 dx$ for $h$ when $v=5$ and $v=12$.


Let's express volume as a function of $h$:


```{julia}
Vol(h) = quadgk(x -> pi * rad(x)^2, 0, h)[1]
```

Then to solve we have:


```{julia}
v₅ = 5
h5 = find_zero(h -> Vol(h) - v₅ * 231 / 128, 4)
```

and


```{julia}
v₁₂ = 12
h12 = find_zero(h -> Vol(h) - v₁₂ * 231 / 128, 4)
```

As a percentage of the total height, these are:


```{julia}
h5/h, h12/h
```

:::{.callout-note}
## Note
Were performance at issue, Newton's method might also have been considered here, as the derivative is easily computed by the fundamental theorem of calculus.

:::

##### Example


By rotating the line segment $x/r + y/h=1$ that sits in the first quadrant around the $y$ axis, we will generate a right-circular cone. The volume of which can be expressed through the above formula by noting the radius, as a function of $y$, will be $R = r(1 - y/h)$. This gives the well-known volume of a cone:


```{julia}
#| hold: true
@syms r h x y
R = r*(1 - y/h)
cone = integrate(pi*R^2, (y, 0, h))
```

It is not unusual to parameterize a cone by the angle $\theta$ it makes and the height. Since $r/h=\tan\theta$, this gives the formula $V = \pi/3\cdot h^3\tan(\theta)^2$.


The frustum of a cone is simply viewed as a cone with its top cut off. If the original height would have been $h_0$ and the actual height $h_1$, then the volume remaining is just $\int_0^{h_1}  \pi r(y)^2 dy$. We can see the formula for the frustum of the right cone:

```{julia}
@syms h0
frustum = integrate(pi*R^2, (y, 0, h0))
```

Simplifying, we can see this volume can be expressed using the ratio $(h_0/h_1)$:

```{julia}
frustum - cone * ( 3h0/h - 3(h0/h)^2 + (h0/h)^3) |> simplify
```



##### Example


[Gabriel's](http://tinyurl.com/8a6ygv) horn is a geometric figure of mathematics - but not the real world - which has infinite height, but not volume! The figure is found by rotating the curve $y=1/x$ around the $x$ axis from $1$ to $\infty$. If the volume formula holds, what is the volume of this "horn?"


```{julia}
radius(x) = 1/x
quadgk(x -> pi*radius(x)^2, 1, Inf)[1]
```

That is a value very reminiscent of $\pi$, which it is as $\int_1^\infty 1/x^2 dx = -1/x\big|_1^\infty=1$.


:::{.callout-note}
## Note
The interest in this figure is that soon we will be able to show that it has **infinite** surface area, leading to the [paradox](http://tinyurl.com/osawwqm) that it seems possible to fill it with paint, but not paint the outside.

:::

##### Example


A movie studio hand is asked to find a prop vase to be used as a [Ming vase](http://en.wikipedia.org/wiki/Chinese_ceramics) in an upcoming scene. The dimensions specified are for the outside diameter in centimeters and are given by


$$
d(h) = \begin{cases}
2 \sqrt{26^2 - (h-20)^2} & 0 \leq h \leq 44\\
20 \cdot e^{-(h - 44)/10} & 44 < h \leq 50.
\end{cases}
$$

If the vase were solid, what would be the volume?


We define `d` using a ternary operator to handle the two cases:


```{julia}
d(h) = h <= 44 ? 2*sqrt(26^2 - (h-20)^2) : 20 * exp(-(h-44)/10)
rad(h) = d(h)/2
```

The volume in cm$^3$ is then:


```{julia}
Vₜ, _ = quadgk(h -> pi * rad(h)^2, 0, 50)
```

For the actual shoot, the vase is to be filled with ash, to simulate a funeral urn. (It will then be knocked over in a humorous manner, of course.) How much ash is needed if the vase has walls that are 1/2 centimeter thick


We need to subtract $0.5$ from the radius and `h` then recompute:


```{julia}
V_int, _ = quadgk(h -> pi * (rad(h) - 1/2)^2, 1/2, 50)
```

A liter of volume is $1000 \text{cm}^3$. So this is about $68$ liters, or more than 15 gallons. Perhaps the dimensions given were bit off.


While we are here, to compute the actual volume of the material in the vase could be done by subtraction.


```{julia}
Vₜ - V_int
```

### The washer method


Returning to the Michelin Man, in our initial back-of-the-envelope calculation we didn't account for the fact that a tire isn't a disc, as it has its center cut out. Returning, suppose $R_i$ is the outer radius and $r_i$ the inner radius. Then each tire has volume


$$
\pi R_i^2 h_i - \pi r_i^2 h_i = \pi (R_i^2 - r_i^2) h_i.
$$

Rather than use $\pi r(x)^2$ for a cross section, we would use $\pi (R(x)^2 - r(x)^2)$.


In general we call a shape like the tire a "washer" and use this formula for a washer's cross section $A_{xc}(x) = \pi(R(x)^2 - r(x)^2)$.


Then the volume for the solid of revolution whose cross sections are washers would be:


$$
V = \int_a^b \pi \cdot (R(x)^2 - r(x)^2) dx.
$$

::: {#fig-washer-illustration}
```{julia}
#| echo: false
plt = let
    gr()
	# Follow lead of # https://github.com/SigurdAngenent/WisconsinCalculus/blob/master/figures/221/09surf_of_rotation2.py
	# plot surface of revolution around x axis between [0, 3]
	# best if r(t) descreases

	rad(x) = 2/(1 + exp(x))
    trange = (0, 3)
    θrange = (0, 2pi)
    viewp = [2,-2,1]

	##
	proj(X) = _proj(X, viewp)


	# surface of revolution
    surf(t, z) = [t, rad(t)*cos(z), rad(t)*sin(z)]
	surf2(t, z) = (t, rad(t)*cos(z)/2, rad(t)*sin(z)/2)

	# project the surface at (t, a=theta)
	psurf(t,z) = proj(surf(t,z))
	psurf2(t, z) = proj(surf2(t,z))

	# create shape holding project disc
    drawdiscF(t) = Shape(invert([psurf(t, 2*i*pi/100) for i in 1:101])...)
	drawdiscI(t) = Shape([psurf2(t, 2*i*pi/100) for i in 1:101])

	α = 1.0
	line_style = (; line=(:black, 1))

	plot(; empty_style..., aspect_ratio=:equal)

	# by layering, we get x-axis as desired
	plot!(pline(viewp, [-1,0,0], [0,0,0]); line_style...)
	plot!(drawdiscF(0);   fill =(:lightgray, α))
	plot!(drawdiscI(0); fill=(:white, .5))

	plot!(pline(viewp, [0,0,0], [1,0,0]); line_style...)
	plot!(drawdiscF(1);   fill =(:black, α)) #  black to lightgray gives thickness
	plot!(drawdiscI(1); fill=(:white, .5))
	plot!(drawdiscF(1.1); fill=(:lightgray, α))
	plot!(drawdiscI(1.1); fill=(:white, .5))

    plot!(pline(viewp, [1.1,0,0], [2,0,0]); line_style...)
	plot!(drawdiscF(2);   fill=(:lightgray, α))
	plot!(drawdiscI(2); fill=(:white, .5))

    plot!(pline(viewp, [2,0,0], [3,0,0]); line_style...)
	plot!(drawdiscF(3);   fill=(:lightgray, α))
	plot!(drawdiscI(3); fill=(:white, .5))

    plot!(pline(viewp, [3,0,0], [4,0,0]); line_style..., arrow=true, side=:head)
	plot!(pline(viewp, [0,0,0], [0,0,1.25]); line_style..., arrow=true, side=:head)

    ## bounding curves
    ### main spine
	tt = range(trange..., 30)
	curve = [psurf(t, pi/2) for t in tt]
	plot!(curve; line=(:black, 2))

    ### the folds
	f1 = [(t, _fold(psurf, t, 0, pi)) for t in tt]
	curve = [psurf(f[1], f[2]) for f in f1]
	plot!(curve; line=(:black,))

	f2 = [(t, _fold(psurf, t, pi, 2*pi)) for t in tt]
    curve = [psurf(f[1], f[2]) for f in f2]
	plot!(curve; line=(:black,))

    ## add shading
    ### find bottom edge (t,θ) again
	f1 = [[t, _fold(psurf, t, pi, 2*pi)] for t in range(trange..., 120)]

	### shade bottom by adding bigger density of lines near bottom
	for (i,f) ∈ enumerate(f1)
		λ = iseven(i) ? 6 : 4 # adjust density by have some lines only extend to 6
		isnan(f[1]) || isnan(f[2]) && continue

		curve = [psurf(f[1], θ) for θ in range(f[2] - 0.2*(λ - f[1]), f[2], 20)]

		plot!(curve; line=(:black, 1))
	end

	current()


end
plt
```

```{julia}
#| echo: false
plotly()
nothing
```

Modification of earlier figure to show washer method. The interior volumn would be given by $\int_a^b \pi r(x)^2 dx$, the entire volume by $\int_a^b \pi R(x)^2 dx$. The difference then is the volume computed by the washer method.

:::

##### Example


An artist is working with a half-sphere of material, and wishes to bore out a conical shape. What would be the resulting volume, if the two figures are modeled by


$$
R(x) = \sqrt{1^2 - (x-1)^2}, \quad r(x) = x,
$$

with $x$ ranging from $x=0$ to $1$?


The answer comes by integrating:


```{julia}
#| hold: true
Rad(x) = sqrt(1 - (x-1)^2)
rad(x) = x
V, _ = quadgk(x -> pi*(Rad(x)^2 - rad(x)^2), 0, 1)
```

## Solids with known cross section


The Dart cup company now produces the red solo cup with a [square](http://www.solocup.com/products/squared-plastic-cup/) cross section. Suppose the dimensions are the same: a top diameter of $d_1 = 3 3/4$ inches, a bottom diameter of $d_0 = 2 1/2$ inches and a height of $h = 4 3/4$ inches. What is the volume now?


The difference, of course, is that cross sections now have area $d^2$, as opposed to $\pi r^2$. This leads to some difference, which we quantify, as follows:


```{julia}
#| hold: true
d0, d1, h = 2.5, 3.75, 4.75
d(x) = d0 + (d1 - d0)/h * x
vol, _ = quadgk(x -> d(x)^2, 0, h)
vol / 231 * 128
```

This shape would have more volume - the cross sections are bigger. Presumably the dimensions have changed. Without going out and buying a cup, let's assume the cross-sectional diameter remained the same, not the diameter. This means the largest dimension is the same. The cross section diameter is $\sqrt{2}$ larger. What would this do to the area?


We could do this two ways: divide $d_0$ and $d_1$ by $\sqrt{2}$ and recompute. However, each cross section of this narrower cup would simply be $\sqrt{2}^2$ smaller, so the total volume would change by $2$, or be 13 ounces. We have $26.04$ is too big, and $13.02$ is too small, so some other overall dimensions are used.


##### Example


For a general cone, we use this [definition](http://en.wikipedia.org/wiki/Cone):


> A cone is the solid figure bounded by a base in a plane and by a surface (called the lateral surface) formed by the locus of all straight line segments joining the apex to the perimeter of the base.



Let $h$ be the distance from the apex to the base. Consider cones with the property that all planes parallel to the base intersect the cone with the same shape, though perhaps a different scale. This figure shows an example, with the rays coming from the apex defining the volume.

::: {#fig-generic-cone}
```{julia}
#| echo: false
plt = let
	gr()
	rad(t) = 3/2 - t
    trange = (0, 3/2)
    θrange = (0, 2pi)
    viewp = [2,-1/1.5,1/2+.2]

    ##
	proj(X) = _proj(X, viewp)
	# our surface
    R, r, rho = 1, 1/4, 1/4
    f(t) = (R-r) * cos(t) + rho * cos((R-r)/r * t)
    g(t) = (R-r) * sin(t) - rho * sin((R-r)/r * t)
	surf(t, θ) = (rad(t)*f(θ), rad(t)*g(θ), t)
	psurf(t,θ) = proj(surf(t,θ))


	empty_style = (xaxis=([], false),
                    yaxis=([], false),
                    framestyle=:origin,
                    legend=false)
	axis_style = (arrow=true, side=:head, line=(:gray, 1))


	drawdiscF(t) = Shape([psurf(t, 2*i*pi/100) for i in 1:101])


	plot(; empty_style..., aspect_ratio=:equal)
	for (i,t) in enumerate(range(0, 3/2, 30))
		plot!(drawdiscF(t); fill=(:gray,1), line=(:black,1))
	end


	θ = 0; plot!([psurf(0, θ), psurf(3/2, θ)]; line=(:black, 2))
	θ = pi/2; plot!([psurf(0, θ), psurf(3/2, θ)]; line=(:black, 1))
	θ = 3pi/2; plot!([psurf(0, θ), psurf(3/2, θ)]; line=(:black, 1))


	current()
end
plt
```

```{julia}
#| echo: false
plotly()
nothing
```


A "cone" formed from the parameterized curve
$r(t) = \langle
(R-r) \cdot \cos(t) + \rho \cdot \cos((R-r)/r \cdot t),
(R-r) \cdot \sin(t) - \rho \cdot \sin((R-r)/r \cdot t)
\rangle$ with apex at the point $[0,0,3/2]$ and rays extending down through the origin following $3/2-z$.

:::



```{julia}
#| echo: false
#| eval: false
plt = let
    h = 5
    R, r, rho = 1, 1/4, 1/4
    f(t) = (R-r) * cos(t) + rho * cos((R-r)/r * t)
    g(t) = (R-r) * sin(t) - rho * sin((R-r)/r * t)
    ts = range(0, 2pi, length=100)

    plot(f.(ts), g.(ts), zero.(ts), legend=false)
    for t ∈ range(0, 2pi, length=25)
        plot!([0,f(t)], [0,g(t)], [h, 0], linecolor=:red)
    end
    current()
end
plt
```

A right circular cone is one where this shape is a circle. This definition can be more general, as a square-based right pyramid is also such a cone. After possibly reorienting the cone in space so the base is at $u=0$ and the apex at $u=h$ the volume of the cone can be found from:


$$
V = \int_0^h A_{xc}(u) du.
$$

The cross sectional area $A_{xc}(u)$ satisfies a formula in terms of $A_{xc}(0)$, the area of the base:


$$
A_{xc}(u)  = A_{xc}(0) \cdot (1 - \frac{u}{h})^2
$$

So the integral becomes:


$$
V = \int_0^h A_{xc}(u) du = A_{xc}(0) \int_0^h (1 - \frac{u}{h})^2 du = A_{xc}(0) \int_0^1 v^2 h dv = A_{xc}(0) \frac{h}{3}.
$$

This gives a general formula for the volume of such cones.

::: {#fig-cross-sections}

```{julia}
#| echo: false
plt = let
    gr()
	# sections
	# https://github.com/SigurdAngenent/WisconsinCalculus/blob/master/figures/221/09Xsections.py

	x(h,z) =  0.3*h^2+(0.6-0.2*h)*cos(z)
	y(h,z) =  h+(0.3-0.2*h)*sin(z)+0.05*sin(4*z)
	r(h,z) = (x(h,z), y(h,z))
	r1(h,z) = (2,0) .+ r(h,z)

    empty_style = (xaxis=([], false),
                   yaxis=([], false),
                   framestyle=:origin,
                    legend=false)

	Nh=30
	heights = range(-1/2, 1/2, Nh)
	h0=heights[Nh ÷ 2]
	h1=heights[Nh ÷ 2 + 1]
	hs = [heights[1], h0, h1, heights[end]]
	ts = range(0, 2pi, 300)

	plot(; empty_style..., aspect_ratio=:equal)

	# stack the curves
	for h in heights
		curve = r.(h, ts)
		plot!(Shape(curve); fill=(:white, 1.0), line=(:black, 1))
	end

	# shape pull outs; use black to give thickness
	for (h, color) in zip(hs, (:white, :black, :white, :white))
	    curve = r1.(h,ts)
		plot!(Shape(curve); fill=(color,1.0), line=(:black, 1,))
	end

	# axis with marked points
	plot!([(-1,-1), (-1, 1)]; axis_style...)
	pts = [(-1, y(h, pi)) for h in hs]
	scatter!(pts, marker=(5, :circle))

	# connect with dashes
	for h in hs
		plot!([(-1, y(h, pi)), r(h,pi)]; line=(:black, 1, :dash))
		plot!([r(h,0), r1(h,pi)]; line=(:black, 1, :dash))
	end





	current()


end
plt
```

This figure shows the volume of a figure being comprised of slices. A discrete approximation would be found by estimating the volume of each slice by the cross sectional area times a small $\Delta h$. This leads to a formula
$V = \int_a^b A(h)dh$, where $A$ computes the cross sectional area.
(This figure was ported from @Angenent.)
:::

```{julia}
#| echo: false
plotly()
nothing
```

### Cavalieri's method


[Cavalieri's](http://tinyurl.com/oda9xd9) Principle is "Suppose two regions in three-space (solids) are included between two parallel planes. If every plane parallel to these two planes intersects both regions in cross-sections of equal area, then the two regions have equal volumes." (Wikipedia).


::: {#fig-Cavalieris-first}

```{julia}
#| echo: false
plt = let
    gr()

    x(h,z) =  (0.6-0.2*h) * cos(z)
    y(h,z) =  h + (0.2-0.15*h) * sin(z) + 0.01 * sin(4*z)
    xa(h,z) =  2 + 0.1 * cos(7*pi*h) + (0.6-0.2*h)*cos(z)


    heights = range(-1/2, 1/2, 50)
    ts = range(0, 2pi, 300)
    h0 = heights[25]
    h1 = heights[26]

    plot(; empty_style..., aspect_ratio=:equal)

    for h in heights
        curve=[(x(h, t), y(h, t)) for t in ts]
        plot!(Shape(curve); fill=(:white,), line=(:black,1))

        curve=[(xa(h, t), y(h, t)) for t in ts]
        plot!(Shape(curve); fill=(:white,), line=(:black,1))
    end
    current()
end
plt
```

```{julia}
#| echo: false
plotly()
nothing
```

Illustration of Cavalieri's first principle. The discs from the left are moved around to form the left volume, but as the volumes of each  cross-sectional disc remains the same, the two valumes are equally approximated. (This figure ported from @Angenent.)

:::

With the formula for the volume of solids based on cross sections, this is a trivial observation, as the functions giving the cross-sectional area are identical. Still, it can be surprising.


Consider a sphere with an interior cylinder bored out of it.  (The [Napkin](http://tinyurl.com/o237v83) ring problem.) The bore has height $h$ - for larger radius spheres this means very wide bores.

::: {#fig-napkin-ring-1}
```{julia}
#| echo: false
plt = let
	# Follow lead of # https://github.com/SigurdAngenent/WisconsinCalculus/blob/master/figures/221/09surf_of_rotation2.py
	# plot surface of revolution around x axis between [0, 3]
	# best if r(t) descreases

	rad(t) = (t = clamp(t, -1, 1); sqrt(1 - t^2))
	rad2(t) = 1/2
    viewp = [2,-2,1]

    ##
    function _proj(X, v)
        # a is ⟂ to v and b is v × a
	    vx, vy, vz = v
	    a = [-vy, vx, 0]
	    b = cross([vx,vy,vz], a)
	    a′, b′ = a/norm(a), b/norm(b)
	    return 	(a′ ⋅ X, b′ ⋅ X)
    end

    # project a curve in R3 onto R2
    pline(viewp, ps...) = [_proj(p, viewp) for p in ps]

    # determinant of Jacobian; area multiplier
    # det(J); used to identify folds
    function jac(X, u, v)
	    return det(ForwardDiff.jacobian(xs -> collect(X(xs...)), [u,v]))
    end

    function _fold(F, t, θmin, θmax)
	    λ = θ -> jac(F, t, θ) # F is projected surface, psurf
	    iszero(λ(θmin)) && return θmin
	    iszero(λ(θmax)) && return θmax

	    return solve(ZeroProblem(λ, (θmin, θmax)))
    end

	##
    proj(X) = _proj(X, viewp)

	# surface of revolution about the z axis
    surf(t, z) = (rad(t)*cos(z), rad(t)*sin(z), t)
	surf2(t, z) = (rad2(t)*cos(z), rad2(t)*sin(z), t)

	# project the surface at (t, a=theta)
	psurf(t,z) = proj(surf(t,z))
	psurf2(t, z) = proj(surf2(t,z))

	bisect(f, a, b) = find_zero(f, (a,b), Bisection())

	# create shape holding project disc
    drawdiscF(t) = Shape([psurf(t, 2*i*pi/100) for i in 1:101])
	drawdiscI(t) = Shape([psurf2(t, 2*i*pi/100) for i in 1:101])


	α = 1.0
	line_style = (; line=(:black, 1))

	plot(; empty_style..., aspect_ratio=:equal)

	# washer
	t0 = sqrt(3/4)

	Δ = .03
	δ = 0.785398 + 0.05

	x₀ = -.25
	plot!(drawdiscF(x₀-Δ); fill=(:black,), line=(:black,1))
	plot!(drawdiscF(x₀); fill=(:orange,), line=(:black,1))
	plot!(drawdiscI(x₀); fill=(:white,1.0), line=(:black,1))

	x₀ = 0.35
	plot!(drawdiscF(x₀-Δ); fill=(:black,), line=(:black,1))
	plot!(drawdiscF(x₀); fill=(:orange,), line=(:black,1))
	plot!(drawdiscI(x₀); fill=(:white,1.0), line=(:black,1))

	z0 = 3pi/2 - δ
	plot!(pline(viewp, surf(t0, z0), surf(-t0, z0)); line=(:black, 1))
	plot!(pline(viewp, surf(t0, z0+pi), surf(-t0, z0+pi)); line=(:black, 1))

	# caps
	curve = [psurf(t0, θ) for θ in range(0, 2pi, 100)]
	plot!(curve, line=(:black, 2))
	curve = [psurf(-t0, θ) for θ in range(0, 2pi, 100)]
	plot!(curve, line=(:black, 2))


    ## folds
	tθs = [(t, _fold(psurf, t, 0,pi)) for t in range(-t0, t0, 50)]
	curve = [psurf(t, θ) for (t,θ) ∈ tθs]
	plot!(curve, line=(:black, 3))

	tθs = [(t, _fold(psurf, t, pi, 2pi)) for t in range(-t0, t0, 50)]
	curve = [psurf(t, θ) for (t,θ) ∈ tθs]
	plot!(curve, line=(:black, 3))



	# Shade lines
	δ = pi/6
	Δₜ =  (4pi/2 - (3pi/2 - δ))/(2*25)
	for θ ∈ range(3pi/2-δ, 4pi/2, 25)
		curve = [psurf(t, θ) for t in
				 range(-t0, max(-t0, -t0 + 1/2*sin(θ+δ+pi/2 + pi/2)), 20)]
		plot!(curve, line=(:black, 1))
		curve = [psurf(t, θ+Δₜ) for t in
				 range(-t0, max(-t0, -t0 + 1/3*sin(θ+δ+pi/2 + pi/2)), 20)]
		plot!(curve, line=(:black, 1))
	end

	#=
    f1 = [[t, _fold(psurf, t, 0, pi/2)] for t in range(-0.5, -0.1, 26)]
	for f in f1
  		plot!([psurf( f[1], f[2]-k*0.01*(6-f[1]) )
            for k in 1:21]; line=(:black, 1))
	end
	=#

	current()
end
plt
```

Figure showing sphere with interior cylinder bored out.
:::


This cross-sectional figure is used to better understand the key dimensions.

::: {#fig-napkin-ring-2}

```{julia}
#| hold: true
#| echo: false
#The following illustrates $R=1$ and $h=2sqrt(3/4$.
plt = let
    gr()
    empty_style = (xaxis=([], false),
                    yaxis=([], false),
                    framestyle=:origin,
                   legend=false)

    R = 1; h = 2*sqrt(3/4)

    θ = theta = asin(h/2/R)
    thetas = range(-theta, stop=theta, length=100)
    ts = range(-pi, stop=pi, length=100)
    y = h/4

    plot(; empty_style..., aspect_ratio=:equal)
    plot!(R*cos.(ts), R*sin.(ts); line=(:black,));
    plot!(R*cos.(thetas), R*sin.(thetas), line=(:orange,1));

    plot!([R*cos.(theta), R*cos.(theta)], [h/2, -h/2]; color=:orange);
    plot!([R*cos.(theta), sqrt(R^2 - y^2)], [y, y]; line=(:orange,3))

    plot!([0, R*cos.(theta)], [0,0], color=:red);
    plot!([ 0, R*cos.(theta)], [0,h/2], color=:red);

    x₀ = sqrt(R^2 - (h/2)^2)
    annotate!( [
        (x₀/2, 0, text(L"\sqrt{R^2- (\frac{h}{2})^2}",10, :top)),
        (x₀, h/4, text(L"\frac{h}{2}",:right)),
        (R/2*cos(θ),R/2*sin(θ), text(L"R", :bottom; rotation=rad2deg(θ)))
    ])
    current()
end
plt

```

```{julia}
#| echo: false
plotly()
nothing
```

Side view illustrating key dimensions of napkin ring problem with $R$ being the radius of the sphere and $h$ being the height of the resulting interior cylinder.

:::


The small orange line is rotated, so using the washer method we get the cross sections given by $\pi(r_0^2 - r_i^2)$, the outer and inner radii, as a function of $y$.


The outer radii has points $(x,y)$ satisfying $x^2 + y^2 = R^2$, so is $\sqrt{R^2 - y^2}$. The inner radii has a constant value, and as indicated in the figure, is $\sqrt{R^2 - (h/2)^2}$, by the Pythagorean theorem.


Thus the cross sectional area is


$$
\pi( (\sqrt{R^2 - y^2})^2 - (\sqrt{R^2 - (h/2)^2})^2 )
= \pi ((R^2 - y^2) - (R^2 - (h/2)^2))
= \pi ((\frac{h}{2})^2 - y^2)
$$

As this does not depend on $R$, and the limits of integration would always be $-h/2$ to $h/2$ by Cavalieri's principle, the volume of the solid will be independent of $R$ too.


To actually compute this volume, we take $R=h/2$, so that the bore hole is just a line of no volume, the resulting volume is then that of a sphere with radius $h/2$, or $4/3\pi(h/2)^3 = \pi h^3/6$.


## The second theorem of Pappus


The second theorem of [Pappus](http://tinyurl.com/l43vw4) says that if a plane figure $F$ is rotated around an axis to form a solid of revolution, the total volume can be written as $2\pi r A(F)$, where $r$ is the distance the centroid is from the axis of revolution, and $A(F)$ is the area of the plane figure. In short, the distance traveled by the centroid times the area.


This can make some computations trivial. For example, we can make a torus (or donut) by rotating the circle $(x-2)^2 + y^2 = 1$ about the $y$ axis. As the centroid is clearly $(2, 0)$, with $r=2$ in the above formula, and the area of the circle is $\pi 1^2$, the volume of the donut is $2\pi(2)(\pi) = 4\pi^2$.


##### Example


Above, we found the volume of a cone, as it is a solid of revolution, through the general formula. However, parameterizing the cone as the revolution of a triangle with vertices $(0,0)$, $(r, 0)$, and $(0,h)$ and using the formula for the center of mass in the $x$ direction of such a triangle, $r/3$, we get that the volume of a cone with height $h$ and radius $r$ is  $2\pi (r/3)\cdot (rh/2) = \pi r^2 h/3$, in agreement with the calculus based computation.


## Questions


###### Question


Consider this big Solo cup:


```{julia}
#| hold: true
#| echo: false
imgfile ="figures/big-solo-cup.jpg"
caption = " Big solo cup. "
# ImageFile(:integrals, imgfile, caption)
nothing
```

![Big solo cup.](./figures/big-solo-cup.jpg)

It has approximate dimensions: smaller radius 5 feet, upper radius 8 feet and height 15 feet. How many gallons is it? At $8$ pounds a gallon this would be pretty heavy!


Two facts are useful:


  * a cubic foot is 7.48052 gallons
  * the radius as a function of height is $r(h) = 5  + (3/15)\cdot h$


```{julia}
#| hold: true
#| echo: false
gft = 7.48052
rad(h) = 5 + (3/15)*h
a,err = quadgk(h -> pi*rad(h)^2, 0, 15)
val = a*gft
numericq(val, 1e1)
```

###### Question


In *Glass Shape Influences Consumption Rate* for Alcoholic [Beverages](http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0043007) the authors demonstrate that the shape of the glass can have an effect on the rate of consumption, presumably people drink faster when they aren't sure how much they have left. In particular, they comment that people have difficulty judging the half-finished-by-volume mark.


This figure shows some of the wide variety of beer-serving glasses:


```{julia}
#| hold: true
#| echo: false
imgfile ="figures/beer_glasses.jpg"
caption = "A variety of different serving glasses for beer."
#ImageFile(:integrals, imgfile, caption)
nothing
```

![A variety of different serving glasses for beer.](./figures/beer_glasses.jpg)

We work with metric units, as there is a natural relation between volume in cm$^3$ and liquid measure ($1$ liter = $1000$ cm$^3$, so a $16$-oz pint glass is roughly $450$ cm$^3$.)


Let two glasses be given as follows. A typical pint glass with linearly increasing radius:


$$
r(h) = 3 + \frac{1}{5}h, \quad 0 \leq h \leq b;
$$

and a curved-edge one:


$$
s(h) = 3 + \log(1 + h), \quad 0 \leq h \leq b
$$

The following functions find the volume as a function of height, $h$:


```{julia}
r1(h) = 3 + h/5
s1(h) = 2 + log(1 + h)
r_vol(h) = quadgk(x -> pi*r1(x)^2, 0, h)[1]
s_vol(h) = quadgk(x -> pi*s1(x)^2, 0, h)[1]
```

  * For the straight-sided glass find $h$ so that the volume is $450$.


```{julia}
#| echo: false
h450 = find_zero(h -> r_vol(h) - 450, 10)
numericq(h450)
```

  * For the straight-sided glass find $h$ so that the volume is $225$ (half full).


```{julia}
#| echo: false
h225 = find_zero(h -> r_vol(h) - 225, 10)
numericq(h225)
```

  * For the straight-sided glass, what is the percentage of the total height when the glass is half full. (For a cylinder it would just be 50.)


```{julia}
#| echo: false
numericq(h225/h450 * 100, 2, units="percent")
```

  * People often confuse the half-way by height amount for the half way by volume, as it is for the cylinder. Take the height for the straight-sided glass filled with $450$ mm, divide it by $2$, then compute the percentage of volume at the half way height to the original.


```{julia}
#| echo: false
numericq(r_vol(h450/2)/450*100, 2, units="percent")
```

---


  * For the curved-sided glass find $h$ so that the volume is $450$.


```{julia}
#| echo: false
h_450 = find_zero(h -> s_vol(h) - 450, 10)
numericq(h_450)
```

  * For the curved-sided glass find $h$ so that the volume is $225$ (half full).


```{julia}
#| echo: false
h_225 = find_zero(h -> s_vol(h) - 225, 10)
numericq(h_225)
```

  * For the curved-sided glass, what is the percentage of the total height when the glass is half full. (For a cylinder it would just be 50.)


```{julia}
#| echo: false
numericq(h_225/h_450 * 100, 2, units="percent")
```

  * People often confuse the half-way by height amount for the half way by volume, as it is for the cylinder. Take the height for the curved-sided glass filled with $450$ mm, divide it by $2$, then compute the percentage of volume at the half way height to the original.


```{julia}
#| echo: false
numericq(s_vol(h_450/2)/450*100, 2, units="percent")
```

###### Question


A right [pyramid](http://en.wikipedia.org/wiki/Pyramid_%28geometry%29) has its apex (top point) above the centroid of its base, and for our purposes, each of its cross sections. Suppose a pyramid has square base of dimension $w$ and height of dimension $h$.


Will this integral give the volume:


$$
V = \int_0^h w^2 (1 - \frac{y}{h})^2 dy?
$$

```{julia}
#| hold: true
#| echo: false
yesnoq("yes")
```

What is the volume?


```{julia}
#| hold: true
#| echo: false
choices = [
"``1/3 \\cdot b\\cdot h``",
"``1/3 \\cdot w^2\\cdot h``",
"``l\\cdot w \\cdot h/ 3``"
]
answ = 2
radioq(choices, answ)
```

###### Question


An ellipsoid is formed by rotating the region in the first and second quadrants bounded by the ellipse $(x/2)^2 + (y/3)^2=1$ and the $x$ axis around the $x$ axis. What is the volume of this ellipsoid? Find it numerically.


```{julia}
#| hold: true
#| echo: false
f(x) = 3*sqrt( 1 - (x/2)^2 )
val, _ = quadgk(x -> pi * f(x)^2, -2, 2)
numericq(val)
```

###### Question


An ellipsoid is formed by rotating the region in the first and second quadrants bounded by the ellipse $(x/a)^2 + (y/b)^2=1$ and the $x$ axis around the $x$ axis. What is the volume of this ellipsoid? Find it symbolically.


```{julia}
#| hold: true
#| echo: false
choices = [
"``4/3 \\cdot \\pi a b^2``",
"``4/3 \\cdot \\pi a^2 b``",
"``\\pi/3 \\cdot a b^2``"
]
answ = 1
radioq(choices, answ)
```

###### Question


A solid is generated by rotating the region enclosed by the graph $y=\sqrt{x}$, the lines $x=1$, $x=2$, and $y=1$ about the $x$ axis. Find the volume of the solid.


```{julia}
#| hold: true
#| echo: false
Ra(x) = sqrt(x)
ra(x) = 1
a,b=1,2
val, _ = quadgk(x -> pi * (Ra(x)^2 - ra(x)^2), a,b)
numericq(val)
```

###### Question


The region enclosed by the graphs of $y=x^3 - 1$ and $y=x-1$ are rotated around the $y$ axis. What is the volume of the solid?


```{julia}
#| hold: true
@syms x
plot(x^3 - 1, 0, 1, legend=false)
plot!(x-1)
```

```{julia}
#| hold: true
#| echo: false
Ra(y) = cbrt(y+1)
ra(y) = y + 1
a,b = -1, 0
val, _ = quadgk(x -> pi * (Ra(x)^2 - ra(x)^2), a,b)
numericq(val)
```

###### Question


Rotate the region bounded by $y=e^x$, the line $x=\log(2)$ and the first quadrant ($x,y \geq 0$) about the line $x=\log(2)$.

```{julia}
#| echo: false
let
    p = plot(exp, -0.1, log(2.2); legend=false, ylim = (-0.3, 2.3))
    hline!(p, [0], color=:black)
    vline!(p, [0], color=:black)
    vline!(p, [log(2)], linestyle=:dash)
    xs = range(0, log(2), length=50)
    ys = exp.(xs)
    plot!(p, [0, xs..., log(2), 0], [0, ys..., 0, 0], linewidth=2)
    ts = pi/2 .+ range(pi/16, 2pi - pi/16, length=50)
    a = 0.15; b = a/(1.5)
    as(t) = log(2) + a*cos(t)
    bs(t) = (exp(log(2)) - 4b) + b*sin(t)
    plot!(p, as.(ts), bs.(ts), linewidth=1, color=:black, arrow=(:closed, 2.0))

    p

end
```

(Be careful, the radius in the formula $V=\int_a^b \pi r(u)^2 du$ is from the line $x=\log(2)$, further, the constraint $x \geq 0$ needs attention.)


```{julia}
#| hold: true
#| echo: false
a, b = 0, exp(log(2))
ra(y) = log(2) - max(0, log(y))
val, _ = quadgk(y -> pi * ra(y)^2, a, b)
numericq(val)
```

###### Question


Find the volume of rotating the region bounded by the line $y=x$, $x=1$ and the $x$-axis around the line $y=x$. (The Theorem of Pappus is convenient and the fact that the centroid of the triangular region lies at $(2/3, 1/3)$.)


```{julia}
#| hold: true
#| echo: false
cm=[2/3, 1/3]
c = [1/2, 1/2]
rr = norm(cm - c)
A = 1/2 * 1 * 1
val = 2pi * rr * A
numericq(val)
```

###### Question


Rotate the region bounded by the line $y=x$ and the function $f(x) = x^2$ about the line $y=x$. What is the resulting volume?


You can integrate in the length along the line $y=x$ ($u$ from $0$ to $\sqrt{2}$). The radius then can be found by intersecting the line perpendicular line to $y=x$ at $u$ to the curve $f(x)$. This will do so:


```{julia}
theta = pi/4  ## we write y=x as y = x * tan(pi/4) for more generality, as this allows other slants.

f(x) = x^2
𝒙(u) = find_zero(x -> u*sin(theta) - 1/tan(theta) * (x - u*cos(theta)) - f(x), (u*cos(theta), 1))
𝒓(u) = sqrt((u*cos(theta) - 𝒙(u))^2 + (u*sin(theta) - f(𝒙(u)))^2)
```

(Though in this case you can also find `r(u)` using the quadratic formula.)


With this, find the volume.


```{julia}
#| hold: true
#| echo: false
a, b = 0, sqrt(2)
val, _ = quadgk(u -> pi*𝒓(u)^2, a, b)
numericq(val)
```

---


Repeat (find the volume) only this time with the function $f(x) = x^{20}$.


```{julia}
#| hold: true
#| echo: false
a, b = 0, sqrt(2)
f(x) = x^20
xval(u) = find_zero(x -> u*sin(theta) - 1/tan(theta) * (x - u*cos(theta)) - f(x), (0,sqrt(2)))
rad(u) = sqrt((u*cos(theta) - xval(u))^2 + (u*sin(theta) - f(xval(u)))^2)
val, _ = quadgk(u -> pi*rad(u)^2, a, b)
numericq(val)
```
