I think a large part of the difficulty here is the variation in
the magnitudes of your "normals" in general.
When you give a "normal" in the form $\langle -f_x, - f_y, 1\rangle$,
basically what you have is a radial projection of
the hemisphere parameterized by
$\langle \sin\phi \cos\theta, \sin\phi \sin\theta, \cos\phi\rangle$
(for $0 \leq \phi \leq \frac\pi2$) onto the plane $z = 1$.
That is, you have a projection that takes a point $P$ on the surface
of the hemisphere to a projected point $P'$ on the plane such that $P$
and $P'$ are collinear with the center of the sphere.
(This is also called a gnomonic projection of the sphere.)
The coordinates of your "normals" are the coordinates of the
projected points.
As long as $0 \leq \phi < \frac\pi2$,
each point of your sphere does in fact project onto that plane,
although as $\phi$ approaches $\frac\pi2$ the magnitudes of your
"normals" grow without bound.
And obviously the boundary of your hemisphere, where $\phi = \frac\pi2$,
does not project on to the plane $z = 1$ at all.
If you change your choice of $f$ as recommended in the answer by H.R.,
this problem goes away. All your normals then will have the same magnitude,
and they will be defined at all points of the hemisphere.
When you give your "normals" in the form
$\langle\sin^2\phi\cos\theta, \sin^2\phi\sin\theta, \sin\phi\cos\phi\rangle$,
you again have non-uniform magnitudes, but this time the magnitudes
go to zero as $\phi$ approaches zero.
In effect, you are radially projecting the hemisphere
onto a kind of degenerate torus given by $r = \sin\phi$.
But notice that the three components of your "normals" all have the
common factor $\sin\phi$. You can normalize the magnitudes of all your
"normals" (except for the case $\phi = 0$)
by multiplying by the scalar $1/\sin\phi$.
If you do this, you get the vectors
$\langle\sin\phi\cos\theta, \sin\phi\sin\theta, \cos\phi\rangle$,
that is, the coordinates of each vector (for $0 < \phi \leq \frac\pi2$)
are simply the coordinates of the point on the sphere.
If you define these vectors as the normals for all points such
that $0 < \phi \leq \frac\pi2$, you may see that you can use the same
formula to define the normal for $\phi = 0$ as well, and it works very nicely.
The formula ${\bf r}_\phi \times {\bf r}_\theta$
works nicely for the measurement of area (which is what it is used for
in the example you quoted) precisely because its magnitude does go to zero
as $\phi$ goes to zero and does so in the same
way that the area element $r \,d\phi\,d\theta$ goes to zero.
But as you observed, this property is not so desirable when you are
trying to construct a set of normal vectors rather than trying to
integrate some scalar function over area.
I suppose you got to $\displaystyle{\frac{-c}{\tan(\theta)}=a\,\cos(\phi)+b\,\sin(\phi)}$
Now this is equivalent to $\displaystyle{\frac{-c}{\sqrt{a^2+b^2}}\cot(\theta)=\frac{a}{\sqrt{a^2+b^2}}\,cos(\phi)+\frac{b}{\sqrt{a^2+b^2}}\,\sin(\phi)}$
The two numbers before $\cos(\phi),\sin(\phi)$ represent a point on the unit sphere, since its norm is $1$.
In other words it we can define a unique $\phi_0$ such that $\begin{cases}\cos\phi_0=\frac{a}{\sqrt{a^2+b^2}}\\\sin\phi_0=\frac{b}{\sqrt{a^2+b^2}}\end{cases}$
The equation becomes $\displaystyle{\frac{-c}{\sqrt{a^2+b^2}}\cot(\theta)=\cos(\phi_0)\cos(\phi)+\sin(\phi_0)\sin(\phi)}=\cos(\phi-\phi_0)$
Best Answer
The way I see it is you need to parametrize the solid bounded by the sphere of radius $5$ and the plane $y=-4$. Here is how I would proceed :
Use spherical coordinates as follows : let $y=\rho \cos \phi$, $z=\rho \sin \phi \cos \theta$ and $x=\rho \sin \phi \sin \theta$, such that the sphere has equation $$ \rho =5, $$ and the plane $y=-4$ has equation $$ \rho \cos \phi=-4 $$
Now things become easy. The projection of the solid in the $yz$ plane is the domain $$ D= \{(\rho,\phi)\;|\; \frac{4}{ \cos \phi} \le \rho \le 5 , \cos^{-1}\left(\frac{-4}{5}\right)\le \phi \le \pi\} $$ And it follows that $$ E = \{(\rho,\phi,\theta )\;|\;(\rho,\phi) \in D, 0 \le \theta \le 2\pi \} $$