From 27772ba64bb8f49a9d9057acaed2e08d0f7cec68 Mon Sep 17 00:00:00 2001 From: We-unite <3205135446@qq.com> Date: Fri, 15 Sep 2023 23:21:24 +0800 Subject: Formalized this program. - Make src dir, put all code into it. - Devide old source code into several parts of files with makefile. - Delete some useless files. --- src/Julian.cpp | 89 +++++++++++++++++++++++ src/calendar.cpp | 120 ++++++++++++++++++++++++++++++ src/calendar.h | 71 ++++++++++++++++++ src/parameters.cpp | 210 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 490 insertions(+) create mode 100644 src/Julian.cpp create mode 100644 src/calendar.cpp create mode 100644 src/calendar.h create mode 100644 src/parameters.cpp (limited to 'src') diff --git a/src/Julian.cpp b/src/Julian.cpp new file mode 100644 index 0000000..1d3c4b7 --- /dev/null +++ b/src/Julian.cpp @@ -0,0 +1,89 @@ +#include "calendar.h" + +double Julian::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 Julian::dt_ext(int y,double jsd) +{ + double dy=(double)(y-1820)/100.0; + return -20+jsd*dy*dy; +} + +//计算力学时与世界时之差,传入年份 +double Julian::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=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: "<>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: "<tm_year+1900<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.sun_longitude(t)>target) + { + time-=86400; + t_end=t; + t=julian.getJulianKiloYear(time); + } + while(p.sun_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.sun_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; +} diff --git a/src/calendar.h b/src/calendar.h new file mode 100644 index 0000000..c87531d --- /dev/null +++ b/src/calendar.h @@ -0,0 +1,71 @@ +#ifndef _CALENDAR_H_ +#define _CALENDAR_H_ + + +#include +#include +#include +#include +#include +#include +using namespace std; + +typedef struct tm Date; + +extern double pi; +extern double delta; + +class Julian +{ + private: + static double d[23][5]; + + double dt_ext(int y,double jsd); + + //计算力学时与世界时之差,传入年份 + double delta_t(int y); + public: + // 计算儒略日 + double getJulianDay(time_t time); + + // 计算儒略千年数 + double getJulianKiloYear(time_t time); + + //儒略千年数转时间戳 + time_t kiloYearToTime(double t,int year); +}; + + +class parameter{ + private: + + // 计算地球日心黄经 + double get_longitude(vector l,double t); + + // 计算地球日心黄纬 + double get_latitude(vector b,double t); + + double get_R(vector r,double t); + + // 转换FK5误差，返回弧度制 + double delta_FK5(double L,double B,double T); + + //获取章动有关角 + vector get_angle(double T); + + + //章动修正 + double nution(double T); + + //光行差修正，返回弧度制 + double aberration(double R); + + // 获取地日运行参数，L为地球日心黄经，B为地球日心黄纬，R为地日距离 + vector> get_parameters(double t); + + public: + + double sun_longitude(double t); +}; + +#endif 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 l,double t) +{ + double L=l[0]; + for(int i=1;i b,double t) +{ + double B=b[0]; + for(int i=1;i r,double t) +{ + double R=r[0]; + for(int i=1;i2*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 parameter::get_angle(double T) +{ + vector 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 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> 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> parameters; + vector 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> p=get_parameters(t); + + //取出各类参数 + vector 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> p=get_moon_parameters(t); +// vector 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); + +//} -- cgit v1.2.3-70-g09d2