1 files changed, 424 insertions, 0 deletions

diff --git a/calendar.cpp b/calendar.cpp
new file mode 100644
index 0000000..d869d5d
--- /dev/null
+++ b/calendar.cpp
@@ -0,0 +1,424 @@
+// Encoding: UTF-8
+#include <iostream>
+#include <string>
+#include <vector>
+#include <fstream>
+#include <cmath>
+#include <ctime>
+using namespace std;
+
+double pi=3.14159265358979323846;
+double delta=5e-9;
+
+typedef struct tm Date;
+
+
+class Julian
+{
+    private:
+        double d[23][5]={
+            {-4000,  108371.700000,  -13036.800000,  392.000000,   0.000000},
+            { -500,   17201.000000,    -627.820000,   16.170000,  -0.341300},
+            { -150,   12200.600000,    -346.410000,    5.403000,  -0.159300},
+            {  150,    9113.800000,    -328.130000,   -1.647000,   0.037700},
+            {  500,    5707.500000,    -391.410000,    0.915000,   0.314500},
+            {  900,    2203.400000,    -283.450000,   13.034000,  -0.177800},
+            { 1300,     490.100000,     -57.350000,    2.085000,  -0.007200},
+            { 1600,     120.000000,      -9.810000,   -1.532000,   0.140300},
+            { 1700,      10.200000,      -0.910000,    0.510000,  -0.037000},
+            { 1800,      13.400000,      -0.720000,    0.202000,  -0.019300},
+            { 1830,       7.800000,      -1.810000,    0.416000,  -0.024700},
+            { 1860,       8.300000,      -0.130000,   -0.406000,   0.029200},
+            { 1880,      -5.400000,       0.320000,   -0.183000,   0.017300},
+            { 1900,      -2.300000,       2.060000,    0.169000,  -0.013500},
+            { 1920,      21.200000,       1.690000,   -0.304000,   0.016700},
+            { 1940,      24.200000,       1.220000,   -0.064000,   0.003100},
+            { 1960,      33.200000,       0.510000,    0.231000,  -0.010900},
+            { 1980,      51.000000,       1.290000,   -0.026000,   0.003200},
+            { 2000,      63.870000,       0.100000,    0.000000,   0.000000},
+            { 2005,      64.700000,       0.210000,    0.000000,   0.000000},
+            { 2012,      66.800000,       0.220000,    0.000000,   0.000000},
+            { 2018,      69.000000,       0.360000,    0.000000,   0.000000},
+            { 2028,      72.600000,       0.000000,    0.000000,   0.000000}
+        };
+
+        double dt_ext(int y,double jsd)
+        {
+            double dy=(double)(y-1820)/100.0;
+            return -20+jsd*dy*dy;
+        }
+
+        //计算力学时与世界时之差,传入年份
+        double delta_t(int y)
+        {
+            int y0=d[22][0];//表中最后一行的年份
+            double t0=d[22][1];//表中最后一行的delta_t
+            double jsd;
+
+            if(y>y0)
+            {
+                jsd=31;
+                if(y-y0>100)
+                {
+                    return dt_ext(y,jsd);
+                }
+
+                double v=dt_ext(y,jsd);
+                double dv=dt_ext(y0,jsd)-t0;
+                return v-dv*(y0+100-y)/100;
+            }
+
+            double res;
+            for(int i=0;i<22;i++)
+            {
+                if(y<d[i+1][0])
+                {
+                    break;
+                }
+                double t1=(y-d[i][0])/(d[i+1][0]-d[i][0])*10;
+                double t2=t1*t1;
+                double t3=t2*t1;
+                res=d[i][1]+d[i][2]*t1+d[i][3]*t2+d[i][4]*t3;
+            }
+            return res;
+
+
+        }
+    public:
+        // 计算儒略日
+        double getJulianDay(time_t time)
+        {
+            double t=(double)time;
+            return (double)t/86400.0+2440587.5;
+        }
+
+        // 计算儒略千年数
+        double getJulianKiloYear(time_t time)
+        {
+            double jd = getJulianDay(time);
+            return (double)(jd-2451545.0)/365250.0;
+        }
+
+        //儒略千年数转时间戳
+        time_t kiloYearToTime(double t,int year)
+        {
+            double jd=t*365250.0+2451545.0;
+            time_t time=(time_t)((jd-2440587.5)*86400.0-delta_t(year));
+            return time;
+        }
+};
+
+class parameter{
+    private:
+
+        // 计算地球日心黄经
+        double 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 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 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 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> 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 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 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>> 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;
+        }
+
+    public:
+
+        double 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;
+        }
+};
+
+void Plus(Date &date)
+{
+    if(date.tm_mday<10)
+    {
+        date.tm_mday=20;
+    }
+    else
+    {
+        date.tm_mday=6;
+        date.tm_mon++;
+    }
+    return;
+}
+
+double fuck(double a,double b)
+{
+    while(b<0)
+        b+=360;
+    while(b>=360)
+        b-=360;
+    if(b==0)
+    {
+        return min(fabs(b),fabs(360-b));
+    }
+    else
+    {
+        return fabs(a-b);
+    }
+}
+
+int main()
+{
+    Date* date=new Date;
+    cout<<"Please input the year: "<<endl;
+
+    //时间确认
+    cin>>date->tm_year;
+    date->tm_year-=1900;//年份减去1900
+                        //当年元旦零时
+    date->tm_mon=0;
+    date->tm_mday=1;
+    date->tm_hour=0;
+    date->tm_min=0;
+    date->tm_sec=0;
+
+    //输出年份
+    cout<<"Year: "<<date->tm_year+1900<<endl;
+
+    //计算当年二十四节气准确时刻
+    parameter p;
+    Julian julian;
+    time_t time;
+    //double dL;//dL为L对t导数
+    double t,t_end,t_middle;
+    double target;
+    char JieQi[24][10]={"小寒","大寒","立春","雨水","惊蛰","春分","清明","谷雨","立夏","小满","芒种","夏至","小暑","大暑","立秋","处暑","白露","秋分","寒露","霜降","立冬","小雪","大雪","冬至"};
+    //date->tm_mday=6;
+    for(int i=0;date->tm_mon<12;Plus(*date),i++)
+    {
+        //转换为时间戳
+        time=mktime(date);
+        //暂时放弃牛顿迭代法
+        //t=julian.getJulianKiloYear(time);
+        //while(fuck(p.longitude(t),(double)i*15-75)>=1e-8)
+        //{
+        //    dL=(p.longitude(t+delta)-p.longitude(t-delta))/(2*delta);
+        //    t=t-p.longitude(t)/dL;
+        //}
+
+
+        //计算儒略千年数
+        t=julian.getJulianKiloYear(time);
+        t_end=julian.getJulianKiloYear(time+86400);
+        target=i*15-75;
+        if(target<0)
+        {
+            target+=360;
+        }
+
+        //二分法
+        while(p.longitude(t)>target)
+        {
+            time-=86400;
+            t_end=t;
+            t=julian.getJulianKiloYear(time);
+        }
+        while(p.longitude(t_end)<=target)
+        {
+            time+=86400;
+            t=t_end;
+            t_end=julian.getJulianKiloYear(time);
+        }
+        while(t_end-t>=1e-10)
+        {
+            t_middle=(t+t_end)/2;
+            if(p.longitude(t_middle)>target)
+            {
+                t_end=t_middle;
+            }
+            else
+            {
+                t=t_middle;
+            }
+        }
+        //double t转为时间戳
+        time=julian.kiloYearToTime(t,date->tm_year+1900);
+
+        //date换新
+        date=localtime(&time);
+        printf("%d-%02d-%02d  %02d:%02d:%02d  %s\n",date->tm_year+1900,date->tm_mon+1,date->tm_mday,date->tm_hour,date->tm_min,date->tm_sec,JieQi[i]);
+    }
+    printf("\n\n");
+    return 0;
+}