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

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

string radToDMS(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);

    // 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;
    if (degrees < 0) {
        min = -min;
        sec = -sec;
    }

    char buf[100];
    sprintf(buf, "%4d°%2d'%06.3f''", deg, min, sec);
    return string(buf);
}

// 计算地球日心黄经
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::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 += nutation(T) * cos(epsilon);
    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;
}