Crossing counting in `phicurve` failing in an edge case

[williameclee]

The issue

The phicurve function counts the number of crossings with the number of sign changes (line 40). Therefore, when a vertex sits on the colatitude at where crossings are evaluated, the number of sign changes erroneously counts one additional sign changes, resulting in non-even number of crossings (the error in line 108).

Reproducing

To reproduce:

X = [-15, -15, -15, 15, 15, -15];
Y = [15, 0, -15, -15, 15, 15];
XY = [X', Y'];
collon = [90 - XY(:, 2), XY(:, 1)];
th = 70:10:110;
[phint, thp, php] = phicurve(collon, th);

In the case above, XY defines an ordinary square, but the second vertex has the colatitude of 90, which is one of the values in th, and the function cannot correctly count the number of crossings.

Why this matters

This problem is usually not relevant for landmasses, since the Gauss-Legendre points (the inputting th array) are almost always non-rational and do not share colatitude with vertex points. But when the geometry has northern and southern bound symmetrical to the equator, like oceans with symmetric polar caps (or the square above), one of the Gauss-Legendre points can have the colatitude of 90. The function therefore fails when a vertex also have the colatitude of 0 (which is likely for manually drawn geometries).

Some possible workarounds

I don't think it is worth it to come up with a new crossing-counting algorithm, but there are some simple workarounds to prevent this problem. A simple solution is to make sure no vertices share the same colatitude with any Gauss-Legendre points, which are usually irrational. Therefore by manually offsetting the colatitudes by a very small amount, e.g. by changing line 39 to

xmth = repmat(thph(:, 1), 1, length(th)) + (1e-14) - repmat(th(:)', length(thph(:, 1)), 1);

the problem can be bypassed. This is not elegant, but it works. A slightly more sophisticated solution is to only offset the vertices causing the problem. Another possible solution is to simply remove such vertices, if there will be enough vertices remaining.

[williameclee]

The more sophisticated solution mentioned above can be implemented like:

xmth = repmat(thph(:, 1) + any(thph(:,1) == th,2) * 1e-14, 1, length(th)) - repmat(th(:)', length(thph(:, 1)), 1);

[williameclee]

To avoid potentially altering the behaviour of other functions that call phicurve, it is also possible to only recompute xmth after the counting fails (i.e. in the if loop at line 50).