Conformal mapping circle onto square (and back)

Question

I'm programming an implementation of the Peirce quincuncial map projection. The projection involves a stereographic projection of a hemisphere of the globe onto a circle (I've got that part), then mapping points on that circle onto a square with a conformal mapping.

Wikipedia describes the relationship between a point $(p, \theta)$ on the circle and a point $(x, y)$ on the square as

$$\tan \left( \frac{p}{2} \right) e^{i \theta} = \mathrm{cn} \left( z, \frac{1}{2} \right), \text{ where } z = x + i y.$$

I don't understand the notation $\mathrm{cn} \left( z, \frac{1}{2} \right)$. Can it be written using algebraic and trigonometric functions?

That is, can you rewrite the above like this?

$\tan \left( \frac{p}{2} \right) \cos\theta = $ something in terms of $x$ and $y$

$\tan \left( \frac{p}{2} \right) \sin\theta = $ something in terms of $x$ and $y$

Thanks.

Answer 1

The $\mathrm{cn}$ function is a Jacobi elliptic function of $z=x+iy$; to get $z$, you require the corresponding inverse Jacobi elliptic function. By analogy with trig functions, these go by $\mathrm{arccn}$, $\mathrm{arcsn}$, etc. For your case, we have $$ z = \mathrm{arccn}\left(\tan\left(\frac{p}{2}\right)e^{i\theta},\frac{1}{2}\right) $$

The inverse Jacobi elliptic function $\mathrm{arccn}$ has an expression in terms of an elliptic integral. In general, it is given by

$$ \mathrm{arccn}\left(a,k\right) = \int_{a}^1 \frac{1}{\sqrt{(1-t^2)(k'^2+k^2t^2)}} dt, $$

where $k' = \sqrt{1-k^2}$ (see here). In your case, $k=1/2$ and so $k'=\sqrt{3}/2$, and $a=\tan\left(\frac{p}{2}\right)e^{i\theta}$. So,

$$ z = \int_{\tan\left(\frac{p}{2}\right)e^{i\theta}}^1 \frac{1}{\sqrt{(1-t^2)\left(\frac{3}{4}+\frac{1}{4}t^2\right)}} dt $$

So, for a given $p$ and $\theta$, you have to crunch that integral numerically. Or, you maybe be able to find a good elliptical integral package that does this for you, for complex arguments. The real part would be your $x$, and the imaginary part your $y$. As a final note, I'm not sure about the value $k=1/2$ from the Wikipedia article; this article seems to use the value of $k=k'=\frac{1}{\sqrt{2}}$.

Answer 2

I've found this implementation to work, using the scipy package in python. The implementation is actually just the 'peirce_map' function at the bottom; the rest of it is to obtain a cn() function that works on complex numbers.

import numpy
from scipy.special import ellipj, ellipk

#
# The scipy implementation of ellipj only works on reals;
# this gives cn(z,k) for complex z. 
# It works with array or scalar input.
#
def cplx_cn( z, k):
    z = numpy.asarray(z)
    if z.dtype != complex:
        return ellipj(z,k)[1]

    # note, 'dn' result of ellipj is buggy, as of Mar 2015
    # see https://github.com/scipy/scipy/issues/3904
    ss,cc = ellipj( z.real, k )[:2]
    dd = numpy.sqrt(1-k*ss**2)   # make our own dn
    s1,c1 = ellipj( z.imag, 1-k )[:2]
    d1 = numpy.sqrt(1-k*s1**2)

    # UPDATE: scipy bug seems to have been fixed mid 2016, so 
    # four lines above could be done as these two, if you have that.
    #   ss,cc,dd = ellipj( z.real, k )
    #   s1,c1,d1 = ellipj( z.imag, 1-k )

    ds1 = dd*s1
    den = (1-ds1**2)
    rx = cc*c1/den
    ry = ss*ds1*d1/den
    return rx - 1j*ry
#
# Kval is the first solution to cn(x,1/2) = 0
# This is K(k) (where 4*K(k) is the period of the function).
Kval = ellipk(0.5) # 1.8540746773013719

#######################################################
# map a complex point in unit square to unit circle
# The following points are the corners of the square (and map to themselves):
#     1    -1     j    -j
#  The origin also maps to itself.
# Points which are in : abs( re(z)) <=1, abs(im(z)) <=1, but outside the square, will map to
# points outside the unit circle, but are still consistent with mapping a full-sphere
# peirce projection to a full-sphere stereographic projection; however that means that
# the corners 1+j, 1-j, -1+j -1-j all map to the 'south pole' at infinity. You will get
# a divide-by-zero, or near to it, at or near those points.
# It works with array or scalar input.
#
def peirce_map( z ):
    return cplx_cn( Kval*(1-z), 0.5 )