How do I apply the divergence theorem to this region? I am trying to find the flux of a given field ($\vec{F} = (x+yz^3)\hat{\textbf{i}} + (xz^3)\hat{\textbf{j}} -z\hat{\textbf{k}}$ over a surface defined by $x^2 +y^2 + z^2 = 9$ between the planes $z=1$ and $z=2$.
My vector calculus is quite rusty. I think the optimal way to do this is use the divergence theorem over the solid region, and I want to check that I'm setting up the integral correctly and also see if there are other ways to go about this (to build more intuition for the geometry of these situations).
My work:
I think setting up the volume integral will be easiest in cylindrical coordinates. For a given $z$, $x^2+y^2+z^2 = 9 \Rightarrow ~r^2 = 9-z^2 \Rightarrow r = \sqrt{9-z^2}$.
So I should be able to set up my integral as:
$$\int _1 ^2 \int_0 ^{\sqrt{9-z^2}} \int_0^{2\pi} \nabla \cdot \vec{F}~ r d\theta dr dz$$
Then I can simply use the standard substitution $x = r\cos \theta ~ y = r\sin \theta, z = z$ on the vector field to get 
$\vec{F} = (r\cos \theta+ r\sin \theta~ z^3) \hat{\textbf{i}} + r\cos \theta~z^3 \hat{\textbf{j}} - z~\hat{\textbf{k}}$.
a) Is this correct?
b) Is this the optimal way to have done it (in terms of ease)
c) Are there other set ups that would work without complicating things significantly?
 A: Divergence theorem can not be used right away as the surface is not closed.
Using the divergence theorem,
Close the surface by including the $z=1$ and $z=2$ that would make top and bottom parts as parts of the surface. So then the theorem can be applied.
$$\iint_{S} \vec F \cdot d \vec S=\iint \text{div} \vec F dV$$
Clearly $\text{div} F=0$ in this example, so the surface integral is $0$.
But $S=S_{\text{top}} \bigcup S_{\text{interest}} \bigcup S_{\text{bottom}}$.
So,
$$\iint_{S_{\text{top}}}+\iint_{S_{\text{interest}}}+\iint_{S_{\text{bottom}}}=0$$
And hence the surface integral of interest is equal to,
$$\iint_{S_{\text{interest}}}=-\iint_{S_{\text{bottom}}}-\iint_{S_{\text{top}}}$$
Note they all must be computed using an outward orientation, outward relative to the sphere.
For $z=1$ we have a downward outward orientation, so use $d \vec S=\langle 0,0,-1 \rangle dA$. For $z=2$ we have an upward outward orientation, so use $d \vec S= \langle 0,0,1 \rangle dA$. Note for $z=1$ we have $\vec F=\langle ...,...,-1 \rangle$. And for $z=2$ we have $\vec F= \langle ...,....,-2 \rangle$.
So the integral we need to evaluate is
$$-\iint_{D_1} \langle ...,...,-1 \rangle \cdot \langle 0,0,-1 \rangle dA -\iint_{D_2} \langle ...,....,-2 \rangle \cdot \langle 0,0,1 \rangle dA$$
Or,
$$-A(D_1)+2A(D_2)$$
$$=-\pi(\sqrt{8})^2+2\pi(\sqrt{5})^2$$
$$=10\pi-8\pi=2\pi$$
A: On a sphere of radius $a$,
$$\vec n=\frac{1}{a}<x,y,z>$$
Then,
$$\iint_{S} \vec F \cdot d \vec S=\iint_{S} \frac{1}{a} \vec F \cdot <x,y,z> dS$$
In your case we have,
$$\iint_{S} \vec F \cdot d \vec S=\frac{1}{a}\iint_{S} (x^2+2xyz^3-z^2) dS$$
Note $ds=a^2 \sin (\phi) dA$. This can be directly calculated by parametrizing the sphere by $r(\theta, \phi)$ and calculating $|r_\theta \times r_\phi|$. But can also be easily seen by the fact that we have to factor in a Jacobian $\rho^2 \sin (\phi)$ whenever switching to spherical coordinates.
What are the bounds for $\theta$ and $\phi$?
It's pretty clear $\theta \in [0,2\pi]$ 
The bounds for $\phi$ however take some work.
Note we have $r^2+z^2=9$ so when $z=1$ we have $r^2=8$ and hence we have $r=\sqrt{8}$.
Note that have $r^2+2^2=9$ when $z=2$. So $r=\sqrt{5}$ there.
$\frac{r}{z}=\tan (\phi)$
So it follows,
$$\phi \in [\arctan (\frac {\sqrt{5}}{2}), \arctan (\sqrt{8})]$$
Messy but this is how you would set it up using spherical coordinates. 
