How is the spherical elevation angle distributed when $(x,y,z)$ are uniformly and normally chosen? As a follow up to
How the polar coordinate, $\theta$, is distributed when $(x,y) \sim U(-1,1) \times U(-1,1)$ and if $(x,y) \sim N(0,1)\times N(0,1)$?
Assume $(x,y,z) \sim U(-10,10) \times U(-10,10) \times U(-10,10)$ how are $\theta$ and $\phi$ distributed?
It's clearly from the wonderful answers of the previous question that $\theta$ looks like that:

But why the $\phi$ doesn't get maximum likelihood at $\phi = \pi/4$?

If we select $x,y,z$ in normally distribution fashion we get this 2 p.d.f's:

Is there a name for $\theta$ and $\phi$ distributions in both cases? For me it looks like $\beta$ distribution  on the interval $[-90,90]$. 
 A: In my discussion here I am assuming your $\theta$ is effectively a longitude and $\phi$ is effectively a latitude. Perhaps more typical spherical co-ordinates use an angle down from the north pole rather than up from the equator and swap the roles of the two symbols from that - but it's no problem to deal with it either way, so I'll stick with what your notation appears to be.

Note that the distribution of radius is not of interest here, only the angles, so we could project everything onto a unit sphere without changing the angles. This is quite useful in the normal case.
With a spherically symmetric distribution like the three dimensional standard normal, the appearance of the distribution of inclination has to do with the fact that there's a lot more area on the surface of a sphere near the equator than near the poles. 

If you follow through the mathematics (or write a geometric argument in terms of elements of probability similar to the earlier 2D question), you can get that the inclination should have a density proportional to $\cos(\phi)$. Here's a geometrical argument that should motivate it in the "elements of probability" terms:

Since the radius at the equator is 1 and the radius at latitude $\phi$ is $\cos(\phi)$, the circumference at latitude $\phi$ is proportional to $\cos(\phi)$, and so the density at $\phi$ is proportional to $\cos(\phi)$.

Uniform case: With the 3D-uniform normalized to constant radius, you don't have uniformity of density on the sphere for the same reason that we didn't in the 2D case - when you project onto the sphere, there's a lot more "density" on the sphere near the angles where the corners are than where the sides are (with parts near the middle of the edges being in between) -- because there's more of the volume of the cube for angles close to the corners than for angles near the middle of the faces. 
We can see this by generating many random values uniformly in the cube and projecting them onto the sphere. Since there's more volume near the corners than near the faces of the cube, there's greater density looking "inward" from the corners than the faces. If we plot the height (recall this is a projected-z-value, $z^* = z/r$, where $r=\sqrt{x^2+y^2+z^2}$) above the equator against the longitude, we get the top plot below:

That height corresponds to the vertical side of the right triangle in the previous diagram; that height is the $\sin$ of $\phi$ ($z^*=\sin(\phi)$). To convert that to the latitude ($\phi$), we would take the arcsin of that projected vertical height, which is what we see in the lower plot. This "stretches" things more the closer we get to the pole, making the density as a function of latitude drop to 0 at the north and south pole (for both the uniform and for the normal case).
The density for $\phi$ will then be the integral of that bivariate density over $\theta$.

Looking at the marginal for $\theta$ (i.e. strips running down at fixed values of $\theta$) makes for four peaks in the density of $\theta$ as you note - indeed this follows directly from the 2D case, but as we now see, it also makes for a pair of peaks in the density of $\phi$ away from the equator, corresponding to a region on the surface of the unit sphere where the corners and upper/lower edges of the cube project.
A: The complementary cumulative distribution for the spherical latitude $\phi$ gives the chance that a random point in the cube $[-1,1]^3$ will lie above the cone that graphs the function $z = \cot(\phi)\sqrt{x^2+y^2}$.  Because these points are uniformly distributed throughout the cube (which has a volume of $8$), this chance is one-eighth the volume between the cone and the top of the cube.  When the latitude exceeds $\pi/4$, this volume is that of a right cone with height $1$ and base $\cot(\phi)$, equal to
$$F_{+}(\phi) = \frac{1}{8}\frac{\pi}{3}\cot^2(\phi).$$
See the two righthand plots in the figure.
When the latitude is less than $\arctan(1/\sqrt{2})$, this is the volume of the intersection of a semi-infinite cone and the cube. An integration in polar coordinates gives the expression
$$F_{-}(\phi) = \frac{1}{8}\left(4-\frac{4}{3} \tan (\phi ) \left(\sqrt{2}+2 \tanh ^{-1}\left(\tan \left(\frac{\pi }{8}\right)\right)\right)\right).$$
See the two leftmost plots in the figure.

The negative derivatives of these expressions give the density.  Between $\arctan(1/\sqrt{2})\approx \pi/5$ and $\pi/4$ is a transition region where the intersection of the cone with the cube is complicated.  Although an exact expression could be developed, it would be messy.  What we do know is that the density must change continuously from the derivative of $-F_{-}$ to the derivative of $-F_{+}$ as $\phi$ varies between those points.  This is shown in a histogram of a million simulated values (from the upper half of the cube only--the lower half will be a mirror image).  The gold curve is the graph of $-\frac{d}{d\phi}F_{-}$ while the red curve at the right is the graph of  $-\frac{d}{d\phi}F_{+}.$ 

This clarifies why the modes are not at $\phi=\pm \pi/4$, but must lie between these values and $\pm \arctan(1/\sqrt{2})$.
