Crossing counting in phicurve failing in an edge case #22
Description
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.