path: root/src/parameters.cpp

// 凡无特殊说明的函数，全部返回弧度制

#include "parameters.h"
#include "calendar.h"
#include <utility>
using namespace std;

void radToDMS(const char *message, int lowerBound, int upperBound,
              double radians) {
    // Convert radians to degrees
    double degrees = radians * 180.0 / M_PI;

    while (degrees < lowerBound)
        degrees += (upperBound - lowerBound);
    while (degrees >= upperBound)
        degrees -= (upperBound - lowerBound);

    bool flag = degrees < 0;
    if (flag)
        degrees = -degrees;
    // Calculate degrees, minutes, and seconds
    int deg = static_cast<int>(degrees);
    double min_sec = (degrees - deg) * 60.0;
    int min = static_cast<int>(min_sec);
    double sec = (min_sec - min) * 60.0;

    // Output the result
    printf("%s: %c%d°%d'%06.3f''\n", message, flag ? '-' : ' ', deg, min, sec);
}

// 计算地球日心黄经
double parameter::get_longitude(vector<double> l, double t) {
    double L = l[0];
    for (int i = 1; i < l.size(); i++) {
        L += l[i] * pow(t, i);
    }
    //返回弧度制，太阳地心黄经=地球日心黄经+180°
    return L + M_PI;
}

// 计算地球日心黄纬
double parameter::get_latitude(vector<double> b, double t) {
    double B = b[0];
    for (int i = 1; i < b.size(); i++) {
        B += b[i] * pow(t, i);
    }
    //地心黄纬=-太阳黄纬
    return -B;
}

// 计算地日距离，单位为天文单位
double parameter::get_R(vector<double> r, double t) {
    double R = r[0];
    for (int i = 1; i < r.size(); i++) {
        R += r[i] * pow(t, i);
    }
    return R;
}

// 转换FK5误差
double parameter::delta_FK5(double L, double B, double T) {
    // L转角度制
    L *= 180 / M_PI;
    //计算L'，角度制
    double L_p = L - 1.397 * T - 0.00031 * T * T;
    //计算L'，弧度制
    L_p *= M_PI / 180;
    while (L_p < 0)
        L_p += 2 * M_PI;
    while (L_p >= 2 * M_PI)
        L_p -= 2 * M_PI;
    //计算delta_L，单位为角秒
    double delta_L = -0.09033 + 0.03916 * (cos(L_p) + sin(L_p)) * tan(B);
    //转换为弧度制
    delta_L *= M_PI / (180 * 3600);
    return delta_L;
}

//章动有关角，结果均为角度制
vector<double> parameter::get_angle(double T) {
    vector<double> ang;
    //平距角
    ang.push_back(297.85036 + 445267.111480 * T - 0.0019142 * T * T +
                  T * T * T / 189474);
    //日地平近点角
    ang.push_back(357.52772 + 35999.050340 * T - 0.0001603 * T * T -
                  T * T * T / 300000);
    //月球平近点角
    ang.push_back(134.96298 + 477198.867398 * T + 0.0086972 * T * T +
                  T * T * T / 56250);
    //月球纬度参数
    ang.push_back(93.27191 + 483202.017538 * T - 0.0036825 * T * T +
                  T * T * T / 327270);
    //黄道与月球平轨道交点黄经
    ang.push_back(125.04452 - 1934.136261 * T + 0.0020708 * T * T +
                  T * T * T / 450000);
    return ang;
}

//章动修正
double parameter::nutation(double T) {
    vector<double> ang = get_angle(T);
    double sita;
    double s, s1, s2;
    double c, c1, c2;
    double m, n, o, p, q;

    s = c = 0;
    for (int i = 0; i < nutation_size; i++) {
        sita = nutation_parameters[i][0] * ang[0] +
               nutation_parameters[i][1] * ang[1] +
               nutation_parameters[i][2] * ang[2] +
               nutation_parameters[i][3] * ang[3] +
               nutation_parameters[i][4] * ang[4];
        sita *= M_PI / 180;
        s += (nutation_parameters[i][5] + nutation_parameters[i][6] * T) *
             sin(sita);
        // c+=(nutation_parameters[i][7]+nutation_parameters[i][8]*T)*cos(sita);
    }
    // 计算得到的s和c单位0.0001秒，转换为弧度制
    s *= 0.0001;
    s *= M_PI / (180 * 3600);
    return s;
}

//光行差修正
double parameter::aberration(double R) {
    //公式，单位0.0001角秒
    double delta = -20.4898 / R;
    //转换为弧度制
    delta *= M_PI / (180 * 3600);
    return delta;
}

// 获取地日运行参数，L为地球日心黄经，B为地球日心黄纬，R为地日距离
vector<vector<double>> parameter::get_parameters(double t) {
    double l = 0;
    double a, b, c, num_tmp = 0;
    vector<vector<double>> parameters;
    vector<double> tmp;

    //计算L
    tmp.resize(0);
    for (int i = 0; i < L_size; i++) {
        num_tmp = 0;
        for (int j = 0; j < L_length[i]; j++) {
            num_tmp += L[i][j][0] * cos(L[i][j][1] + L[i][j][2] * t);
        }
        tmp.push_back(num_tmp);
    }
    parameters.push_back(tmp);

    //计算B
    tmp.resize(0);
    for (int i = 0; i < B_size; i++) {
        num_tmp = 0;
        for (int j = 0; j < B_length[i]; j++) {
            num_tmp += B[i][j][0] * cos(B[i][j][1] + B[i][j][2] * t);
        }
        tmp.push_back(num_tmp);
    }
    parameters.push_back(tmp);

    //计算R
    tmp.resize(0);
    for (int i = 0; i < R_size; i++) {
        num_tmp = 0;
        for (int j = 0; j < R_length[i]; j++) {
            num_tmp += R[i][j][0] * cos(R[i][j][1] + R[i][j][2] * t);
        }
        tmp.push_back(num_tmp);
    }
    parameters.push_back(tmp);

    /* 最终传参情况：
     * parameters[0]为L系列，用来计算地球日心黄经
     * parameters[1]为B系列，用来计算地球日心黄纬
     * parameters[2]为R系列，用来计算地日距离
     */
    return parameters;
}

// 计算太阳地心黄经，外部调用，返回角度制
pair<double, double> parameter::sun_longitude(double t) {
    //从文件中计算地日运行参数
    vector<vector<double>> p = get_parameters(t);

    //取出各类参数
    vector<double> l = p[0], b = p[1], r = p[2];

    //计算地球日心黄经纬，返回弧度制/天文单位
    double L = get_longitude(l, t);
    double B = get_latitude(b, t);
    double R = get_R(r, t);

    //以下修正需要的都是儒略世纪数而非千年数
    //转FK5误差，返回弧度制
    L += delta_FK5(L, B, 10 * t);
    //章动修正，返回弧度制
    L += nutation(10 * t);
    //光行差修正，返回弧度制
    L += aberration(R);

    // 返回太阳地心黄经和黄纬
    return make_pair(L, -B);
}

// 计算月球地心黄经，外部调用，返回角度制
double parameter::moon_longitude(double t) {
    double T = 10 * t;
    vector<double> angles = get_angle(T);
    double a, b, c, d, moon_longitude_a;

    double L = 218.3164591 + 481267.88134236 * T - 0.0013268 * T * T +
               T * T * T / 538841 - T * T * T * T / 65194000; //月球平黄经
    double EI = 0, sita = 0;
    double E = 1 - 0.002516 * T - 0.0000074 * T * T; //地球离心率

    for (int i = 0; i < moon_size; i++) {
        sita = moon_longitude_parameters[i][0] * angles[0] +
               moon_longitude_parameters[i][1] * angles[1] +
               moon_longitude_parameters[i][2] * angles[2] +
               moon_longitude_parameters[i][3] * angles[3];
        sita *= M_PI / 180;
        EI += moon_longitude_parameters[i][4] * sin(sita) *
              pow(E, fabs(moon_longitude_parameters[i][1]));
    }

    //地外行星影响月球地心黄经修正项
    double a1, a2;
    a1 = 119.75 + 131.849 * T;
    a2 = 53.09 + 479264.290 * T;

    EI += 3958.0 * sin(a1 * M_PI / 180);
    EI += 1962.0 * sin((L - angles[3]) * M_PI / 180);
    EI += 318.0 * sin(a2 * M_PI / 180);

    L += EI / 1000000;
    L += 180 / M_PI * nutation(T); //地球章动修正

    while (L < 0)
        L += 360;
    while (L >= 360)
        L -= 360;
    return L;
}

double parameter::getHourAngle(double julianKiloYear, double longitude,
                               double latitude, double alpha) {
    double T = julianKiloYear * 10;
    double GST = 280.46061837 + 360.98564736629 * 36525 * T +
                 0.000387933 * T * T - T * T * T / 38710000;
    GST *= M_PI / 180;
    return GST + longitude - alpha; // 返回弧度制
}

double parameter::get_epsilon(double T) {
    double epsilon = 23.439291111 - 0.01300416667 * T -
                     (1.638888889e-7) * T * T + (5.036111111e-7) * T * T * T;
    return epsilon * M_PI / 180;
}