天文算法--坐标变换

Ancony_

于 2024-04-02 23:48:05 发布

阅读量605

点赞数 2

分类专栏：中国农历文章标签：算法农历 java

本文链接：https://blog.csdn.net/Ancony_/article/details/137299961

版权

中国农历专栏收录该内容

8 篇文章 0 订阅

订阅专栏

本文给出天文计算中常用的坐标系统之间的转换。

package cn.ancony.chinese_calendar;

import static java.lang.Math.*;

/**
 * 坐标变换。 五个坐标系统：赤道坐标、黄道坐标、银道坐标、本地地平坐标、观测点坐标。
 * 它们使用如下参数表示：
 * 赤道坐标：
 * α = 赤径。一般表达纬时间单位，用时分秒表示。计算时要先转为度，再转为弧度。
 * . α1950 = 涉及 B1950.0 标准分点的赤经。
 * . α2000 = 涉及 J2000.0 标准分点的赤经
 * δ = 赤纬。天赤道以北纬正，以南为负。
 * . δ1950 = 涉及 B1950.0 标准分点的赤纬
 * . δ2000 = 涉及 J2000.0 标准分点的赤纬。
 * 黄道坐标：
 * λ = 黄经。从春风点，沿黄道测量的经度。
 * β = 黄纬。黄道以北为正，以南为负。
 * 银道坐标：
 * l = 银经。
 * b = 银纬。
 * 本地地平坐标：
 * A = 地平经度(方位角)。由南向西测量
 * h = 地平纬度，地平线以上为正，以下为负
 * 观测点坐标参数：
 * H = 本地时角，从南向西测量。
 * φ = 观测者(站)纬度,北半球为正,南半球为负。
 * 如果θ是本地恒星时，θo是格林尼治恒星时，L 是观者站经度（从格林尼治向西为正，东为负），
 * 那么本地时角计算如下：(如果α含章动效果，那么 H 也含章动)
 * H = θ - α 或 H =θo - L – α
 * 其他参数：
 * ε = 黄赤交角。黄道与天赤道的夹角。如果使用视赤经及视赤纬（受光行差及章动影响），那么计算时就要用到真黄赤交角ε+Δε
 * . ε2000 = 23°26′21".448 = 23°.439291 J2000.0 标准分点坐标
 * . ε1950 = 23°.4457889 历元 B1950.0 标准分点坐标
 */
public class CoordinateTransformation {
    public static double
            epsilon2000 = toRadians(23.439291),//ε2000
            epsilon1950 = toRadians(23.4457889);//ε1950

    /**
     * 从赤道坐标转到黄道坐标的公式
     * <p>
     * tan(λ)=(sin(α)*cos(ε)+tan(δ)*sin(ε))/cos(α)  (12.1)
     * <p>
     * sin(β)=sin(δ)*cos(ε)-cos(δ)*sin(ε)*sin(α)    (12.2)
     *
     * @param alpha   α 赤经,单位:弧度
     * @param delta   δ 赤纬,单位:弧度
     * @param epsilon ε 黄赤交角,单位:弧度
     * @return 返回 [λ,β] 单位:弧度
     */
    public static double[] func1(double alpha, double delta, double epsilon) {
        double sinA = sin(alpha), sinE = sin(epsilon), cosE = cos(epsilon);
        return new double[]{
                atan2(sinA * cosE + tan(delta) * sinE, cos(alpha)),
                asin(sin(delta) * cosE - cos(delta) * sinE * sinA)};
    }

    /**
     * 黄道转赤道
     * <p>
     * tan(α)=(sin(λ)*cos(ε)-tan(β)*sin(ε))/cos(λ)      (12.3)
     * <p>
     * sin(δ)=sin(β)*cos(ε)+cos(β)*sin(ε)*sin(λ)        (12.4)
     *
     * @param lambda  λ,单位:弧度
     * @param beta    β,单位:弧度
     * @param epsilon ε,单位:弧度
     * @return [α, δ],单位:弧度
     */
    public static double[] func2(double lambda, double beta, double epsilon) {
        double sinL = sin(lambda), sinE = sin(epsilon), cosE = cos(epsilon);
        return new double[]{
                atan2(sinL * cosE - tan(beta) * sinE, cos(lambda)),
                asin(sin(beta) * cosE + cos(beta) * sinE * sinL)};
    }

    /**
     * 由观测点计算本地地平坐标：
     * <p>
     * tan(A)=sin(H)/(cos(H)*sin(φ)-tan(δ)*cos(φ))  (12.5)
     * <p>
     * sin(h)=sin(φ)*sin(δ)+cos(φ)*cos(δ)*cos(H)    (12.6)
     *
     * @param h     H 本地时角，从南向西测量。,单位:弧度
     * @param phi   φ,单位:弧度
     * @param delta δ  赤纬,单位:弧度
     * @return [A, h],单位:弧度。A-从南点起算。如果希望从北点计算，那么将A加上π即可。
     * h没有考虑大气折射的影响，也没有考虑行星视差。
     */
    public static double[] func3(double h, double phi, double delta) {
        double cosH = cos(h), sinP = sin(phi), cosP = cos(phi);
        return new double[]{
                atan2(sin(h), cosH * sinP - tan(delta) * cosP),
                asin(sinP * sin(delta) + cosP * cos(delta) * cosH)};
    }

    /**
     * 地平坐标转到赤道坐标
     * <p>
     * tan(H) = sin(A)/( cos(A)*sin(φ) +tan(h)*cos(φ))
     * <p>
     * sin(δ)= sin(φ)*sin(h) - cos(φ)*cos(h)*cos(A)
     *
     * @param a   A,单位:弧度
     * @param phi φ,单位:弧度
     * @param h   H 本地时角，从南向西测量。,单位:弧度
     * @return [H, δ],单位:弧度
     */
    public static double[] func4(double a, double phi, double h) {
        double cosA = cos(a), sinP = sin(phi), cosP = cos(phi);
        return new double[]{
                atan2(sin(a), cosA * sinP + tan(h) * cosP),
                asin(sinP * sin(h) - cosP * cos(h) * cosA)};
    }

    /**
     * 银河(银河系)北极的坐标。(B1950.0 标准赤道系统,IAU于1959年定义)
     * 银经的原点在银道上，该点距银道与 B1950.0 赤道升点 33 度。
     */
    public static double
            alpha1950 = toRadians(192.25),//α1950 = 12h 49m = 192°.25
            delta1950 = toRadians(27.4);//δ1950 = +27°.4

    /**
     * 从 B1950.0 标准分点赤道坐标转到银道坐标
     * <p>
     * 如果给定恒星 2000.0 的平位置而不是 1950.0 平位置，那么应将α2000 转为α1950，将δ2000 转为δ1950
     * <p>
     * tan(x)=sin(192.25-α)/(cos(192.25-α)*sin(27.4)-tan(δ)*cos(27.4))  (12.7)
     * <p>
     * l = 303° - x
     * <p>
     * sin(b)=sin(δ)*sin(27.4)+cos(δ)*cos(27.4)*cos(192.25- α)          (12.8)
     *
     * @param a1950 α B1950.0 标准分点,单位:弧度
     * @param d1950 δ B1950.0 标准分点,单位:弧度
     * @return [l, b],单位:弧度
     */
    public static double[] func5(double a1950, double d1950) {
        double diffA = alpha1950 - a1950, sinDelta1950 = sin(delta1950),
                cosDelta1950 = cos(delta1950), cosDA = cos(diffA);
        return new double[]{
                toRadians(303) - atan2(sin(diffA), cosDA * sinDelta1950 - tan(d1950) * cosDelta1950),
                asin(sin(d1950) * sinDelta1950 + cos(d1950) * cosDelta1950 * cosDA)};
    }

    /**
     * 从银道坐标转到 B1950.0 标准分点坐标：
     *
     * @param l 银经,单位:弧度
     * @param b 银纬,单位:弧度
     * @return [α, δ] 均为B1950.0 标准分点坐标,单位:弧度
     */
    public static double[] func6(double l, double b) {
        double diff = l - toRadians(123), cosD = cos(diff),
                sinDelta1950 = sin(delta1950), cosDelta1950 = cos(delta1950);
        return new double[]{
                toRadians(125) + atan2(sin(diff), cosD * sinDelta1950 - tan(b) * cosDelta1950),
                asin(sin(b) * sinDelta1950 + cos(b) * cosDelta1950 * cosD)};
    }

    /**
     * 黄道和地平线
     * <p>
     * λ:在地平线(地面与天球交的大圆)与黄道相交的两点的黄经。λ是其中一点，另一点和它相距π。
     * <p>
     * tan(λ)=(-cos(θ))/(sin(ε)tan(φ)+cos(ε)*sin(θ))        (12.9)
     * <p>
     * I: 黄道面与地平面的夹角。未考虑视差。不是太阳周日运动的平面与地平面的夹角。而是周年运动平面(黄道面)与地平的夹角。
     * <p>
     * cos(I)=cos(ε)*sin(φ)-sin(ε)*cos(φ)*sin(θ)            (12.10)
     * <p>
     * 在一个恒星日期间，角度 I 在两个极限值之间。例如，纬度 48°00'N，ε=23°26'，那么 I 的两个极限值是：
     * <p>
     * θ = 90°时， I = 90°- φ + ε = 65°26′
     * <p>
     * θ = 270°时，I = 90°- φ - ε = 18°34′
     *
     * @param theta   θ是本地恒星时,单位:弧度
     * @param phi     φ是观测站的纬度,单位:弧度
     * @param epsilon ε是黄赤交角,单位:弧度
     * @return [λ, I] 单位:弧度
     */
    public static double[] func7(double theta, double phi, double epsilon) {
        double sinE = sin(epsilon), cosE = cos(epsilon), sinT = sin(theta);
        return new double[]{
                atan2(-cos(theta), sinE * tan(phi) + cosE * sinT),
                acos(cosE * sin(phi) - sinE * cos(phi) * sinT)};
    }
}

坐标类，可以使用时间表示，也可以使用度分秒表示。

package cn.ancony.chinese_calendar;

import lombok.Getter;

import java.text.DecimalFormat;

import static java.lang.Math.toRadians;

/**
 * celestial sphere longitude:天球经度
 * Coordinate：坐标，两种表示方法，一种是使用时分秒，一种是使用度。
 * todo 1,支持精度的表示. 2,使用转换系统. 3,为负数的情况
 */
@Getter
public class Coordinate {
    private final int hour, minute, degree, arcMinute;
    private final double second, arcSecond, degrees;

    private static final LevelConvert lc = LevelConvert.DegreeMinuteSecond;
    private static final int lastScale = 3;

    //从复合量构建
    private Coordinate(int a, int b, double c, boolean angle) {
        double high = lc.toHighest(a, b, c);
        if (angle) {
            degree = a;
            arcMinute = b;
            arcSecond = c;
            degrees = high;
            Number[] numbers = lc.fromHighest(degrees / 15, lastScale);
            hour = (int) numbers[0];
            minute = (int) numbers[1];
            second = numbers[2].doubleValue();
        } else {
            hour = a;
            minute = b;
            second = c;
            degrees = high * 15;
            Number[] numbers = lc.fromHighest(degrees, lastScale);
            degree = (int) numbers[0];
            arcMinute = (int) numbers[1];
            arcSecond = numbers[2].doubleValue();
        }
    }

    //从单一量构建
    private Coordinate(double v, boolean angle) {
        if (angle) {
            degrees = v;
            Number[] numbers = lc.fromHighest(v, lastScale);
            degree = (int) numbers[0];
            arcMinute = (int) numbers[1];
            arcSecond = numbers[2].doubleValue();
            Number[] numbers2 = lc.fromHighest(v / 15, lastScale);
            hour = (int) numbers2[0];
            minute = (int) numbers2[1];
            second = numbers2[2].doubleValue();
        } else {
            degrees = v * 15;
            Number[] numbers = lc.fromHighest(degrees, lastScale);
            degree = (int) numbers[0];
            arcMinute = (int) numbers[1];
            arcSecond = numbers[2].doubleValue();
            Number[] numbers2 = lc.fromHighest(v, lastScale);
            hour = (int) numbers2[0];
            minute = (int) numbers2[1];
            second = numbers2[2].doubleValue();
        }
    }

    //使用时分秒表示
    public String toString() {
        int s = (int) second;
        DecimalFormat df = new DecimalFormat("#.###"); // 设置小数点后三位
        return String.format("%dh%dm%ds.%s", hour, minute, s, df.format(second - s).substring(2));
    }

    //使用度表示
    public String toDegreeString() {
        int d = (int) degrees;
        DecimalFormat df = new DecimalFormat("#.###"); // 设置小数点后三位
        return String.format("%d°.%s", d, df.format(degrees - d).substring(2));
    }

    //得到弧度
    public double radian() {
        return toRadians(degrees);
    }

    public static Coordinate ofTime(int hour, int minute, double second) {
        return new Coordinate(hour, minute, second, false);
    }

    //单一量以最高单位计。
    public static Coordinate fromTime(double v) {
        return new Coordinate(v, false);
    }

    public static Coordinate ofAngle(int degree, int arcMinute, double arcSecond) {
        return new Coordinate(degree, arcMinute, arcSecond, true);
    }

    public static Coordinate fromAngle(double v) {
        return new Coordinate(v, true);
    }
}

测试方法：

package cn.ancony.javafx.chinese_calendar;

import cn.ancony.chinese_calendar.Coordinate;
import org.junit.jupiter.params.ParameterizedTest;
import org.junit.jupiter.params.provider.CsvSource;

import static cn.ancony.chinese_calendar.CoordinateTransformation.*;
import static java.lang.Math.*;
import static org.junit.jupiter.api.Assertions.assertTrue;

public class TestCoordinateTransformation {
    /**
     * 例 12.a——计算恒星 Pollux(β Gem)的黄道坐标，它的赤道坐标是：
     * α2000 = 7h 45m 18s.946， δ2000 = +28°.026183
     * 黄道坐标为：λ=113°.215630 β = +6°.684170
     * 请进行正算和反算。
     */
    @ParameterizedTest
    @CsvSource({"7, 45, 18.946,   28.026183,   113.215630, 6.684170"})
    public void testA(int h, int m, double s, double delta, double lambda, double beta) {
        double[] res = func1(Coordinate.ofTime(h, m, s).radian(), toRadians(delta), epsilon2000);
        double precision = 1e-7;//精度
        assertTrue(
                abs(toRadians(lambda) - res[0]) < precision
                        && abs(toRadians(beta) - res[1]) < precision);
        double[] res2 = func2(res[0], res[1], epsilon2000);
        Coordinate coordinate = Coordinate.fromAngle(toDegrees(res2[0]));
        assertTrue(
                h == coordinate.getHour()
                        && m == coordinate.getMinute()
                        && s == coordinate.getSecond()
                        && toDegrees(res2[1]) - delta < precision);
    }

    /**
     * 例 12.c——ε=23°.44，φ=+51°，θ = (5h 00m) =75°，求λ和I。
     * λ =169°21′和 λ = 349°21′。I = 62°
     */
    @ParameterizedTest
    @CsvSource({"5, 51, 23.44,   169, 21,  62"})
    public void testC(int hour, int phi, double epsilon, int degree, int minute, int d2) {
        double[] res = func7(Coordinate.ofTime(hour, 0, 0).radian(),
                toRadians(phi), toRadians(epsilon));
    }
}