本文给出天文计算中常用的坐标系统之间的转换。
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));
}
}
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
