--------------------------------------------------------
#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 double rad2rrad(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(Vector3 xyz) {
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 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面坐标;
//本含数是通用的球面坐标中心平移函数,行星计算中将反复使用;
Vector3 h2g(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)));
}
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[dt_at_length - 2]; //表中最后一年
long double t0 = dt_at[dt_at_length - 1]; //表中最后一年的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[i + 5]) break;
}