1 files changed, 210 insertions, 0 deletions

diff --git a/src/parameters.cpp b/src/parameters.cpp
new file mode 100644
index 0000000..704f083
--- /dev/null
+++ b/src/parameters.cpp
@@ -0,0 +1,210 @@
+#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);
+    }
+    //返回弧度制
+    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);
+    }
+    //返回弧度制
+    //地心黄纬=-太阳黄纬
+    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/pi;
+    //计算L'，角度制
+    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;
+    //计算delta_L，单位为角秒
+    double delta_L=-0.09033+0.03916*(cos(L_p)+sin(L_p))*tan(B);
+    //转换为弧度制
+    delta_L*=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::nution(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))
+    {
+        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*=pi/180;
+        //计算黄经章动，单位为0.0001角秒
+        s+=(s1+s2*T)*sin(sita);
+        //计算交角章动，单位为0.0001角秒
+        //c+=(c1+c2*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;
+    //转换为弧度制
+    delta*=pi/(180*3600);
+    return delta;
+}
+
+// 获取地日运行参数，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))
+    {
+        //当分割符为/时，表明一组参数结束，将参数存入parameters中，并清空tmp
+        if(s[0]=='/')
+        {
+            parameters.push_back(tmp);
+            tmp.resize(0);
+            num_tmp=0;
+        }
+        //当分割符为=时，表明一组参数中的一类结束，将参数存入tmp中，并清空num_tmp
+        else if(s[0]=='=')
+        {
+            tmp.push_back(num_tmp);
+            num_tmp=0;
+        }
+        else
+        {
+            sscanf(s.c_str(),"%lf%lf%lf",&a,&b,&c);
+            num_tmp+=a*cos(b+c*t);
+        }
+    }
+    //关闭文件
+    fr.close();
+    /* 最终传参情况：
+     * parameters[0]为L系列，用来计算地球日心黄经
+     * parameters[1]为B系列，用来计算地球日心黄纬
+     * parameters[2]为R系列，用来计算地日距离
+     */
+    return parameters;
+}
+
+
+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+=nution(10*t);
+    //光行差修正，返回弧度制
+    L+=aberration(R);
+
+
+    //L转角度制
+    L*=180/pi;
+    while(L<0)
+        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 L=get_moon_longitude(l,t);
+//    double B=get_moon_latitude(b,t);
+//    double R=get_moon_R(r,t);
+
+//}