package com.tourism.common.core.util;

public class MathCalculateUtils {
    public static double calculateDistance(double lat1, double lon1, double lat2, double lon2) {
        double a = 6378137.0; // 第一偏心率
        double b = 6356752.3142; // 第二偏心率
        double f = 1 / 298.257223563; // 扁率

        double L = Math.toRadians(lon2 - lon1);
        double U1 = Math.atan((1 - f) * Math.tan(Math.toRadians(lat1)));
        double U2 = Math.atan((1 - f) * Math.tan(Math.toRadians(lat2)));
        double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
        double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);

        double lambda = L;
        double lambdaP = 2 * Math.PI;
        int iterLimit = 100;

        double sinSigma = 0, cosSigma = 0, sigma = 0, sinAlpha = 0;
        double cosSqAlpha = 0, cos2SigmaM = 0;

        while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0) {
            double sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);

            sinSigma = Math.sqrt((Math.pow(cosU2 * sinLambda, 2)) +
                    (Math.pow(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda, 2)));
            cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
            sigma = Math.atan2(sinSigma, cosSigma);
            sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
            cosSqAlpha = 1 - Math.pow(sinAlpha, 2);
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;

            double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
            lambdaP = lambda;
            lambda = L + (1 - C) * f * sinAlpha *
                    (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * Math.pow(cos2SigmaM, 2))));
        }

        double uSq = cosSqAlpha * (Math.pow(a, 2) - Math.pow(b, 2)) / Math.pow(b, 2);
        double A = 1 + uSq / 16384 *
                (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
        double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));

        double deltaSigma = B * sinSigma *
                (cos2SigmaM + B / 4 *
                        (cosSigma * (-1 + 2 * Math.pow(cos2SigmaM, 2)) -
                                B / 6 * cos2SigmaM *
                                        (-3 + 4 * Math.pow(sinSigma, 2)) *
                                        (-3 + 4 * Math.pow(cos2SigmaM, 2))));
        double distance = b * A * (sigma - deltaSigma);

        return distance;
    }
}
