package com.huashe.park.common.geo;

import org.apache.commons.math3.linear.*;

public class CoordinateConverter {
    // WGS84椭球体参数
    private static final double A = 6378137.0; // 长半轴（米）
    private static final double F = 1 / 298.257223563; // 扁率
    private static final double B = A * (1 - F); // 短半轴
    private static final double E = Math.sqrt(1 - (B * B) / (A * A)); // 第一偏心率

    // 将经纬度转换为平面坐标（墨卡托投影）
    public static double[] latLonToXY(double latitude, double longitude) {
        double x = A * Math.toRadians(longitude);
        double y = A * Math.log(Math.tan(Math.PI / 4 + Math.toRadians(latitude) / 2));
        return new double[]{x, y};
    }

    // 将平面坐标转换为经纬度（墨卡托投影逆变换）
    public static double[] xyToLatLon(double x, double y) {
        double lon = Math.toDegrees(x / A);
        double lat = Math.toDegrees(2 * Math.atan(Math.exp(y / A)) - Math.PI / 2);
        return new double[]{lat, lon};
    }

    // 计算基准经度和比例因子
    public static double[] calculateParameters(double[] pointA, double[] pointB) {
        double lambda0 = (pointA[0] + pointB[0]) / (2 * A);
        double scaleFactor = (pointA[1] - pointB[1]) / (A * Math.log(Math.tan(Math.PI / 4 + Math.toRadians(90) / 2) / Math.tan(Math.PI / 4 + Math.toRadians(-90) / 2)));
        return new double[]{lambda0, scaleFactor};
    }

    public static void main(String[] args) {
        // 定义已知的平面坐标和对应的经纬度
        double[][] points = {{8, 14.63, 118.86839482667, 32.0131180416999}, {6, 0, 118.869042069359, 32.0131156239015}, {0, 15.3, 118.868399686907, 32.0133391503578}};

        // 提取x, y和对应的经纬度
        double[] x = new double[points.length];
        double[] y = new double[points.length];
        double[] lon = new double[points.length];
        double[] lat = new double[points.length];

        for (int i = 0; i < points.length; i++) {
            x[i] = points[i][0];
            y[i] = points[i][1];
            lon[i] = points[i][2];
            lat[i] = points[i][3];
        }

        // 构建线性方程组Ax=B
        RealMatrix A = new Array2DRowRealMatrix(new double[][]{x, y,new double[]{1,1,1}});
        DecompositionSolver solverLon = new QRDecomposition(A).getSolver();
        DecompositionSolver solverLat = new QRDecomposition(A).getSolver();

        // 解线性方程组得到系数
        RealVector coeffsLon = solverLon.solve(new ArrayRealVector(lon));
        RealVector coeffsLat = solverLat.solve(new ArrayRealVector(lat));

        // 计算第四点的经纬度
        double xNew = 4;
        double yNew = 4;
        double lonNew = coeffsLon.getEntry(0) * xNew + coeffsLon.getEntry(1) * yNew + coeffsLon.getEntry(2);
        double latNew = coeffsLat.getEntry(0) * xNew + coeffsLat.getEntry(1) * yNew + coeffsLat.getEntry(2);

        System.out.println("第四点的经纬度为: " + lonNew + ", " + latNew);
    }
}
