Formula for Nicolosi Globular Projection What is the formula for the Nicolosi Globular projection? In other words, what is the function that maps the coordinate (φ, λ) on the sphere to the coordinate (x, y) on the plane? 
You can assume that there are two such functions, one for each of the two hemispheres of the projection, such that the origin (x = 0, y = 0) in both cases is in the center of the given hemisphere.
Example of a Nicolosi Globular projection:
http://map-projections.net/img/rot-w/nicolosi-globular-110w.png

Let $(\lambda,\phi)$ be the relative longitude/latitude of points on the globe with respect to some reference point $X$ on equator. Let
$$u = \frac{2\lambda}{\pi},\; v = \frac{2\phi}{\pi},\; c = \cos\phi,\; s = \sin\phi$$
Under the Nicolosi Globular probjection which maps the hemisphere "centered" at $X$ to unit disk, the parallels (locus of constant $\phi$) and meridians (locus of constant $\lambda$) get mapped to circular arcs.
For simplicity of discussion, let us first assume $0 < |\lambda|,|\phi| < \frac{\pi}{2}$.


*

*For the parallels, their intersections with the central meridian ($\lambda = 0$) and boundary meridians ($\lambda = \pm \frac{\pi}{2}$) are equally spaced. This means
for the parallel of a given $\phi$, the circular arc passes through three specific points $(0,v)$ and $(\pm c, s)$. Let $(0,q)$ be the center for this arc. The arc is determined by the equations:
$$x^2 + (y-q)^2 = 1 - 2sq + q^2 = (q-v)^2$$
This implies
$$q = \frac{1-v^2}{2(s-v)}\quad\text{ and }\quad y = s - \frac{\Delta}{2q}\tag{*1a}$$
where $\Delta = 1 - x^2 - y^2$.

*For the meridians, their intersections with the equator ($\phi = 0$) are equally spaced too. This means for the meridian of a given $\lambda$, the circular arc
passes through the point $(u,0)$ and the poles $(0,\pm 1)$. Let $(-p,0)$ be the center for this arc. The arc is determined by the equations:
$$(x+p)^2 + y^2 = 1 + p^2 = (p+u)^2$$
This implies
$$p = \frac{1-u^2}{2u}\quad\text{ and }\quad x = \frac{\Delta}{2p}\tag{*1b}$$
Substitute the expression of $x$ and $y$ into expression of $\Delta$, we get
$$\begin{align}
x^2 + y^2 + \Delta = 1
&\iff \left(\frac{\Delta}{2p}\right)^2 + \left(s - \frac{\Delta}{2q}\right)^2 + \Delta = 1\\
&\iff (p^2+q^2)\Delta^2 + 4p^2q(q-s) \Delta -4p^2q^2c^2 = 0
\end{align}
$$
Under the assumption $0 < |\lambda|,|\phi| < \frac{\pi}{2}$, this is a quadratic equation in $\Delta$ with a negative constant term. It has one positive and one negative root. Since we want the intersection of the parallel and meridian inside the unit disk, 
$0 < \Delta < 1$ and hence $\Delta$ is the positive root of this quadratic equation. This means
$$\Delta = -A + \sqrt{A^2 + B}\quad\text{ where }\quad
A = \frac{2p^2q(q-s)}{p^2+q^2},\;
B = \frac{4p^2q^2c^2}{p^2+q^2}\tag{*2}
$$
Combine all these arguments, we find for the general case where $|\lambda|,|\phi| \le \frac{\pi}{2}$, we can determine the intersection of the corresponding parallel and meridian by following procedure:


*

*if $\lambda = 0$ or $\phi = \pm \frac{\pi}{2}$, then  $(x,y) = (0,v)$  

*else if $\phi = 0$, then $(x,y) = (u,0)$  

*else if $\lambda = \pm \frac{\pi}{2}$, then  $(x,y) = (\text{sign}(\lambda) c,s)$

*else 
first compute $p,q$ by $(*1a/1b)$,
next compute $A,B,\Delta$ by $(*2)$,
and finally, compute $x,y$ by $(*1a/1b)$ again.
