天文算法--地球形状

地球相关参数计算。

package cn.ancony.chinese_calendar;

import lombok.AllArgsConstructor;
import lombok.Data;

import static java.lang.Math.*;

/**
 * 地球形状。
 * <p>
 * 大地测量学的地表不能用任何一个可定义的固体严格的描
 * 述。一种可足够地理学及天文学使用的近似假设是：把它看
 * 作回转椭球。
 * <p>
 * 计算相关地球参数：
 * <p>
 * ρ*sin(φ′)和 ρ*cos(φ′)  用于计算周日视差、日月食、星蚀所需要的量
 * <p>
 * Rp                     指定纬度φ的圆的半径
 * <p>
 * perLongitude           指定纬度φ处经度1度的长度
 * <p>
 * lineSpeed              指定纬度φ处线速度
 * <p>
 * Rm                     指定纬度φ处地球子午圈的曲率半径
 * <p>
 * perLatitude            指定纬度φ处纬度1度的长度
 * <p>
 * 地表两点的最短距离
 */
public class EarthShape {
    private static final double
            a = 6378140.0d,//地球赤道半径，单位米 (1976IAU)
            b = 6356755.0d,//地球极半径，单位米. b = a*(1-f)
            f = 1 / 298.257,//地球偏率 (1976IAU)
            e = 0.08181922d,// 地球子午圈离心率。 e = sqrt(2f-f*f)
            w = 7.292115018e-5;//地球旋转的角速度（相对于恒星，不是相对于春风点）单位：(弧度/秒)

    private static final LevelConvert lc = LevelConvert.of("度角分角秒", 60, 60);

    /**
     * 已知度分秒和高度，求 ρ*sin(φ′)和ρ*cos(φ′)
     * <p>
     * tan(u) = b / a * tan(φ)
     * ρ*sin(φ′) = b / a * sin(u) + H / 6378140 * sin(φ)
     * ρ*cos(φ′) = cos(u) + H / 6378140 * cos(φ)
     * <p>
     * ρ*sin(φ′)值在北半球为正，南半球为负，而ρ*cos(φ′)始终为正。
     * ρ表示观测者到地心的距离，地球的赤道半径看作 1 个单位
     * <p>
     *
     * @param degree 度，单位：角度
     * @param minute 分，单位：角度
     * @param second 秒，单位：角度
     * @param height 高度，米
     * @return 含两个值的double数组，第一个值表示ρ*sin(φ′)，第二个值表示ρ*cos(φ′)
     */
    public static double[] baseInfo(int degree, int minute, int second, int height) {
        double fai = toRadians(lc.toHighest(degree, minute, second));
        double ba = b / a, ha = height / a, u = atan(ba * tan(fai));
        return new double[]{ba * sin(u) + ha * sin(fai), cos(u) + ha * cos(fai)};
    }

    /**
     * 纬度 φ 的圆半径：Rp = a*cos(φ)/sqrt(1-(e*sin φ)^2)    e:子午圈椭圆的半径
     * 在同一纬度φ上，经度变化 1 度，相应的长度变化为：(π/180)Rp
     * 地球子午圈的曲率半径，在纬度φ：Rm = a(1-e^2)/((1-(e^2)*(sin φ)*(sin φ))^(3/2))
     * 地球纬度变化1度，长度变化为 (π/180)Rm
     * 在赤道时，Rm 达到最大值，值为 a(1-e^2)=6335.44km
     * 在极点时达到最小值，值为 a/sqrt(1-e^2)=6399.60km
     *
     * @param lat 纬度，单位：弧度
     * @return 由纬度可以获得的信息
     */
    public static LatitudeInfo info(double lat) {
        double sinF = sin(lat), e2 = e * e, sinF2 = sinF * sinF, base = 1 - e2 * sinF2, ra = PI / 180,
                rp = a * cos(lat) / sqrt(base), milePerLongitude = ra * rp,
                lineSpeed = w * rp, rm = a * (1 - e2) / pow(base, 1.5), milePerLatitude = ra * rm;
        return LatitudeInfo.of(rp, milePerLongitude, lineSpeed, rm, milePerLatitude);
    }

    /**
     * 纬度信息。即可以从纬度中获得的信息。
     */
    @AllArgsConstructor
    public static class LatitudeInfo {
        double
                rp,                 //Rp，单位：米
                milePerLongitude,   //经度1度的长度，单位：米
                lineSpeed,          //线速度，单位: 米/秒
                rm,                 //Rm，单位：米
                milePerLatitude;    //纬度1度的长度，单位：米

        public static LatitudeInfo of(double rp, double milePerLongitude, double lineSpeed,
                                      double rm, double milePerLatitude) {
            return new LatitudeInfo(rp, milePerLongitude, lineSpeed, rm, milePerLatitude);
        }
    }

    /**
     * 已知地球表面两点(x,y)的地理坐标，求地球表面两点间的距离。(粗略算法)
     * <p>
     * 设第 1 点的经度和纬度是分别是 L1 和 φ1，第二点的是 L2 和 φ2，我们假设这两点在海平面。
     * 如果精度要求不高，可以把地球看作球形，平均半径为 6371km。
     * <p>
     * 使用下式可得到两点间的角距离：
     * <p>
     * cos(d) = sin(φ1) * sin(φ2) + cos(φ1) * cos(φ2) * cos(L1-L2)
     * <p>
     * 最短距离 s = (6371*π*d)/180，单位：千米,d 的单位为度。
     * <p>
     * 使用条件，d不能很小。
     *
     * @param l1   x点的经度，单位：弧度
     * @param fai1 x点的纬度，单位：弧度
     * @param l2   y点的经度，单位：弧度
     * @param fai2 y点的纬度，单位：弧度
     * @return 单位：米
     */
    public static int distance(double l1, double fai1, double l2, double fai2) {
        double d = toDegrees(acos(sin(fai1) * sin(fai2) + cos(fai1) * cos(fai2) * cos(l1 - l2)));
        return (int) ((6371 * PI * d) / 180 * 1000);
    }

    /**
     * 已知地球表面两点(x,y)的地理坐标，求地球表面两点间的距离。(高精度算法,H.Andoyer)
     * <p>
     * 设赤道半径为 a，扁率为 f，则：
     * <p>
     * F = (φ1 + φ2) /2, G = (φ1 - φ2) / 2, λ = (L1 - L2) / 2
     * <p>
     * S = ((sin(G))^2) * ((cos(λ))^2) + ((cos(F))^2)*((sin(λ))^2)
     * <p>
     * C = ((cos(G))^2) * ((cos(λ))^2) + ((sin(F))^2)*((sin(λ))^2)
     * <p>
     * tan(ω) = sqrt(S / C)
     * <p>
     * R = sqrt(S * C) / ω      (ω为弧度)
     * <p>
     * D = 2 * ω * a, H1 = (3 * R-1) / (2 * C), H2 = (3 * R + 1) / (2 * S)
     * <p>
     * 最短距离 s = D * (1 + f * H1 * ((sin(F))^2) * ((cos(G))^2) -f * H2 * ((cos(F))^2) * ((sin(G))^2))
     *
     * @param l1   x点的经度，单位：弧度
     * @param fai1 x点的纬度，单位：弧度
     * @param l2   y点的经度，单位：弧度
     * @param fai2 y点的纬度，单位：弧度
     * @return 返回两点间距离。 单位：米。 该公式计算误差为 50米(天文计算)
     */
    public static int distance2(double l1, double fai1, double l2, double fai2) {
        double
                F = (fai1 + fai2) / 2, G = (fai1 - fai2) / 2, lambda = (l1 - l2) / 2,
                sinG = sin(G), cosG = cos(G), sinF = sin(F),
                cosF = cos(F), sinL = sin(lambda), cosL = cos(lambda),
                sinF2 = sinF * sinF, cosF2 = cosF * cosF, sinG2 = sinG * sinG,
                cosG2 = cosG * cosG, cosL2 = cosL * cosL, sinL2 = sinL * sinL,
                S = sinG2 * cosL2 + cosF2 * sinL2, C = cosG2 * cosL2 + sinF2 * sinL2,
                tanW = sqrt(S / C), w = atan(tanW), R = sqrt(S * C) / w,
                _3R = 3 * R, _2C = 2 * C, _2S = 2 * S,
                D = 2 * w * a, H1 = (_3R - 1) / _2C, H2 = (_3R + 1) / _2S;
        return (int) (D * (1 + f * H1 * sinF2 * cosG2 - f * H2 * cosF2 * sinG2));
    }

    @Data
    @AllArgsConstructor
    public static class Dms {
        private boolean positive;
        private int degree, minute, second;

        public static double of(int degree, int minute, int second, boolean positive) {
            double radians = of(degree, minute, second);
            return positive ? radians : -radians;
        }

        //默认弧度为正值
        public static double of(int degree, int minute, int second) {
            return toRadians(lc.toHighest(degree, minute, second));
        }
    }

}