Java通过两点经纬度和范围距离计算范围四个顶点的经纬度

   参考博文：https://blog.csdn.net/zengmingen/article/details/68490497

private static final double mile = 1609.344;//米-》英里换算
//通过起点、终点的经纬度求方位角
double azimuth = gps2d(originLat, originLong, destinationLat, destinationLong);
double range = 50;//范围 单位：英里mile
double dist = (range+1)*mile;//由于计算顶点经纬度的方法存在误差，所以距离需要+1
//求范围四个顶点的经纬度
Map<String,Object> origin1 = computerThatLonLat(originLong, originLat, azimuth+90, dist);
Map<String,Object> origin2 = computerThatLonLat(originLong, originLat, azimuth-90, dist);
Map<String,Object> destination1 = computerThatLonLat(destinationLong, destinationLat, azimuth+90, dist);
Map<String,Object> destination2 = computerThatLonLat(destinationLong, destinationLat, azimuth-90, dist);

    /* 
     * 大地坐标系资料WGS-84 长半径a=6378137 短半径b=6356752.3142 扁率f=1/298.2572236 
     */  
    /** 长半径a=6378137 */  
    private static double a = 6378137;  
    /** 短半径b=6356752.3142 */  
    private static double b = 6356752.3142;  
    /** 扁率f=1/298.2572236 */  
    private static double f = 1 / 298.2572236;  
  
    /** 
     * 经纬度换成弧度 （°）  
     * @param d 度 
     * @return 弧度 
     */  
    private static double rad(double d) {  
        return d * Math.PI / 180.0;  
    } 
  
    /** 
     * 弧度换成度  
     * @param x 弧度 
     * @return 度 
     */  
    private static double deg(double x) {  
        return x * 180 / Math.PI;  
    } 
    /** 
     * 通过一个点的经纬度、距离、方位角，计算另一点经纬度
     * @param lon 经度 
     * @param lat 纬度 
     * @param lonlat 已知点经纬度 
     * @param brng 方位角 
     * @param dist 距离（米） 
     */  
    public static Map<String,Object> computerThatLonLat(double lon, double lat, double brng, double dist) {
    	Map<String,Object> map = new HashMap<String,Object>();
    	double alpha1 = rad(brng);  
        double sinAlpha1 = Math.sin(alpha1);  
        double cosAlpha1 = Math.cos(alpha1);    
        double tanU1 = (1 - f) * Math.tan(rad(lat));  
        double cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1));  
        double sinU1 = tanU1 * cosU1;  
        double sigma1 = Math.atan2(tanU1, cosAlpha1);  
        double sinAlpha = cosU1 * sinAlpha1;  
        double cosSqAlpha = 1 - sinAlpha * sinAlpha;  
        double uSq = cosSqAlpha * (a * a - b * b) / (b * b);  
        double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));  
        double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));    
        double cos2SigmaM=0;  
        double sinSigma=0;  
        double cosSigma=0;  
        double sigma = dist / (b * A), sigmaP = 2 * Math.PI;  
        while (Math.abs(sigma - sigmaP) > 1e-12) {  
            cos2SigmaM = Math.cos(2 * sigma1 + sigma);  
            sinSigma = Math.sin(sigma);  
            cosSigma = Math.cos(sigma);  
            double deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)  
                    - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));  
            sigmaP = sigma;  
            sigma = dist / (b * A) + deltaSigma;  
        }  
  
        double tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;  
        double lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,  
                (1 - f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));  
        double lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);  
        double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));  
        double L = lambda - (1 - C) * f * sinAlpha  
                * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));  
  
        double revAz = Math.atan2(sinAlpha, -tmp); // final bearing  
        double new_lng = lon+deg(L);
        double new_lat = deg(lat2);
        map.put("longitude", new_lng);
        map.put("latitude", new_lat);
        return map;
    }    
    
    /**
     * 经纬度方位角计算
     * @param lat_a 起点纬度
     * @param lng_a 起点经度
     * @param lat_b 终点纬度
     * @param lng_b 终点经度
     * @return 方位角度数
     */
    public static double gps2d(double lat_a, double lng_a, double lat_b, double lng_b) {
        double d = 0;
        lat_a=lat_a*Math.PI/180;
        lng_a=lng_a*Math.PI/180;
        lat_b=lat_b*Math.PI/180;
        lng_b=lng_b*Math.PI/180;     
        d=Math.sin(lat_a)*Math.sin(lat_b)+Math.cos(lat_a)*Math.cos(lat_b)*Math.cos(lng_b-lng_a);
        d=Math.sqrt(1-d*d);
        d=Math.cos(lat_b)*Math.sin(lng_b-lng_a)/d;
        d=Math.asin(d)*180/Math.PI;
        return d;
     }