天文算法--章动和黄赤交角

章动和黄赤交角的计算方法

package cn.ancony.chinese_calendar;

import java.util.Arrays;

import static java.lang.Math.*;

/**
 * 章动和黄赤交角
 * <p>
 * 地球章动：地球自转轴围绕其平位置周期振动。
 * <p>
 * 影响：地球自转轴围绕黄极做岁差运动的同时，自转轴还会在岁差运动的轨迹上波动。
 * <p>
 * 章动主要是月球运动引起的，可以描述为一些周期项的和。主要项的周期是 6798.4 日(18.6 年)，
 * 但其它项是一些短周期项(小于 10 天)。
 * <p>
 * 章动可以很容易的分解为黄道的水平分量和的垂直分量。
 * <p>
 * 黄道上的水平分量记为Δψ，称为黄经章动；它影响了天球上所有天体的经度。
 * 黄道的垂直分量记为Δε，称为交角章动，它影响了黄赤交角。
 * <p>
 * 计算所有天体的“视位置”及“恒星时”，都需要计算章动。
 * <p>
 * 黄赤交角：是地球自转轴的倾角，它也是黄道面与赤道面的夹角。
 * <p>
 * 有平黄赤交角(黄道与平赤道的夹角，用εo表示)与真黄赤交角(黄道与真赤道的夹角，用ε表示)之分。
 * <p>
 * 真黄赤交角是ε=εo+Δε，Δε是交角章动
 */
public class NutationAndObliquity {
    /**
     * 以下单位为度： IAU1980 章动理论(忽略了系数小于 0".0003 的项)。系数的单位是 0".0001
     * <p>
     * 平距角(日月对地心的角距离)：D = 297.85036 +455267.111480*T - 0.0019142*T^2 + T^3/189474
     * <p>
     * 太阳（地球）平近点角：M = 357.52772 + 35999.050340*T - 0.0001603*T^2 - T^3/300000
     * <p>
     * 月球平近点角：M′= 134.96298 + 477198.867398*T + 0.0086972*T^2 + T^3/56250
     * <p>
     * 月球纬度参数：F = 93.27191 + 483202.017538*T - 0.0036825*T^2 + T^3/327270
     * <p>
     * 黄道与月球平轨道升交点黄经，从 Date 黄道平分点开始测量：Ω= 125.04452 - 1934.136261*T + 0.0020708*T^2 + T^3/450000
     * <p>
     * L = 280°.4665 + 36000.7698*T  月球太阳的平黄经
     * <p>
     * L′= 218°.3165 + 481267°.8813*T 太阳的平黄经
     * <p>
     * Δψ精度: 0".5，Δε精度: 0".1 (忽略以上表达式中的 T 平方项及 T 三次方项),则：
     * <p>
     * Δψ=-17".20sin(Ω)-1".32sin(2L)-0".23sin(2L′)+0".21sin(2Ω)
     * <p>
     * Δε=+9".20cos(Ω)+0".57cos(2L)+0".10cos(2L′)-0".09cos(2Ω)
     *
     * @param t J2000.0 起算的儒略世纪数：即 (JDE -2451545)/36525  JDE 是历书儒略日数
     * @return [Δψ, Δε],即[黄经章动,交角章动]。单位：秒
     */
    public static double[] nutation(double t) {
        double t2 = t * t, t3 = t2 * t;
        double d = 297.85036 + 455267.111480 * t - 0.0019142 * t2 + t3 / 189474,
                m = 357.52772 + 35999.050340 * t - 0.0001603 * t2 - t3 / 300000,
                m2 = 134.96298 + 477198.867398 * t + 0.0086972 * t2 + t3 / 56250,
                f = 93.27191 + 483202.017538 * t - 0.0036825 * t2 + t3 / 327270,
                omega = 125.04452 - 1934.136261 * t + 0.0020708 * t2 + t3 / 450000,
                l = 280.4665 + 36000.7698 * t,
                l2 = 218.3165 + 481267.8813 * t;
        double rOmega = toRadians(omega), r2l = toRadians(2 * l),
                r2l2 = toRadians(2 * l2), r2omega = toRadians(2 * omega);
        double ans1 = -17.20 * sin(rOmega) - 1.32 * sin(r2l) - 0.23 * sin(r2l2) + 0.21 * sin(r2omega),
                ans2 = 9.20 * cos(rOmega) + 0.57 * cos(r2l) + 0.10 * cos(r2l2) - 0.09 * cos(r2omega);
        return new double[]{ans1, ans2};
    }

    private static final double start = 84381.448;//23°26′21".448 转换为秒

    /**
     * 平黄赤交角的计算方法一(IAU 提供)：
     * 2000年为εo 误差 1"，4000 年εo误差为 10"。
     * εo = 23°26'21".448 - 46".8150*T - 0".00059*T^2 + 0".001813*T^3
     *
     * @param t J2000.0 起算的儒略世纪数
     * @return 黄赤交角，单位：秒
     */
    public static double epsilonO1(double t) {
        double t2 = t * t, t3 = t2 * t;
        return start - 46.8150 * t - 0.00059 * t2 + 0.001813 * t3;
    }

    /**
     * 平黄赤交角的计算方法二(改良的的表达式，比方法一精度更高，Laskar提供)：
     * 精度是：1000 年后误差 0".01(公元 1000 到 3000)。10000 年后误差数个角秒。
     * 适用于|U|<1，即 J2000.0 起算前后各 10000 年的范围内。
     *
     * @param u J2000.0 起算的儒略万年数，或U=T/100
     * @return 黄赤交角, 单位：秒
     */
    public static double epsilonO2(double u) {
        double u2 = u * u, u3 = u2 * u, u4 = u3 * u, u5 = u4 * u, u6 = u5 * u,
                u7 = u6 * u, u8 = u7 * u, u9 = u8 * u, u10 = u9 * u;
        return start
                - u * 4680.93 - u2 * 1.55 + u3 * 1999.25 - u4 * 51.38 - u5 * 249.67
                - u6 * 39.05 + u7 * 7.12 + u8 * 27.87 + u9 * 5.79 + u10 * 2.45;
    }

}