2点の緯度経度の距離を計算する

#経緯
専門学校時代にマップを扱うことがありその時に作った関数をこの機会に残しておこうと思い書きました。

#解説
山だらけにて解説されているヒュベニの公式というものを使いました。式は

\sqrt{(d_yM)^2+(d_zN \cos µ_y)^2}

です。

わかりやすく言葉に直すと

\sqrt{(緯度の差 × 子午線曲線半径)^2 + (経度の差 × 卯酉線曲率半径 × cos緯度の平均値)^2}

何をしてるんだろうって私も思いました。
よく見ると

\sqrt{(c−a)^2+(d−b)^2}

という二点間の距離を求める公式になります。
そこに地球を楕円と仮定した場合の誤差を補正した式らしいです。

#コード

HubenyDistance.java

public class HubenyDistance {

    //世界観測値系
    public static final double GRS80_A = 6378137.000;//長半径 a(m)
    public static final double GRS80_E2 = 0.00669438002301188;//第一遠心率  eの2乗

    public static double deg2rad(double deg){
        return deg * Math.PI / 180.0;
    }

    public static double calcDistance(double lat1, double lng1, double lat2, double lng2){
        double my = deg2rad((lat1 + lat2) / 2.0); //緯度の平均値
        double dy = deg2rad(lat1 - lat2); //緯度の差
        double dx = deg2rad(lng1 - lng2); //経度の差

        //卯酉線曲率半径を求める(東と西を結ぶ線の半径)
        double sinMy = Math.sin(my);
        double w = Math.sqrt(1.0 - GRS80_E2 * sinMy * sinMy);
        double n = GRS80_A / w;

        //子午線曲線半径を求める(北と南を結ぶ線の半径)
        double mnum = GRS80_A * (1 - GRS80_E2);
        double m = mnum / (w * w * w);

        //ヒュベニの公式
        double dym = dy * m;
        double dxncos = dx * n * Math.cos(my);
        return Math.sqrt(dym * dym + dxncos * dxncos);
    }
}