Confused about your implementation of Lee's solution to Dixon's elliptic functions

[carbonox-infernox]

Recently I've become interested in Lee's conformal world in a tetrahedron projection. I've also been wanting to start getting experience with java, so I was drawn to this repo. I noticed that you wrote the wiki article for this projection, so you seem to be the right person to talk to about this.

For several reasons, I can't quite wrap my head around your implementation of Lee's approximation to Dixon's elliptic functions. I was hoping you would help me understand some things.

I've reproduced this projection with MapDesignerRaster.jar and obviously your implementation succeeds. My confusion probably stems from unfamiliarity with Java, or some math misunderstanding.

Things I'm confused about:

coefficients are all positive (not alternating sign)

In Lee's paper, the coefficients to the MacLaurin series alternate in sign like this:

+0.625000
-0.223214
+0.069754
-0.020121
...

but in your implementation, they're all positive. I'm not seeing how negatives could be alternatingly introduced in this line of code:

sum = sum.times(w3).plus(COEF[i]);

w³ taken initially and not updated after each iteration of converging series

In the following block of code:

public static Complex leeFunc(Complex w) {
    Complex w3 = w.pow(3);
    Complex sum = new Complex(0);
    for (int i = COEF.length-1; i >= 0; i --)
        sum = sum.times(w3).plus(COEF[i]);
    return sum.times(w);
}

w³ is computed before the iterative process, but then never updated. Although Lee didn't say this explicitly, he seemed to heavily imply that every iteration should add the coefficient times the current best estimate of w³

Iterative process starts with least significant term and goes backwards.

In Lee's paper, higher order coefficients (further down on the list) accumulate more powers of w, but your process (above code, for loop) seems to start at the end of the coefficients and then tack on powers of w for lower order (higher up on the list) terms. What I mean is that yours will look like this (I'm ignoring powers of 10^-1 for now):

+0.625000 w¹³ -0.223214 w¹⁰ +0.069754 w⁷ -0.020121 w⁴

when it should look like this:

+0.625000 w⁴ -0.223214 w⁷ +0.069754 w¹⁰ -0.020121 w¹³

Lee says to multiple by w³ * 10^-1 in each iteration, so the powers of 10^-1 are taken care of. I wrote an implementation in python that should be correct:

import numpy as np
COEF = [
    .625000,
   -.223214,
    .069754,
   -.020121,
    .005539,
   -.001477,
    .000385,
   -.000099,
    .000025
    ]
def leeFunc(alpha, sigma):
    w = 1.791797 * np.tan(sigma/2) * np.exp(alpha * 1j)
    for coef in COEF:
        print(w)
        w = w + coef * w**3 * 1e-1
    print(w)
# or alternatively, take just w as input, but this is where the initial w is supposed to come from

This converges for me for values of w within one spherical triangle.

sigma and alpha description on wikipedia does not match Lee

The Wikipedia article says that α is the longitude and σ is the angular distance from the pole, whereas Lee defines α as the azimuth and σ as the arc length "of a great circle from the centroid to any other point within the spherical triangle." It seems that these only align when a pole is the centroid and a vertex lies on the prime meridian. In Lee's projection, only one triangle has a pole (south) as a centroid, and the vertices are at 20W, 140W, and 100E.

I haven't dug deep enough into this repo to see how you're getting the initial w, but you must be doing a transformation somewhere.

[jkunimune]

Glad to help!

The alternating sign thing isn't terribly important, because the powers of w also alternate between even and odd. So if my w is defined as Lee's w times negative one (I can't remember why that would be, but it seems to be so), then all of the coefficients in front of an odd power of w change sign, and all of the coefficients are positive. i.e. if w3 is negative, then there are alternating negatives. Multiplying by negative one is just a 180° rotation in the complex plane, so it doesn't make a substantial difference as long as I made sure my sign conventions were consistent in my code.

The leeFunc(w) function you're looking at isn't the iterative process. That's just evaluating the polynomial. The iterative part is the next function down, invFunc(z). You can see that that one takes z as the initial guess of w and then repeatedly calls leeFunc(w) and uses that to update the value of w.

Note that i starts at the end of the coefficient array and decrements down to 0. So the coefficients do match the powers of w, just in the opposite order as Lee has them.

+0.020121* w^13
-0.069754 * w^10
+0.223214 * w^7
-0.625000 * w^4

The reason I do it in that order is just that ((a x + b) x + c) x + d is more efficient to compute than d + c x + b x^2 + a x^3.

You're right that those definitions match up when the centroid is also the pole. Lee considers his projection as a combination of four projections onto four triangles, but by the periodic nature of the Dixon functions, it can also be considered as a single projection on a single large triangle. If one simply extends the borders of the center triangle to the rest of the sphere, the other three triangles fall into place (though convergence isn't as good, so for computational purposes, I also implemented it as four separate projections). So when I wrote that Wikipedia page, I just referred to it as a single equation on a single face centered on a pole to simplify the explanation somewhat; it's equivalent.

The transformation is in src/maps/Tetrahedral.java. That's where I define the arrangement (Configuration.TETRAHEDRON_WIDE_FACE and Configuration.TRIANGLE_FACE), define my initial w for the sphere-to-plane projection (lines 48\~49), and extract latitude and longitude from the polynomial value for the plane-to-sphere projection (lines 58\~60). Note that I just use "lat" and "lon" instead of "distance from centroid" and "azimuth", because in the general PolyhedralProjection class, I already rotate the centroid onto the pole before passing those transformed values to faceProject(lat, lon) (line 348).

I hope that helps!