地球相关参数计算。
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));
}
}
}