I'm running into a problem trying to understand the work in the two following papers: Feasibility Study of a Quadrilateralized Spherical Cube Earth Data Base. and An Icosahedron-based Method for Pixelizing the Celestial Sphere.
Both papers are about performing an equal-area projection from the surface of a sphere to one face of a polyhedron. The first paper uses a cube, the second an icosahedron, but the general procedure is the same: You take a point on the sphere $(\theta,\phi)$, linearly project it to a tangent plane $(u,v)$, then remap it to another point on the same tangent plane $(x,y)$ such that the area element in the remapping ($dA_p=dxdy$) is equivalent to the area element on the surface of the sphere ($dA_s$).
I trust the first paper's results because it's NASA and as far as I can tell they've been using the system successfully for some time, however the mapping function is given by minimizing the residuals of a power-series approximation. The second paper seemed more attractive because it gives an exact solution to the problem, but I'm having trouble working through how they got it.
To begin, we have a unit sphere centered on the origin $$X=\sin\theta\cos\phi, \\Y=\sin\theta\sin\phi, \\Z=\cos\theta,$$ and we project a point on its upper hemisphere to the plane $Z=1$ (ie. a Gnomonic projection at the N pole) $$u=X/Z=\tan\theta\cos\phi, \\v=Y/Z=\tan\theta\sin\phi.$$ The relationship between an area element on the plane, $dA_p=du\,dv$, and the sphere, $dA_s=d\theta\,d\phi$, is given by $$du\,dv = \left|J\left(\frac{u,v}{\theta,\phi}\right)\right|d\theta\,d\phi, \\du\,dv = \left|\begin{matrix}\partial u/\partial\theta & \partial u/\partial\phi\\\partial v/\partial\theta & \partial v/\partial\phi\end{matrix}\right|d\theta\,d\phi=\frac{\sin\theta}{\cos^3\theta}d\theta\,d\phi=\sqrt{u^2+v^2}(1+u^2+v^2)d\theta\,d\phi.$$ (Solved with Mathematica.) And with that step I've deviated from both papers already. They say: $$dA_p=\frac{1}{\cos^3\theta}dA_s=(1+u^2+v^2)^{3/2}dA_s,$$ which is supposedly obvious, but not to me. Where have I gone wrong? Or did I just stumble across a mistake in these research papers?
My code:
u[th_, ph_] := Tan[th] Cos[ph]; v[th_, ph_] := Tan[th] Sin[ph];
D[{u[th, ph], v[th, ph]}, {{th, ph}}] // MatrixForm
J = Simplify[Det[D[{u[th, ph], v[th, ph]}, {{th, ph}}]]]
(*th=ArcTan[1,Sqrt[u^2+v^2]];ph=ArcTan[u,v];*)
J /. {th -> ArcTan[1, Sqrt[u^2 + v^2]]}