1 files changed, 119 insertions, 128 deletions

diff --git a/src/parameters.cpp b/src/parameters.cpp
index 8a99160..57d98c2 100644
--- a/src/parameters.cpp
+++ b/src/parameters.cpp
@@ -1,130 +1,126 @@
-#include "calendar.h"
 #include "parameters.h"
+#include "calendar.h"
 using namespace std;
 
 //凡无特殊说明的函数，全部返回弧度制
 
 // 计算地球日心黄经
-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);
+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+pi;
+    return L + 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);
+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);
+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/pi;
+double parameter::delta_FK5(double L, double B, double T) {
+    // L转角度制
+    L *= 180 / pi;
     //计算L'，角度制
-    double L_p=L-1.397*T-0.00031*T*T;
+    double L_p = L - 1.397 * T - 0.00031 * T * T;
     //计算L'，弧度制
-    L_p*=pi/180;
-    while(L_p<0)
-        L_p+=2*pi;
-    while(L_p>=2*pi)
-        L_p-=2*pi;
+    L_p *= pi / 180;
+    while (L_p < 0)
+        L_p += 2 * pi;
+    while (L_p >= 2 * pi)
+        L_p -= 2 * pi;
     //计算delta_L，单位为角秒
-    double delta_L=-0.09033+0.03916*(cos(L_p)+sin(L_p))*tan(B);
+    double delta_L = -0.09033 + 0.03916 * (cos(L_p) + sin(L_p)) * tan(B);
     //转换为弧度制
-    delta_L*=pi/(180*3600);
+    delta_L *= pi / (180 * 3600);
     return delta_L;
 }
 
 //章动有关角，结果均为角度制
-vector<double> parameter::get_angle(double T)
-{
+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(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(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(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(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);
+    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 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*=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);
+    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 *= 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*=pi/(180*3600);
+    s *= 0.0001;
+    s *= pi / (180 * 3600);
     return s;
 }
 
 //光行差修正
-double parameter::aberration(double R)
-{
+double parameter::aberration(double R) {
     //公式，单位0.0001角秒
-    double delta=-20.4898/R;
+    double delta = -20.4898 / R;
     //转换为弧度制
-    delta*=pi/(180*3600);
+    delta *= 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>> 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);
+    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);
     }
@@ -132,12 +128,10 @@ vector<vector<double>> parameter::get_parameters(double t)
 
     //计算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);
+    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);
     }
@@ -145,18 +139,15 @@ vector<vector<double>> parameter::get_parameters(double t)
 
     //计算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);
+    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系列，用来计算地球日心黄纬
@@ -166,71 +157,71 @@ vector<vector<double>> parameter::get_parameters(double t)
 }
 
 // 计算太阳地心黄经，外部调用，返回角度制
-double parameter::sun_longitude(double t)
-{
+double parameter::sun_longitude(double t) {
     //从文件中计算地日运行参数
-    vector<vector<double>> p=get_parameters(t);
+    vector<vector<double>> p = get_parameters(t);
 
     //取出各类参数
-    vector<double> l=p[0],b=p[1],r=p[2];
+    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);
-
+    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 += delta_FK5(L, B, 10 * t);
     //章动修正，返回弧度制
-    L+=nutation(10*t);
+    L += nutation(10 * t);
     //光行差修正，返回弧度制
-    L+=aberration(R);
-
-
-    //L转角度制
-    L*=180/pi;
-    while(L<0)
-        L+=360;
-    while(L>=360)
-        L-=360;
+    L += aberration(R);
+
+    // L转角度制
+    L *= 180 / pi;
+    while (L < 0)
+        L += 360;
+    while (L >= 360)
+        L -= 360;
     return L;
 }
 
 // 计算月球地心黄经，外部调用，返回角度制
-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*=pi/180;
-        EI+=moon_longitude_parameters[i][4]*sin(sita)*pow(E,fabs(moon_longitude_parameters[i][1]));
+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 *= 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;
+    double a1, a2;
+    a1 = 119.75 + 131.849 * T;
+    a2 = 53.09 + 479264.290 * T;
 
-    EI+=3958.0*sin(a1*pi/180);
-    EI+=1962.0*sin((L-angles[3])*pi/180);
-    EI+=318.0*sin(a2*pi/180);
+    EI += 3958.0 * sin(a1 * pi / 180);
+    EI += 1962.0 * sin((L - angles[3]) * pi / 180);
+    EI += 318.0 * sin(a2 * pi / 180);
 
-    L+=EI/1000000;
-    L+=180/pi*nutation(T);//地球章动修正
+    L += EI / 1000000;
+    L += 180 / pi * nutation(T); //地球章动修正
 
-    while(L<0)
-        L+=360;
-    while(L>=360)
-        L-=360;
+    while (L < 0)
+        L += 360;
+    while (L >= 360)
+        L -= 360;
     return L;
 }