I'm trying to parameterize a sphere so it has **6 faces of equal area**, like this:

But this is the closest I can get (simply jumping $\frac{\pi}{2}$ in $\phi$ azimuth angle for each "slice").

I can't seem to get the $\theta$ elevation parameter correct. Help!

Parametrizing squares with circular coordinates is doomed to bring pain!

Parametrizing a square with a square is much easier:

```
up = NestList[RotateLeft, {s, t, 1}/Sqrt[s^2 + t^2 + 1], 2];
ParametricPlot3D[
Join[up, -up],
{s, -0.9, 0.9}, {t, -0.9, 0.9},
PlotStyle -> {Red, Blue, Yellow, Green, Magenta, Cyan},
Boxed -> False, Axes -> False
]
```

produces

Notice I draw from -0.9 to 0.9 in order to emphasize the patches, but you want to do it from -1 to 1.

(If you want the parametrizations to have their obvious normals going out of the sphere, then you need to be more gentle when constructing the last 3 patches: I just changed the sign of the other three, but that gets them oriented the wrong way..

Rotate the figure so that the $z$-axis passes through two opposite vertices, $N$ and $S$. Then the edges ---three "northern" and three "southern"--- adjacent to either $N$ or $S$ lie on three uniformly-arranged planes through the $z$-axis. Let the complete great circles determined by these planes sub-divide your six "sides" into twelve congruent triangular regions. Then each region's "azimuth range" covers exactly one-sixth of a circle. That's the easy part.

The remaining six of the figure's original edges ---each separating one northern triangular region from one southern region--- form a "belt" that zig-zags above and below the equatorial $xy$-plane; let $\psi$ be the elevation angle at one of its highest points. One extreme of the "elevation range" of a region for a given azimuth is determined by this belt (the other extreme is $\pi/2$ for a northern region, and $-\pi/2$ for a southern region). That's the trickier part.

With an appropriate rotation about the $z$-axis, we can take the belt vertices to be $$P_k(\;\cos\psi \; \cos\frac{k\pi}{3},\;\cos\psi \; \sin\frac{k\pi}{3},\;(-1)^k\sin\psi\;)$$ for $k=0,1,2,3,4,5$. The great circle containing the belt-edge $P_k P_{k+1}$ is the intersection of the sphere with the plane through those vertices and the origin; thus, its the collection of sphere-points orthogonal to the normal to that plane, with azimuth-elevation coordinates $(\phi,\theta)$ satisfying this equation $$ \begin{bmatrix}\cos\theta\cos\phi \\ \cos\theta\sin\phi \\ \sin\theta \end{bmatrix} \cdot \left( P_k \times P_{k+1} \right) = 0 $$

(Tip of the hat to @joriki and his answer to your more-basic version of this question.) That is, $$ \begin{eqnarray*} 0 &=&(-1)^k \cos\psi \sin\psi \left(\sin\frac{k\pi}{3}+\sin\frac{(k+1)\pi}{3}\right) \cos\theta\cos\phi \\ &-& (-1)^k \cos\psi \sin\psi \left(\cos\frac{k\pi}{3}+\cos\frac{(k+1)\pi}{3}\right) \cos\theta\sin\phi \\ &-& \cos^2\psi \left(\sin\frac{(k+1)\pi}{3}\cos\frac{k\pi}{3}-\cos\frac{(k+1)\pi}{3}\sin\frac{k\pi}{3}\right) \sin\theta \end{eqnarray*} $$

This collapses to $$ \tan\theta = (-1)^k \; \left( 2 \tan\psi \right) \; \sin\left(\frac{(2k+1)\pi}{6} - \phi\right) $$ (Sanity check: At the endpoints $\phi_{+} = \frac{k\pi}{3}$ and $\phi_{-} = \frac{(k+1)\pi}{3}$ of a belt-edge, we get $\tan\theta_{\pm} = \pm (-1)^k \tan\psi$, so that $\theta_{\pm} = \pm(-1)^k\psi$, as desired. Also, at the midpoint $\phi_{\star} = \frac{(2k+1)\pi}{6}$, where the edge crosses the $xy$-plane, we get the correct (lack of) elevation, $\theta_{\star} = 0$.)

Now, the edge from $N$ to nearby belt-vertex (like any edge in the original figure) has arc-measure $\mathrm{acos}\frac{1}{3}$; our $\psi$ ---which measures from *equator* to belt-vertex--- is the complement of this angle, so that $\tan\psi = \frac{1}{2\sqrt{2}}$ and the final formula becomes
$$
\tan\theta = \frac{(-1)^k}{\sqrt{2}} \sin\left(\frac{(2k+1)\pi}{6} - \phi\right)
$$

This gives the non-obvious extreme of the elevation range for points in one of the triangular regions. Since these regions merely sub-divide your six "sides", you have a parameterization of those sides. (I'll leave it to you to transform the parameterization if you absolutely require the $z$-axis to pass through opposing sides rather than opposing vertices.)

It's easiest to work from the corners. There are 8 corners, regularly spaced; you can take them to be the centers of the 8 octants. So in Cartesian coordinates they would be $(\pm 1/\sqrt{3}, \pm 1/\sqrt{3}, \pm 1/\sqrt{3})$. In spherical coordinates, $\phi$ is $\pm \pi/4$ or $\pm 3\pi/4$, and $\theta$ is $\cos^{-1} (\pm 1/\sqrt{3})$. You then simply draw the appropriate edges.

The following doesn't have much to do with spherical coordinates, but it might be worth noting that these 6 regions can be seen as the projections of the 6 faces of an enclosing cube.

In other words, each of the 6 regions can be parametrized as the region of the sphere $$S=\{\{x,y,z\}\in\mathbb R^3\mid x^2+y^2+z^2=1\}$$ for which, respectively:

- $x > \max(\lvert y\rvert,\lvert z\rvert)$
- $x < -\max(\lvert y\rvert,\lvert z\rvert)$
- $y > \max(\lvert z\rvert,\lvert x\rvert)$
- $y < -\max(\lvert z\rvert,\lvert x\rvert)$
- $z > \max(\lvert x\rvert,\lvert y\rvert)$
- $z < -\max(\lvert x\rvert,\lvert y\rvert)$

For example, the top region can be parametrized by $z =\sqrt{1-x^2-y^2}$, with region \begin{align*} \lvert x\rvert&<\frac{1}{\sqrt{2}}\\ \lvert y\rvert&<\min\left(\sqrt{\frac{1-x^2}{2}},\sqrt{1-2x^2}\right) \end{align*}

Following on FelixCQ's approach, take the cube to have corners $\pm 1$ and we want to project onto the unit sphere. So one side is $(-1+2t,1,1)$ Projecting that onto the unit sphere gives $\frac{1}{\sqrt{3-4t+4t^2}}(-1+2t,1,1)$, $\theta=\arccos(\frac{1}{\sqrt{3-4t+4t^2}}, \phi=\arctan(-1+2t)$. The other sides can be treated similarly, but you have to worry about the quadrants for $\phi$.