Better way to handle poles?

[HansBrende]

Calculating the distance over a pole, e.g., between [0, 89.99] and [180, 89.99] (with 89.99 as the latitude passed into the cheap ruler constructor) gives me a 57% error relative to the Vincenty calculation. (For reference, I used https://github.com/chrisveness/geodesy which seems to be better maintained than the derivative https://github.com/TankofVines/node-vincenty)

I know that this is already documented as a shortcoming, but food for thought: if you can calculate the maximum error that your shortcut formula would give you in advance (which seems likely, given that you have a % error chart in the readme), could you not simply fall back to the Vincenty (or even haversine) formula if the error is too large? That would fix the pole issue... (for comparison, the haversine formula gives only a 0.4% error for those same two coordinates -- and could probably be improved upon by tweaking the earth radius).

[kkdd]

The Flat-surface formula, used by cheap-ruler, has singularity at the poles.

Use other non-singular formulas instead, such as Gauss mid-latitude method.

The following shows the accuracy test including Gauss mid-latitude method:

[Artoria2e5]

Hmm. We already basically have the formulae for the M and N parts of Gauss. A very basic translation would be

class Ruler {
    constructor(lat, units) {
        const r = RE * (units ? factors[units] : 1);
        const coslat = Math.cos(lat * RAD);
        const w2 = 1 / (1 - E2 * (1 - coslat * coslat));
        const w = Math.sqrt(w2);

        // multipliers for converting longitude and latitude degrees into distance
        this.norm = r * w;                 // normal radius of curvature
        this.meri = r * w * w2 * (1 - E2); // meridonal radius of curvature
        this.coslat = coslat;
    }

    distance([lon1, lon2], [lat1, lat2]) {
        const dx = (lon2 - lon1) * RAD;
        const dy = (lat2 - lat1) * RAD;

        return 2 * this.norm * Math.asin(
            Math.sqrt((Math.sin(dx / 2) * this.coslat) ** 2 +
            (Math.cos(dx / 2) * Math.sin(this.meri * dy / 2 / this.norm)) ** 2);
    }
}

(Maybe this is wrong! I'm not great at Greek algebra.)

Correct or not, this translation looks a lot more expensive than the planar one -- I count one more asin, one more sincos (well, some languages allow calculating them at the same time -- not js), and one more sin. Can we cut some corners somewhere?

(Well, the story is that I'd like to have something that's cheaper than haversine -- maybe even ol' unreliable cosine -- that's accurate to 10%. Yeah, pretty big margin there.)

[kkdd]

Hello, The quite-shorter-range approximation of Gauss mid-latitude method is as follows:

const londiffhalf = (lon1 - lon2) / 2 * degree;
const latdiffhalf = (lat1 - lat2) / 2 * degree;
const lat = (lat1 + lat2) / 2 * degree; // middle
const coslat = cos(lat);
const n2 = 1 / (1 - e2 * (1 - coslat ** 2));
const n = sqrt(n2);  // prime vertical radius of curvature, or normal
const m_by_n = (1 - e2) * n2;  // meridian ratio to normal

const c = tan(londiffhalf / 2);  // tangent half-angle formula
const sin_londiffhalf = 2 * c / (1 + c ** 2);
const cos_londiffhalf = (1 - c ** 2) / (1 + c ** 2);

return 2 * n * RE * hypot(
  sin_londiffhalf * coslat,
  cos_londiffhalf * latdiffhalf * m_by_n);

The above is derived by applying $\arcsin \Delta \approx \Delta$ and $\sin (\frac{\Delta \phi}{2} \frac{M}{N}) \approx \frac{\Delta \phi}{2} \frac{M}{N}$ into the following term:

asin(hypot(
  sin(londiffhalf) * coslat,
  cos(londiffhalf) * sin(latdiffhalf * m_by_n))
))

Note that https://en.wikipedia.org/wiki/Tangent_half-angle_formula is usefull, additionally.