このホームページ「GPSwalking」の地図ページの撮影ポイント間の距離を表示したいと思っていたが計算方法が分からなかった。ネットでヒュベニの公式というものを使うと簡単にできることを偶然見つけた。
https://hp.vector.co.jp/authors/VA002244/yacht/geo.htm

精度は50km以下ではLambert-Andoyerの公式（航海算法）の 10^-5 以下だと書いてある。今回の利用法はGPSで取得した２点間の距離だから、歩きの場合は数メートル、車で移動の場合でも数十メートルぐらいだから、精度は心配なさそうだ。

いろいろややこしい説明図があるが、このヒュベニの公式を理解する頭はなくても公式を使わせてもらって、計算することは出来そうな気がする。
上記のサイトの式を分かり易く書き直すと

GPXデータのポイントの緯度経度は「度」だから、ラジアンに変換しておく必要がある。
長半径（赤道半径）と短半径（極半径）はWikipediaからWGS84のデータを使うことにした。
ここで1番目は日本測地系、2番目は世界測地系、3番目はGPSで使われているWGS84測地系。
2番目と3番目の違いは、僅か0.1mm。この地球の大きさをそこまで精密に調べてあることに驚きだ。

第2離心率e2、およびMの分子a1e2は、長半径・短半径が分かったので、前もって計算できるので、これらを元に関数にまとめるたのが下記。

function hubeny(lat1, lng1, lat2, lng2) {
    function rad(deg) {
        return deg * Math.PI / 180;		//degreeをラジアンに変換
    }
    lat1 = rad(lat1);				//degreeをラジアンに変換した値
    lng1 = rad(lng1);
    lat2 = rad(lat2);
    lng2 = rad(lng2);
    var latDiff = lat1 - lat2;			// Ay：緯度差
    var lngDiff = lng1 - lng2;			// Ax：経度差
    var latAvg = (lat1 + lat2) / 2.0;		// P：緯度の平均

    var a = 6378137.0;				// Rx：赤道半径（WGS84系）
    var b = 6356752.314245;			// Ry：極半径（WGS84系）
    var e2 = 0.00669437999019758;		// Eの２乗：第2離心率　e^2 = (a^2 - b^2) / a^2
    var a1e2 = 6335439.32729246;		// Mの分子：赤道上の子午線曲率半径  a(1 - e^2)

    var sinLat = Math.sin(latAvg);
    var W2 = 1.0 - e2 * (sinLat * sinLat);
    var M = a1e2 / (Math.sqrt(W2) * W2);	// 子午線曲率半径M（分母はW^3）
    var N = a / Math.sqrt(W2);			// 卯酉線（ぼうゆうせん）曲率半径N

    t1 = M * latDiff;
    t2 = N * Math.cos(latAvg) * lngDiff;
    return Math.sqrt((t1 * t1) + (t2 * t2));
}

参考にしたサイト
・ヒュベニの式を用いた、緯度・経度と距離・方位の相互変換の解説
・TRAIL NOTE