1 files changed, 88 insertions, 62 deletions

diff --git a/src/parameters.cpp b/src/parameters.cpp
index 704f083..8a99160 100644
--- a/src/parameters.cpp
+++ b/src/parameters.cpp
@@ -1,6 +1,9 @@
 #include "calendar.h"
+#include "parameters.h"
 using namespace std;
 
+//凡无特殊说明的函数，全部返回弧度制
+
 // 计算地球日心黄经
 double parameter::get_longitude(vector<double> l,double t)
 {
@@ -9,7 +12,7 @@ double parameter::get_longitude(vector<double> l,double t)
     {
         L+=l[i]*pow(t,i);
     }
-    //返回弧度制
+    //返回弧度制，太阳地心黄经=地球日心黄经+180°
     return L+pi;
 }
 
@@ -21,11 +24,11 @@ double parameter::get_latitude(vector<double> b,double t)
     {
         B+=b[i]*pow(t,i);
     }
-    //返回弧度制
     //地心黄纬=-太阳黄纬
     return -B;
 }
 
+// 计算地日距离，单位为天文单位
 double parameter::get_R(vector<double> r,double t)
 {
     double R=r[0];
@@ -36,7 +39,7 @@ double parameter::get_R(vector<double> r,double t)
     return R;
 }
 
-// 转换FK5误差，返回弧度制
+// 转换FK5误差
 double parameter::delta_FK5(double L,double B,double T)
 {
     //L转角度制
@@ -47,7 +50,7 @@ double parameter::delta_FK5(double L,double B,double T)
     L_p*=pi/180;
     while(L_p<0)
         L_p+=2*pi;
-    while(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);
@@ -56,61 +59,49 @@ double parameter::delta_FK5(double L,double B,double T)
     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::nution(double T)
+double parameter::nutation(double T)
 {
     vector<double> ang=get_angle(T);
-    fstream fr;
-    fr.open("data/nutation.txt",ios::in);
-    string str;
     double sita;
     double s,s1,s2;
     double c,c1,c2;
     double m,n,o,p,q;
 
     s=c=0;
-    while(getline(fr,str))
+    for(int i=0;i<nutation_size;i++)
     {
-        sscanf(str.c_str(),"%lf%lf%lf%lf%lf%lf%lf%lf%lf",&m,&n,&o,&p,&q,&s1,&s2,&c1,&c2);
-        sita=m*ang[0]+n*ang[1]+o*ang[2]+p*ang[3]+q*ang[4];
-
-        //sita为角度制，转换为弧度制
+        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;
-        //计算黄经章动，单位为0.0001角秒
-        s+=(s1+s2*T)*sin(sita);
-        //计算交角章动，单位为0.0001角秒
-        //c+=(c1+c2*T)*cos(sita);
+        s+=(nutation_parameters[i][5]+nutation_parameters[i][6]*T)*sin(sita);
+        //c+=(nutation_parameters[i][7]+nutation_parameters[i][8]*T)*cos(sita);
     }
-    //关闭文件
-    fr.close();
     // 计算得到的s和c单位0.0001秒，转换为弧度制
     s*=0.0001;
     s*=pi/(180*3600);
     return s;
 }
 
-//光行差修正，返回弧度制
+//光行差修正
 double parameter::aberration(double R)
 {
-    //return -20.4898/R/3600.0*pi/180.0;
     //公式，单位0.0001角秒
     double delta=-20.4898/R;
     //转换为弧度制
@@ -121,36 +112,51 @@ double parameter::aberration(double R)
 // 获取地日运行参数，L为地球日心黄经，B为地球日心黄纬，R为地日距离
 vector<vector<double>> parameter::get_parameters(double t)
 {
-    fstream fr;
-    fr.open("data/earth.txt",ios::in);
-    string s;
     double l=0;
     double a,b,c,num_tmp=0;
     vector<vector<double>> parameters;
     vector<double> tmp;
-    while(getline(fr,s))
+
+    //计算L
+    tmp.resize(0);
+    for(int i=0;i<L_size;i++)
     {
-        //当分割符为/时，表明一组参数结束，将参数存入parameters中，并清空tmp
-        if(s[0]=='/')
+        num_tmp=0;
+        for(int j=0;j<L_length[i];j++)
         {
-            parameters.push_back(tmp);
-            tmp.resize(0);
-            num_tmp=0;
+            num_tmp+=L[i][j][0]*cos(L[i][j][1]+L[i][j][2]*t);
         }
-        //当分割符为=时，表明一组参数中的一类结束，将参数存入tmp中，并清空num_tmp
-        else if(s[0]=='=')
+        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++)
         {
-            tmp.push_back(num_tmp);
-            num_tmp=0;
+            num_tmp+=B[i][j][0]*cos(B[i][j][1]+B[i][j][2]*t);
         }
-        else
+        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++)
         {
-            sscanf(s.c_str(),"%lf%lf%lf",&a,&b,&c);
-            num_tmp+=a*cos(b+c*t);
+            num_tmp+=R[i][j][0]*cos(R[i][j][1]+R[i][j][2]*t);
         }
+        tmp.push_back(num_tmp);
     }
-    //关闭文件
-    fr.close();
+    parameters.push_back(tmp);
+
+
     /* 最终传参情况：
      * parameters[0]为L系列，用来计算地球日心黄经
      * parameters[1]为B系列，用来计算地球日心黄纬
@@ -159,7 +165,7 @@ vector<vector<double>> parameter::get_parameters(double t)
     return parameters;
 }
 
-
+// 计算太阳地心黄经，外部调用，返回角度制
 double parameter::sun_longitude(double t)
 {
     //从文件中计算地日运行参数
@@ -178,7 +184,7 @@ double parameter::sun_longitude(double t)
     //转FK5误差，返回弧度制
     L+=delta_FK5(L,B,10*t);
     //章动修正，返回弧度制
-    L+=nution(10*t);
+    L+=nutation(10*t);
     //光行差修正，返回弧度制
     L+=aberration(R);
 
@@ -189,22 +195,42 @@ double parameter::sun_longitude(double t)
         L+=360;
     while(L>=360)
         L-=360;
-    if(360-L<=5)
-    {
-        L-=360;
-    }
     return L;
 }
 
-//double moon_longitude(double t)
-//{
-//    //文件中获取地月运行参数
-//    vector<vector<double>> p=get_moon_parameters(t);
-//    vector<double> l=p[0],b=p[1],r=p[2];
+// 计算月球地心黄经，外部调用，返回角度制
+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 L=get_moon_longitude(l,t);
-//    double B=get_moon_latitude(b,t);
-//    double R=get_moon_R(r,t);
+    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);//地球章动修正
+
+    while(L<0)
+        L+=360;
+    while(L>=360)
+        L-=360;
+    return L;
+}