这是C语言的不？能封装成万年历等？
--------------------------------------------------------
#include "eph.h"
#include "const.h"
#include <math.h>
#include <string>
//=================================数学工具=========================================
//对超过0-2PI的角度转为0-2PI;
long double rad2mrad(long double v)
{
v = fmod(v, (2 * PI));
if (v < 0) return v + 2 * PI;
return v;
}
//对超过-PI到PI的角度转为-PI到PI;
long doublerad2rrad(long double v)
{
v = fmod(v, (2 * PI));
if (v <= -PI) return v + 2 * PI;
if (v > PI) return v - 2 * PI;
return v;
}
//临界余数(a与最近的整倍数b相差的距离);
long double mod2(long double a, long double b)
{
long double c = a / b;
c -= floor(c);
if (c > 0.5) c -= 1;
return c*b;
}
//球面转直角坐标;
Vector3 llr2xyz(Vector3 JW)
{
Vector3 r;
long double J = JW.x, W = JW.y, R = JW.z;
r.x = R*cos(W)*cos(J);
r.y = R*cos(W)*sin(J);
r.z = R*sin(W);
return r;
}
//直角坐标转球;
Vector3 xyz2llr(Vector3xyz) {
Vector3 r;
long double x = xyz.x, y = xyz.y, z = xyz.z;
r.z = sqrt(x*x + y*y + z*z);
r.y = asin(z / r.z);
r.x = rad2mrad(atan2(y, x));
return r;
}
//球面坐标旋转;
Vector3 llrConv(Vector3 JW, long double E) {
Vector3 r;
long double J = JW.x, W = JW.y;
r.x = atan2(sin(J)*cos(E) - tan(W)*sin(E), cos(J));
r.y = asin(cos(E)*sin(W) + sin(E)*cos(W)*sin(J));
r.z = JW.z;
r.x = rad2mrad(r.x);
return r;
}
//赤道坐标转为地平坐标;
Vector3 CD2DP(Vector3 z, long double L, long double fa, long double gst)
{
//转到相对于地平赤道分点的赤道坐标;
Vector3 a(z.x + PI / 2 - gst - L, z.y, z.z);
a = llrConv(a, PI / 2 - fa);
a.x = rad2mrad(PI / 2 - a.x);
return a;
}
//求角度差;
long double j1_j2(long double J1, long double W1, long double J2, long double W2) {
long double dJ = rad2rrad(J1 - J2), dW = W1 - W2;
if (fabs(dJ) < 1 / 1000 && fabs(dW) < 1 / 1000) {
dJ *= cos((W1 + W2) / 2);
return sqrt(dJ*dJ + dW*dW);
}
return acos(sin(W1)* sin(W2) + cos(W1)*cos(W2)*cos(dJ));
}
//日心球面转地心球面,Z星体球面坐标,A地qiuqiu面坐标;
//本含数是通用的球面坐标中心平移函数,行星计算中将反复使用;
Vector3h2g(Vector3 z, Vector3 a) {
a = llr2xyz(a); //地球
z = llr2xyz(z); //星体
z.x -= a.x; z.y -= a.y; z.z -= a.z;
return xyz2llr(z);
}
//视差角(不是视差);
long double shiChaJ(long double gst, long double L, long double fa, long double J, long double W) {
long double H = gst + L - J; //天体的时角;
return rad2mrad(atan2(sin(H), tan(fa)*cos(W) - sin(W)*cos(H)));
}
//=================================deltat T计算=====================================
static constexpr long double dt_at[] = { // TD - UT1 计算表
-4000, 108371.7, -13036.80, 392.000, 0.0000,
-500, 17201.0, -627.82, 16.170, -0.3413,
-150, 12200.6, -346.41, 5.403, -0.1593,
150, 9113.8, -328.13, -1.647, 0.0377,
500, 5707.5, -391.41, 0.915, 0.3145,
900, 2203.4, -283.45, 13.034, -0.1778,
1300, 490.1, -57.35, 2.085, -0.0072,
1600, 120.0, -9.81, -1.532, 0.1403,
1700, 10.2, -0.91, 0.510, -0.0370,
1800, 13.4, -0.72, 0.202, -0.0193,
1830, 7.8, -1.81, 0.416, -0.0247,
1860, 8.3, -0.13, -0.406, 0.0292,
1880, -5.4, 0.32, -0.183, 0.0173,
1900, -2.3, 2.06, 0.169, -0.0135,
1920, 21.2, 1.69, -0.304, 0.0167,
1940, 24.2, 1.22, -0.064, 0.0031,
1960, 33.2, 0.51, 0.231, -0.0109,
1980, 51.0, 1.29, -0.026, 0.0032,
2000, 63.87, 0.1, 0, 0,
2005, 64.7, 0.4, 0, 0, //一次项记为x,则 10x=0.4秒/年*(2015-2005),解得x=0.4
2015, 69
};
long double dt_ext(long double y, int jsd) { long double dy = (y - 1820) / 100; return -20 + jsd*dy*dy; } //二次曲线外推
long double dt_calc(long double y) { //计算世界时与原子时之差,传入年
int dt_at_length = sizeof(dt_at) / sizeof(long double);
long double y0 = dt_at; //表中最后一年
long double t0 = dt_at; //表中最后一年的deltatT
if (y >= y0) {
int jsd = 31; //sjd是y1年之后的加速度估计。瑞士星历表jsd=31,NASA网站jsd=32,skmap的jsd=29
if (y > y0 + 100) return dt_ext(y, jsd);
long double v = dt_ext(y, jsd); //二次曲线外推
long double dv = dt_ext(y0, jsd) - t0;//ye年的二次外推与te的差
return v - dv*(y0 + 100 - y) / 100;
}
int i;
const long double* d = dt_at;
for (i = 0; i < dt_at_length; i += 5)
{
if (y < d) break;
}
long double t1 = (y - d) / (d - d) * 10, t2 = t1*t1, t3 = t2*t1;
return d + d * t1 + d * t2 + d * t3;
}
//传入儒略日(J2000起算),计算TD-UT(单位:日)
long double dt_T(long double t)
{
return dt_calc(t / 365.2425 + 2000) / 86400.0;
}
//=================================岁差计算=========================================
//IAU1976岁差表
static constexpr long double preceTab_IAU1976[] = {
0, 5038.7784, -1.07259, -0.001147, //fi
84381.448, 0, +0.05127, -0.007726, //w
0, +4.1976, +0.19447, -0.000179, //P
0, -46.8150, +0.05059, +0.000344, //Q
84381.448, -46.8150, -0.00059, +0.001813, //E
0, +10.5526, -2.38064, -0.001125, //x
0, 47.0028, -0.03301, +0.000057, //pi(自导出),根据P=sin(pi)*sin(II)
629554.886, -869.8192, +0.03666, -0.001504, //II(自导出),根据Q=sin(pi)*cos(II)
0, 5029.0966, +1.11113, +0.000006, //p
0, 2004.3109, -0.42665, -0.041833, //th
0, 2306.2181, +0.30188, +0.017998, //Z
0, 2306.2181, +1.09468, +0.018203//z
};
//IAU2000岁差表
static constexpr long double preceTab_IAU2000[] = {
0, 5038.478750, -1.07259, -0.001147, 0, 0, //fi
84381.448, -0.025240, +0.05127, -0.007726, 0, 0, //w,
0, +4.1976, +0.19447, -0.000179, 0, 0, //P,
0, -46.8150, +0.05059, +0.000344, 0, 0, //Q,
84381.448, -46.84024, -0.00059, +0.001813, 0, 0, //E,
0, +10.5526, -2.38064, -0.001125, 0, 0, //x
0, 47.0028, -0.03301, +0.000057, 0, 0, //pi(自导出)
629554.886, -869.8192, +0.03666, -0.001504, 0, 0, //II(自导出)
0, 5028.79695, +1.11113, +0.000006, 0, 0, //p
0, 2004.1917476, -0.4269353, -0.0418251, -0.0000601, -0.0000001, //th
+2.5976176, 2306.0809506, +0.3019015, +0.0179663, -0.0000327, -0.0000002, //Z
-2.5976176, 2306.0803226, +1.0947790, +0.0182273, +0.0000470, -0.0000003//z
};
这是c++的，可以直接在火山代码里内嵌吧
页:
[1]