天文算法--分点和至点

本文给出二分点和二至点的大约时间的算法。参考《天文算法》。适用年份为-1000年~+3000年，精度在代码中给出。后续再补全高精度算法。

package cn.ancony.chinese_calendar;

import lombok.AllArgsConstructor;
import lombok.Data;

import static java.lang.Math.*;

/**
 * 分点和至点
 * <p>
 * 分点和至点时刻是指：太阳的地心视黄经(含光行差和章动)为 90 的整位数时对应的时刻。
 * <p>
 * 大约时间的计算方法：(适用范围：-1000年~+3000年)
 * 1. 使用代码中jde0开头的方法计算对应的平分点或平至点，得到JDEo。
 * 2. 计算T，W,
 * T = (JDEo - 2451545.0)/36525
 * W = 35999°.373T - 2°.47
 * Δλ= 1 + 0.0334 cos W + 0.0007 cos 2W
 * 3. 使用T和sTable计算S。
 * S = Σ(A * cos(B + C * T))
 * 4. 要计算的分点或至点时刻(儒略历书时,即力学时）表达为：
 * JDE = JDEo + 0.00001*S/Δλ
 * <p>
 * 这种方法，在公元1951——2050的精度如下：
 * |-------|------------------|------------------|------|
 * |       |误差小于20秒的年份个数|误差小于40秒的年份个数|最大误差|
 * |-------|------------------|------------------|------|
 * | 3月分点|       76         |        97        |  51  |
 * | 6月分点|       80         |       100        |  39  |
 * | 9月分点|       78         |        99        |  44  |
 * |12月分点|       68         |        99        |  41  |
 * |-------|------------------|------------------|------|
 */
public class EquinoxAndSolstice {

    //计算大致时间的公式的使用范围。单位：年
    private static final int low = -1000, mid = 1000, up = 3000;

    /**
     * 计算某一个年的分点和至点
     *
     * @param year 年份。范围：-1000 ~ +3000。
     * @return 返回四个分点/至点。分别为3月分点，6月至点，9月分点，12月至点。单位：JDE
     */
    public static double[] jde(int year) {
        checkYear(year);
        return year > mid ?
                new double[]{jde(jde010(year)), jde(jde011(year)), jde(jde012(year)), jde(jde013(year))} :
                new double[]{jde(jde000(year)), jde(jde001(year)), jde(jde002(year)), jde(jde003(year))};
    }

    //返回春分点
    public static double springEquinoxJde(int year) {
        checkYear(year);
        return year > mid ? jde(jde010(year)) : jde(jde000(year));
    }

    //返回夏至点
    public static double summerSolsticeJde(int year) {
        checkYear(year);
        return year > mid ? jde(jde011(year)) : jde(jde001(year));
    }

    //返回秋分点
    public static double autumnalEquinoxJde(int year) {
        checkYear(year);
        return year > mid ? jde(jde012(year)) : jde(jde002(year));
    }

    //返回冬至点
    public static double winterSolsticeJde(int year) {
        checkYear(year);
        return year > mid ? jde(jde013(year)) : jde(jde003(year));
    }

    private static void checkYear(int year) {
        if (year < low || year > up)
            throw new RuntimeException(String.format("计算范围是：%s ~ %s", low, up));
    }

    /**
     * 传入JDEo，得到分点或至点时刻(JDE)
     *
     * @param jde0 JDEo,分或至的时刻
     * @return JDE(儒略历书时, 即力学时)
     */
    private static double jde(double jde0) {
        double t = (jde0 - 2451545.0) / 36525, w = 35999.373 * t - 2.47,
                deltaLambda = 1 + 0.0334 * cos(toRadians(w)) + 0.0007 * cos(toRadians(2 * w));
        //计算S
        double s = 0;
        for (int i = 0; i < sTable.length; i += 3) {
            double a = sTable[i], b = sTable[i + 1], c = sTable[i + 2];
            s += a * cos(toRadians(b + c * t));
        }
        return jde0 + 0.00001 * s / deltaLambda;
    }

    // ==================== -1000~+1000 JDE0的计算 ====================
    private static double year0(int year) {
        return (double) year / 1000;
    }

    //-1000~+1000 3月平分点时刻JDEo
    private static double jde000(int year) {
        return c0(year0(year), new byte[]{1, 1, 1, 1, 0},
                new double[]{1721139.29189, 365242.13740, 0.06134, 0.00111, 0.00071});
    }

    //-1000~+1000 6月平至点时刻JDEo
    private static double jde001(int year) {
        return c0(year0(year), new byte[]{1, 1, 0, 1, 1},
                new double[]{1721233.25401, 365241.72562, 0.05323, 0.00907, 0.00025});
    }

    //-1000~+1000 9月平分点时刻JDEo
    private static double jde002(int year) {
        return c0(year0(year), new byte[]{1, 1, 0, 0, 1},
                new double[]{1721325.70455, 365242.49558, 0.11677, 0.00297, 0.00074});
    }

    //-1000~+1000 12月平至点时刻JDEo
    private static double jde003(int year) {
        return c0(year0(year), new byte[]{1, 1, 0, 0, 0},
                new double[]{1721414.39987, 365242.88257, 0.00769, 0.00933, 0.00006});
    }

    // ==================== +1000~+3000 JDEo的计算 ====================
    private static double year1(int year) {
        return (double) (year - 2000) / 1000;
    }

    //+1000~+3000 3月平分点时刻JDEo
    private static double jde010(int year) {
        return c0(year1(year), new byte[]{1, 1, 1, 0, 0},
                new double[]{2451623.80984, 365242.37404, 0.05169, 0.00411, 0.00057});
    }

    //+1000~+3000 6月平至点时刻JDEo
    private static double jde011(int year) {
        return c0(year1(year), new byte[]{1, 1, 1, 1, 0},
                new double[]{2451716.56767, 365241.62603, 0.00325, 0.00888, 0.00030});
    }

    //+1000~+3000 9月平分点时刻JDEo
    private static double jde012(int year) {
        return c0(year1(year), new byte[]{1, 1, 0, 1, 1},
                new double[]{2451810.21715, 365242.01767, 0.11575, 0.00337, 0.00078});
    }

    //+1000~+3000 12月平至点时刻JDEo
    private static double jde013(int year) {
        return c0(year1(year), new byte[]{1, 1, 0, 0, 1},
                new double[]{2451900.05952, 365242.74049, 0.06223, 0.00823, 0.00032});
    }

    /**
     * 计算形如 ratio1 * y^0 + ratio2 * y^1 + ratio3 * y^2 + ... + ratioN * y^(N-1)方程的结果
     *
     * @param y        y
     * @param positive 是否为正号。只能取1(正号)和0(负号).
     * @param ratios   系数
     * @return 结果
     */
    private static double c0(double y, byte[] positive, double[] ratios) {
        int len1 = ratios.length, len2 = positive.length;
        if (len1 != len2) throw new RuntimeException("数组大小必须等长");
        double[] yes = new double[len1];
        yes[0] = 1;
        double ans = ratios[0];
        for (int i = 1; i < len1; i++) {
            yes[i] = y * yes[i - 1];
            double temp = ratios[i] * yes[i];
            if (positive[i] == 1) ans += temp;
            else ans -= temp;
        }
        return ans;
    }

    /**
     * 3个一组，分别表示A,B,C。共24组。用来计算S。B和C的单位是"度"。(对应天文算法中的表26.C)
     * S = Σ(A * cos(B + C * T))
     */
    public static final double[] sTable = {
            485, 324.96, 1934.136, 203, 337.23, 32964.467, 199, 342.08, 20.186, 182, 27.85, 445267.112,
            156, 73.14, 45036.886, 136, 171.52, 22518.443, 77, 222.54, 65928.934, 74, 296.72, 3034.906,
            70, 243.58, 9037.513, 58, 119.81, 33718.147, 52, 297.17, 150.678, 50, 21.02, 2281.226,
            45, 247.54, 29929.562, 44, 325.15, 31555.956, 29, 60.93, 4443.417, 18, 155.12, 67555.328,
            17, 288.79, 4562.452, 16, 198.04, 62894.029, 14, 199.76, 31436.921, 12, 95.39, 14577.848,
            12, 287.11, 31931.756, 12, 320.81, 34777.259, 9, 227.73, 1222.114, 8, 15.45, 16859.074};

}

以下是测试类：

package cn.ancony.javafx.chinese_calendar;

import cn.ancony.chinese_calendar.JulianDay;
import cn.ancony.chinese_calendar.LevelConvert;
import cn.ancony.chinese_calendar.YearMonthDay;
import org.junit.jupiter.params.ParameterizedTest;
import org.junit.jupiter.params.provider.CsvSource;

import static cn.ancony.chinese_calendar.EquinoxAndSolstice.*;
import static org.junit.jupiter.api.Assertions.assertTrue;

public class TestEquinoxAndSolstice {
    /**
     * 测试数据为1991-2000年二分二至时刻，VSOP87理论计算，力学时。精度：接近于秒。
     * 参数含义依次为：年份，3月分点，6月至点，9月分点，12月分点。
     * 分点和至点均由5个参数组成，前两个表示月份和日期，后三个表示时分秒。
     */
    @ParameterizedTest
    @CsvSource({
            "1991,  3,21,  3, 2,54,  6,21, 21,19,46,  9,23, 12,49, 4,  12,22,  8,54,38",
            "1992,  3,20,  8,49, 2,  6,21,  3,15, 8,  9,22, 18,43,46,  12,21, 14,44,14",
            "1993,  3,20, 14,41,38,  6,21,  9, 0,44,  9,23,  0,23,29,  12,21, 20,26,49",
            "1994,  3,20, 20,29, 1,  6,21, 14,48,33,  9,23, 12,14, 1,  12,22,  8,17,50",
            "1995,  3,21,  2,15,27,  6,21, 20,35,24,  9,23, 12,14, 1,  12,22,  8,17,50",
            "1996,  3,20,  8, 4, 7,  6,21,  2,24,46,  9,22, 18, 1, 8,  12,21, 14, 6,56",
            "1997,  3,20, 13,55,42,  6,21,  8,20,59,  9,22, 23,56,49,  12,21, 20, 8, 5",
            "1998,  3,20, 19,55,42,  6,21, 14, 3,38,  9,23,  5,38,15,  12,22,  1,57,31",
            "1999,  3,21,  1,46,53,  6,21, 19,50,11,  9,23, 11,32,34,  12,22,  7,44,52",
            "2000,  3,20,  7,36,19,  6,21,  1,48,46,  9,22, 17,28,40,  12,21, 13,38,30"})
    public void test1(int year,
                      int mn1, int d1, int h1, int m1, int s1,
                      int mn2, int d2, int h2, int m2, int s2,
                      int mn3, int d3, int h3, int m3, int s3,
                      int mn4, int d4, int h4, int m4, int s4) {
        double[] jde = jde(year);
        assertTrue(check(year, mn1, d1, h1, m1, s1, jde[0], 0));
        assertTrue(check(year, mn2, d2, h2, m2, s2, jde[1], 1));
        assertTrue(check(year, mn3, d3, h3, m3, s3, jde[2], 2));
        assertTrue(check(year, mn4, d4, h4, m4, s4, jde[3], 3));
    }

    //大致时间的最大误差，单位：秒
    private static final int[] MAX_ERROR = {51, 39, 44, 41};
    private static final int checkSupportedMinYear = 1951, checkSupportedMaxYear = 2050;

    private static boolean check(int year, int month, int date, int hour, int minute, int second,
                                 double julianActual, int no) {
        if (year < checkSupportedMinYear || year > checkSupportedMaxYear)
            throw new RuntimeException(String.format("仅支持检查%d年~%d年之间的误差。",
                    checkSupportedMinYear, checkSupportedMaxYear));
        double intervalDays = JulianDay.intervalDays(
                YearMonthDay.of(year, month, date, hour, minute, second),
                JulianDay.convert(julianActual));
        return (int) (intervalDays * LevelConvert.DayHourMinuteSecond.getProduct()) < MAX_ERROR[no];
    }

    /**
     * 一些历元的四个天文季节的时间长度,单位：天
     * 数据格式：年份，春，夏，秋，冬
     */
    @ParameterizedTest
    @CsvSource({
            "-4000, 93.54, 89.18, 89.08, 93.43",
            "-3500, 93.82, 89.53, 88.82, 93.07",
            "-3000, 94.04, 89.92, 88.62, 92.67",
            "-2500, 94.19, 90.33, 88.48, 92.24",
            "-2000, 94.28, 90.76, 88.40, 91.81",
            "-1500, 94.30, 91.20, 88.38, 91.37",
            "-1000, 94.25, 91.63, 88.42, 90.94",
            " -500, 94.14, 92.05, 88.53, 90.52",
            "    0, 93.96, 92.45, 88.70, 90.14",
            "  500, 93.73, 92.82, 88.92, 89.78",
            " 1000, 93.44, 93.15, 89.18, 89.47",
            " 1500, 93.12, 93.42, 89.50, 89.20",
            " 2000, 92.76, 93.65, 89.84, 88.99",
            " 2500, 92.37, 93.81, 90.22, 88.84",
            " 3000, 91.97, 93.92, 90.61, 88.74",
            " 3500, 91,57, 93.96, 91.01, 88.71",
            " 4000, 91.17, 93.93, 91.40, 88.73",
            " 4500, 90.79, 93.84, 91.79, 88.82",
            " 5000, 90.44, 93.70, 92.15, 88.96",
            " 5500, 90.11, 93.50, 92.49, 89.14",
            " 6000, 89.82, 93.25, 92.79, 89.38",
            " 6500, 89.58, 92.97, 93.04, 89.65"})
    public void testLength(int year, double spring, double summer, double autumn, double winter) {

    }
}