From 9b3951efb56b3e1cf3b115c0a2fbfc3b90e28330 Mon Sep 17 00:00:00 2001 From: We-unite <3205135446@qq.com> Date: Thu, 7 Mar 2024 17:48:56 +0800 Subject: Input time, output sun longitude and latitude --- src/CMakeLists.txt | 2 +- src/List.cpp | 179 ---------------------------------------- src/calendar.cpp | 234 ++++++++--------------------------------------------- src/calendar.h | 32 +------- src/parameters.cpp | 44 +++++----- 5 files changed, 55 insertions(+), 436 deletions(-) delete mode 100644 src/List.cpp (limited to 'src') diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c48a058..53ddadc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -10,5 +10,5 @@ set(LIBRARY_OUTPUT_PATH ${PROJECT_SOURCE_DIR}/lib) # 将Julian.cpp与parameters.cpp编译为一个动态库 add_library(parameters SHARED Julian.cpp parameters.cpp) # 将calendar.cpp与List.cpp编译为可执行文件，链接动态库 -add_executable(calendar calendar.cpp List.cpp) +add_executable(calendar calendar.cpp) target_link_libraries(calendar PRIVATE parameters) diff --git a/src/List.cpp b/src/List.cpp deleted file mode 100644 index b4600bc..0000000 --- a/src/List.cpp +++ /dev/null @@ -1,179 +0,0 @@ -#include "calendar.h" - -List::List() { - // 为方便初始化，链表头结点设为空节点，注意不要使用 - head = new point(0, 0, 0, 0, 0, 0, false, false, false, -1, -1, false, 0); - tail = head; -} - -List::~List() { - point *p = head; - while (p != NULL) { - point *q = p->next; - delete p; - p = q; - } -} - -// 链表尾部插入 -void List::append(Date *date, bool isShuo, bool isJieqi, bool isZhongqi, - int JieqiIndex, double time) { - point *p = new point(date->tm_year + 1900, date->tm_mon + 1, date->tm_mday, - date->tm_hour, date->tm_min, date->tm_sec, isShuo, - isJieqi, isZhongqi, JieqiIndex, 0, false, time); - tail->next = p; - tail = p; - - if (isShuo) { - shuori_vec.push_back(p); - } else if (isJieqi && isZhongqi) { - zhongqi_vec.push_back(p); - } -} - -/* 对链表进行“排序” - * 先按照时间先后排序 - * 而后检查是否有中气和朔日在同一日但中气在前的情况 - * 如果出现，将中气人为地放在朔日之后，否则会影响置闰的判断 - * 这是今历沿用的顺治时宪历之规定， - * 即朔日只要和中气同日，不论先后，该中气即算作下一月 - */ -void List::lunar_sort() { - point *cur, *pre; - bool flag; - do { - flag = false; - pre = head; - do { - cur = pre->next; - if (later(cur, cur->next)) { - flag = true; - pre->next = cur->next; - cur->next = cur->next->next; - pre->next->next = cur; - } - pre = pre->next; - } while (pre->next != tail); - } while (flag); - - for (pre = head; pre->next != tail; pre = pre->next) { - cur = pre->next; - if (cur->isZhongqi && cur->next->isShuo && - cur->next->year == cur->year && cur->next->mon == cur->mon && - cur->next->day == cur->day) { - pre->next = cur->next; - cur->next = cur->next->next; - pre->next->next = cur; - } - } -} - -bool List::later(point *a, point *b) { - // 比较a是否比b晚 - // 比较年 - if (a->year > b->year) { - return true; - } else if (a->year < b->year) { - return false; - } - // 比较月 - else if (a->mon > b->mon) { - return true; - } else if (a->mon < b->mon) { - return false; - } - // 比较日 - else if (a->day > b->day) { - return true; - } else if (a->day < b->day) { - return false; - } - // 比较时 - else if (a->hour > b->hour) { - return true; - } else if (a->hour < b->hour) { - return false; - } - // 比较分 - else if (a->min > b->min) { - return true; - } else if (a->min < b->min) { - return false; - } - // 比较秒 - else if (a->sec > b->sec) { - return true; - } else if (a->sec < b->sec) { - return false; - } else { - return false; - } -} - -//置闰月 -void List::Run() { - point *cur; - int i = 0, j = 0; - int mon_index = 11; - bool flag; - while (later(zhongqi_vec[12], shuori_vec[i])) { - i++; - } - - if (i == 13) { - //冬至之间有13个朔日，置闰 - for (j = 0; j < shuori_vec.size() - 1; j++) { - flag = true; - cur = shuori_vec[j]->next; - while (cur != shuori_vec[j + 1]) { - if (cur->isZhongqi) { - //本月有中气，不置闰 - flag = false; - } - cur = cur->next; - } - if (flag) { - shuori_vec[j]->MonthIndex = mon_index; - shuori_vec[j]->RunYue = true; - j++; - break; - } else { - mon_index++; - if (mon_index == 13) { - mon_index = 1; - } - shuori_vec[j]->MonthIndex = mon_index; - shuori_vec[j]->RunYue = false; - } - } - } - //冬至之间有12个朔日，不置闰 - //同时处理置闰完毕之后的朔日 - for (; j < shuori_vec.size() - 1; j++) { - mon_index++; - if (mon_index == 13) { - mon_index = 1; - } - shuori_vec[j]->MonthIndex = mon_index; - shuori_vec[j]->RunYue = false; - } -} - -//输出链表内容 -void List::output() { - for (point *i = head->next; i != NULL; i = i->next) { - printf("%d-%02d-%02d %02d:%02d:%02d ", i->year, i->mon, i->day, - i->hour, i->min, i->sec); - if (i->isJieqi) { - printf("%s\n", jieqi[i->JieqiIndex]); - } else if (i->isShuo) { - printf("朔日"); - if (i->RunYue) { - printf(" 闰%d月", i->MonthIndex); - } else { - printf(" %d月", i->MonthIndex); - } - printf("\n"); - } - } -} diff --git a/src/calendar.cpp b/src/calendar.cpp index dcb96d7..9ad8413 100644 --- a/src/calendar.cpp +++ b/src/calendar.cpp @@ -2,226 +2,56 @@ #include "calendar.h" #include #include +#include /*以下为"calendar.h规定之量*/ -double pi = 3.14159265358979323846; double day = 86400; double delta = 1e-11; -char jieqi[25][10] = {"冬至", "小寒", "大寒", "立春", "雨水", "惊蛰", "春分", - "清明", "谷雨", "立夏", "小满", "芒种", "夏至", "小暑", - "大暑", "立秋", "处暑", "白露", "秋分", "寒露", "霜降", - "立冬", "小雪", "大雪", "冬至"}; /*计算需要用的类的实例化*/ Julian julian; parameter p; -List list; -Date *input(char *argv) { - Date *date = new Date; +void radToDMS(const char *message, int lowerBound, int upperBound, + double radians) { + // Convert radians to degrees + double degrees = radians * 180.0 / M_PI; - // 输入年份，自当年元旦起算 - date->tm_year = atoi(argv) - 1900; - date->tm_mon = 0; - date->tm_mday = 1; - date->tm_hour = 0; - date->tm_min = 0; - date->tm_sec = 0; + while (degrees < lowerBound) + degrees += (upperBound - lowerBound); + while (degrees >= upperBound) + degrees -= (upperBound - lowerBound); - return date; -} - -/*计算当年二十四节气时刻，自去年冬至到本年冬至*/ -void JieQi(Date *date) { - //计算起始时间为去年12月19日0时，开始寻找去年冬至 - time_t time = mktime(date); - time -= 13 * day; - - double t, t_end, t_middle; //采用二分法计算，t与t_end为间隔24h的两个时刻 - double target_angle; //目标角度 - Date *date_tmp; - - //计算每个节气对应太阳地心黄经 - for (int i = 0; i < 25; i++, time += 12 * day) { - target_angle = 15 * i + 270; - if (target_angle >= 360) - target_angle -= 360; - - t = julian.getJulianKiloYear(time); - t_end = julian.getJulianKiloYear(time + day); - - //定春分 - if (target_angle == 0) { - /* 春分太阳地心视黄经为0，一般二分法不便于计算 - * 因而限定范围为[355,360)∪(0,5]，即春分前后各5度 - * 当t_end时刻太阳地心视黄经在[355,360)时，证明春分在后边，t和t_end前进24h - */ - while (p.sun_longitude(t_end) < 360 && - p.sun_longitude(t_end) > 355) { - time += day; - t = t_end; - t_end = julian.getJulianKiloYear(time + day); - } - /* 循环体time每次在已经计算出的上一节气基础上前进12日，理论上春分在t_end之后 - * 当跳出上一while循环时理论上春分时刻正在t与t_end之间 - * 为保险起见，添加这个循环 - * 当t时刻太阳地心视黄经在[0,5]时，证明春分在前边，t和t_end后退24h - */ - while (p.sun_longitude(t) < 5) { - time -= day; - t_end = t; - t = julian.getJulianKiloYear(time); - } - - //二分法计算春分 - while (t_end - t > delta) { - t_middle = (t + t_end) / 2; - if (p.sun_longitude(t_middle) < 360 && - p.sun_longitude(t_middle) > 355) { - t = t_middle; - } else { - t_end = t_middle; - } - } - } - //其他情况 - else { - //当t_end时刻太阳地心视黄经小于目标角度，节气在后，t和t_end前进24h - while (p.sun_longitude(t_end) < target_angle) { - time += day; - t = t_end; - t_end = julian.getJulianKiloYear(time + day); - } + // Calculate degrees, minutes, and seconds + int deg = static_cast(degrees); + double min_sec = (degrees - deg) * 60.0; + int min = static_cast(min_sec); + double sec = (min_sec - min) * 60.0; - //当t时刻太阳地心视黄经大于目标角度，节气在前，t和t_end后退24h - //该循环也是一个保险循环 - while (p.sun_longitude(t) > target_angle) { - time -= day; - t_end = t; - t = julian.getJulianKiloYear(time); - } - - while (t_end - t > delta) { - t_middle = (t + t_end) / 2; - if (p.sun_longitude(t_middle) < target_angle) { - t = t_middle; - } else { - t_end = t_middle; - } - } - } - - //链表记录节气信息 - time = julian.kiloYearToTime(t, i == 0 ? date->tm_year + 1899 - : date->tm_year + 1900); - date_tmp = localtime(&time); - list.append(date_tmp, false, true, i % 2 == 0, i, t); - } -} - -// 用以判断朔日的函数，返回true则太阳地心黄经大于月球地心黄经，即朔日在此后 -bool date_forward(double sun_longitude, double moon_longitude) { - /* 判断逻辑 - * 当moon_longitudesun_longitude时， - * 如果moon_longitude-sun_longitude<20，即度数相差不大，可知朔日在此前 - * 否则，大于20时候，度数相差过大，理论上只可能为moon_longitude已接近360，sun_longitude略微大于0 - * 此时实际上是月亮应当追赶太阳，因而朔日在此后 - */ - return moon_longitude < sun_longitude - ? true - : (moon_longitude - sun_longitude >= 20); -} - -void ShuoRi(int year) { - double t, t_end, end, - t_middle; //采用二分法计算，t与t_end为间隔24h的两个时刻 - time_t time; - Date *date_tmp = new Date; - - //根据zhongqi_vec的最后一项（今年冬至）计算冬至下一日的0时刻 - date_tmp->tm_year = list.zhongqi_vec[12]->year - 1900; - date_tmp->tm_mon = list.zhongqi_vec[12]->mon - 1; - date_tmp->tm_mday = list.zhongqi_vec[12]->day + 1; - date_tmp->tm_hour = 0; - date_tmp->tm_min = 0; - date_tmp->tm_sec = 0; - date_tmp->tm_isdst = -1; - time = mktime(date_tmp); - end = julian.getJulianKiloYear(time); - - // 计算去年冬至时刻的time - date_tmp->tm_year = list.zhongqi_vec[0]->year - 1900; - date_tmp->tm_mon = list.zhongqi_vec[0]->mon - 1; - date_tmp->tm_mday = list.zhongqi_vec[0]->day; - date_tmp->tm_hour = list.zhongqi_vec[0]->hour; - date_tmp->tm_min = list.zhongqi_vec[0]->min; - date_tmp->tm_sec = list.zhongqi_vec[0]->sec; - date_tmp->tm_isdst = -1; - time = mktime(date_tmp); - - delete date_tmp; - - // 定朔日 - for (t = julian.getJulianKiloYear(time); t <= end; time += 25 * day) { - // 每月长度为28-30天，因而每次循环前进25天 - t = julian.getJulianKiloYear(time); - t_end = julian.getJulianKiloYear(time + day); - - // 太阳在月亮之前，月亮追赶，故朔日在后，t和t_end前进24h - while (date_forward(p.sun_longitude(t_end), p.moon_longitude(t_end))) { - time += day; - t = t_end; - t_end = julian.getJulianKiloYear(time + day); - } - // 太阳在月亮之后，月亮落后，故朔日在前，t和t_end后退24h（用以保险） - while (!date_forward(p.sun_longitude(t), p.moon_longitude(t))) { - time -= day; - t_end = t; - t = julian.getJulianKiloYear(time); - } - - while (t_end - t >= delta) { - t_middle = (t + t_end) / 2; - if (date_forward(p.sun_longitude(t_middle), - p.moon_longitude(t_middle))) { - t = t_middle; - } else { - t_end = t_middle; - } - } - // 输出朔日信息 - time = julian.kiloYearToTime(t, year - 1); - date_tmp = localtime(&time); - list.append(date_tmp, true, false, false, -1, t); - } + // Output the result + printf("%s: %d°%d'%06.3f''\n", message, deg, min, sec); } int main(int argc, char *argv[]) { - + Date date; if (argc != 2) { - printf("Usage: %s year\\n\n", argv[0]); - return 0; + printf("Input the time you want to calculate in " + "format: "); + scanf("%d-%d-%d,%d:%d:%d", &date.tm_year, &date.tm_mon, &date.tm_mday, + &date.tm_hour, &date.tm_min, &date.tm_sec); + } else { + sscanf(argv[1], "%d-%d-%d,%d:%d:%d", &date.tm_year, &date.tm_mon, + &date.tm_mday, &date.tm_hour, &date.tm_min, &date.tm_sec); } - Date *date = input(argv[1]); - int year = date->tm_year + 1900; - double begin, end; - - JieQi(date); - ShuoRi(year); - - list.lunar_sort(); - - // 奉旨置闰 - list.Run(); - - // 输出 - list.output(); - - // 释放内存，使用List的析构函数 + date.tm_year -= 1900; + date.tm_mon -= 1; + date.tm_isdst = -1; - delete date; + time_t time = mktime(&date); + double t = julian.getJulianKiloYear(time); + pair sun = p.sun_longitude(t); + radToDMS("Sun longitude: ", 0, 360, sun.first); + radToDMS("Sun latitude: ", -90, 90, sun.second); return 0; } diff --git a/src/calendar.h b/src/calendar.h index 73676f5..ef785c0 100644 --- a/src/calendar.h +++ b/src/calendar.h @@ -9,7 +9,6 @@ #include using namespace std; -extern double pi; extern double delta; extern char jieqi[25][10]; @@ -67,7 +66,7 @@ class parameter { // 计算地球日心黄经 double get_longitude(vector l, double t); - // 计算地球日心黄纬 + // 计算太阳地心黄纬 double get_latitude(vector b, double t); double get_R(vector r, double t); @@ -88,36 +87,9 @@ class parameter { vector> get_parameters(double t); public: - double sun_longitude(double t); + pair sun_longitude(double t); double moon_longitude(double t); }; -class List { - private: - point *head, *tail; //链表头尾指针 - - bool later(point *a, point *b); - - public: - vector shuori_vec; //记录各朔日的地址 - vector zhongqi_vec; //记录各中气的地址 - - List(); - ~List(); - - // 链表尾部插入 - void append(Date *date, bool isShuo, bool isJieqi, bool isZhongqi, - int JieqiIndex, double time); - - //对链表进行“排序” - void lunar_sort(); - - //置闰月 - void Run(); - - //输出链表内容 - void output(); -}; - #endif diff --git a/src/parameters.cpp b/src/parameters.cpp index 57d98c2..2e2758a 100644 --- a/src/parameters.cpp +++ b/src/parameters.cpp @@ -1,5 +1,6 @@ -#include "parameters.h" #include "calendar.h" +#include "parameters.h" +#include using namespace std; //凡无特殊说明的函数，全部返回弧度制 @@ -11,7 +12,7 @@ double parameter::get_longitude(vector l, double t) { L += l[i] * pow(t, i); } //返回弧度制，太阳地心黄经=地球日心黄经+180° - return L + pi; + return L + M_PI; } // 计算地球日心黄纬 @@ -36,19 +37,19 @@ double parameter::get_R(vector r, double t) { // 转换FK5误差 double parameter::delta_FK5(double L, double B, double T) { // L转角度制 - L *= 180 / pi; + L *= 180 / M_PI; //计算L'，角度制 double L_p = L - 1.397 * T - 0.00031 * T * T; //计算L'，弧度制 - L_p *= pi / 180; + L_p *= M_PI / 180; while (L_p < 0) - L_p += 2 * pi; - while (L_p >= 2 * pi) - L_p -= 2 * pi; + L_p += 2 * M_PI; + while (L_p >= 2 * M_PI) + L_p -= 2 * M_PI; //计算delta_L，单位为角秒 double delta_L = -0.09033 + 0.03916 * (cos(L_p) + sin(L_p)) * tan(B); //转换为弧度制 - delta_L *= pi / (180 * 3600); + delta_L *= M_PI / (180 * 3600); return delta_L; } @@ -88,14 +89,14 @@ double parameter::nutation(double T) { nutation_parameters[i][2] * ang[2] + nutation_parameters[i][3] * ang[3] + nutation_parameters[i][4] * ang[4]; - sita *= pi / 180; + sita *= M_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 *= M_PI / (180 * 3600); return s; } @@ -104,7 +105,7 @@ double parameter::aberration(double R) { //公式，单位0.0001角秒 double delta = -20.4898 / R; //转换为弧度制 - delta *= pi / (180 * 3600); + delta *= M_PI / (180 * 3600); return delta; } @@ -157,7 +158,7 @@ vector> parameter::get_parameters(double t) { } // 计算太阳地心黄经，外部调用，返回角度制 -double parameter::sun_longitude(double t) { +pair parameter::sun_longitude(double t) { //从文件中计算地日运行参数 vector> p = get_parameters(t); @@ -177,13 +178,8 @@ double parameter::sun_longitude(double t) { //光行差修正，返回弧度制 L += aberration(R); - // L转角度制 - L *= 180 / pi; - while (L < 0) - L += 360; - while (L >= 360) - L -= 360; - return L; + // 返回太阳地心黄经和黄纬 + return make_pair(L, -B); } // 计算月球地心黄经，外部调用，返回角度制 @@ -202,7 +198,7 @@ double parameter::moon_longitude(double t) { moon_longitude_parameters[i][1] * angles[1] + moon_longitude_parameters[i][2] * angles[2] + moon_longitude_parameters[i][3] * angles[3]; - sita *= pi / 180; + sita *= M_PI / 180; EI += moon_longitude_parameters[i][4] * sin(sita) * pow(E, fabs(moon_longitude_parameters[i][1])); } @@ -212,12 +208,12 @@ double parameter::moon_longitude(double t) { 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 * M_PI / 180); + EI += 1962.0 * sin((L - angles[3]) * M_PI / 180); + EI += 318.0 * sin(a2 * M_PI / 180); L += EI / 1000000; - L += 180 / pi * nutation(T); //地球章动修正 + L += 180 / M_PI * nutation(T); //地球章动修正 while (L < 0) L += 360; -- cgit v1.2.3-70-g09d2