Why are the body axes of an axisymmetric body principal axes of inertia? I am reading "Vector Mechanics for Engineers, Statics | Dynamics", 12th edition, by Beer, Johnston, Mazurek, Cornwell, Self. 
In section 18.3C Motion of an Axisymmetric Body Under No Force, the authors write (you can skip the text and focus only on the bold sentence that I confused about):


We can now analyze the motion of an axisymmetric body about its mass center under no force except its own weight. Examples of such motion are furnished by projectiles (if air resistance is neglected) and by satellites and space vehicles after the burnout of their launching rockets.

The sum of the moments of the external forces about the mass center G of the body is zero, so $\dot H_G$ = 0 . It follows that the angular momentum $H_G$ of the body about G is constant. Thus, the direction of $H_G$ is fixed in space and can be used to define the Z axis, or axis of precession (Fig. 18.20). Let us select a rotating system of axes Gxyz with the z axis along the axis of symmetry of the body, the x axis in the plane defined by the Z and z axes, and the y axis pointing away from you (Fig. 18.21). This gives us: 
$$ H_x = −H_Gsin θ  ;   H_y = 0   ;    H_z = H_G cos θ $$(18.46)
where θ represents the angle formed by the Z and z axes, and $H_G$ denotes the constant magnitude of the angular momentum of the body about G. Because the x, y, and z axes are principal axes of inertia for the body considered, we have
$$ H_x = I′ω_x  ;  H_y = I′ω_y   ;  H_z = Iω_z  $$(18.47)
where I denotes the moment of inertia of the body about its axis of symmetry and I′ denotes its moment of inertia about a transverse axis through G. It follows from Eqs. (18.46) and (18.47) that
$$ \omega_x = \frac{H_Gsin\theta}{I'} ; \omega_y = 0 ; \omega_z = \frac{H_Gcos\theta}{I}$$(18.48)
The second of these relations shows that the angular velocity ω has no component along the y axis; that is, along an axis perpendicular to the Z-z plane. Thus, the angle θ formed by the Z and z axes remains constant and the body is in steady precession about the Z axis. 

I thought principal axes give $I_{xy} = I_{yx} = I_{xz} = I_{zx} = I_{yz} = I_{zy} = 0$, but look at the figure 18.20: we only have $I_{xy} = I_{yx} = I_{yz} = I_{zy} = 0$ and I think $I_{xz} = I_{zx} \neq 0$ So why are the x, y, and z axes principal axes of inertia for the body considered ? 
 A: I have an intuitive understanding of why this is so, but I cannot articulate it. Instead, I am offering a mathematical proof of this, which unfortunately is rather dry in the insight department.
Summary
When you compute the mass moment of inertia tensor in cylindrical coordinates, the off-diagonal terms contain multiples of $\sin \theta \cos \theta$ and the diagonal terms contain multiples of $\sin^2 \theta $ or $\cos^2 \theta$. When integrated over $0 \ldots 2\pi$ the off-diagonal terms vanish, and the diagonal terms remain.
General Example
An axisymmetric body with the axis of symmetry along the z-axis has a boundary defined by a general function $f(z)$, such that the location of each point inside the body has radial distance between $0$ and $f(z)$.
$$\begin{aligned} {\rm pos} & = \pmatrix{x \\ y \\ z} = \pmatrix{ r \cos \theta \\ r \sin \theta \\ z} & 
& \begin{aligned} r & = 0 \ldots f(z) \\ \theta & = 0 \ldots 2\pi \\ z & = 0 \ldots \ell \end{aligned} \end{aligned} $$
The volume is derived by the integral
$$V = \int \limits_0^\ell \int \limits_0^{2\pi} \int \limits_0^{f(z)}\,r\;{\rm d}r \,{\rm d}\theta \,{\rm d}z\; = \pi \int \limits_0^\ell f^2(z)\,{\rm d}z$$
and the density is thus $\rho = {m \over V}$ which is used next in the MMOI integral
$$ \begin{aligned} \mathbf{I} & = \int \limits_0^\ell \int \limits_0^{2\pi} \int \limits_0^{f(z)} \rho \,r \begin{vmatrix} y^2+z^2 & -x y & -x z \\ -x y & x^2+z^2 & -y z\\ -x z & -y z & x^2+y^2 \end{vmatrix} \;{\rm d}r \,{\rm d}\theta \,{\rm d}z \\
  & =  \int \limits_0^\ell \int \limits_0^{f(z)} \rho \,r  \int \limits_0^{2\pi} \begin{vmatrix} z^2+r^2\sin^2\theta & -r^2 \sin\theta & -r z \cos \theta \\ -r^2 \sin\theta\cos\theta & z^2+r^2 \cos^2 \theta & -r z \sin \theta \\ -r z \cos\theta &  -r z \sin \theta & r^2
\end{vmatrix} \;{\rm d}\theta\,{\rm d}r \,{\rm d}z \\
& = \int \limits_0^\ell \int \limits_0^{f(z)} \rho \,r   \begin{vmatrix} \pi ( r^2+2 z^2) & & \\ & \pi ( r^2+2z^2) & \\ & & 2 \pi r^2
\end{vmatrix} \;{\rm d}r \,{\rm d}z 
\end{aligned}$$
and there you have it, the off-diagonal terms vanish in the integration around the shape. This means that the axis of revolution has to be a principal axis of inertia.
But why? 
Basic Reasoning
Think of this in reverse. The axis that happens to equally distribute the mass around has to be a principal axis. There is no lopsidedness in the distribution of mass. And by definition, a rotation axis has to distribute the mass equally around it. 
A: First, because the object is axisymmetric so $I_{zy} = I_{zx}$ (y and x axes are similar). We can not have  $I_{zy} = 0 $ and $I_{zx} \neq 0$.
Second, after seeing closely, similary to $I_{zy}$, I find that $I_{zx} = 0$.
