summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorWe-unite <3205135446@qq.com>2024-03-07 17:48:56 +0800
committerWe-unite <3205135446@qq.com>2024-03-07 17:48:56 +0800
commit9b3951efb56b3e1cf3b115c0a2fbfc3b90e28330 (patch)
tree118d5bfd8e707f4713b9970b265b37d351664cd0
parent7e85ab6f9e7114ba05264005172e41fb1c47532c (diff)
downloadcalendar-9b3951efb56b3e1cf3b115c0a2fbfc3b90e28330.tar.gz
calendar-9b3951efb56b3e1cf3b115c0a2fbfc3b90e28330.zip
Input time, output sun longitude and latitude
Diffstat
-rw-r--r--src/CMakeLists.txt2
-rw-r--r--src/List.cpp179
-rw-r--r--src/calendar.cpp234
-rw-r--r--src/calendar.h32
-rw-r--r--src/parameters.cpp44
5 files changed, 55 insertions, 436 deletions
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)
10# 将Julian.cpp与parameters.cpp编译为一个动态库 10# 将Julian.cpp与parameters.cpp编译为一个动态库
11add_library(parameters SHARED Julian.cpp parameters.cpp) 11add_library(parameters SHARED Julian.cpp parameters.cpp)
12# 将calendar.cpp与List.cpp编译为可执行文件，链接动态库 12# 将calendar.cpp与List.cpp编译为可执行文件，链接动态库
13add_executable(calendar calendar.cpp List.cpp) 13add_executable(calendar calendar.cpp)
14target_link_libraries(calendar PRIVATE parameters) 14target_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 @@
1#include "calendar.h"
2
3List::List() {
4 // 为方便初始化，链表头结点设为空节点，注意不要使用
5 head = new point(0, 0, 0, 0, 0, 0, false, false, false, -1, -1, false, 0);
6 tail = head;
7}
8
9List::~List() {
10 point *p = head;
11 while (p != NULL) {
12 point *q = p->next;
13 delete p;
14 p = q;
15 }
16}
17
18// 链表尾部插入
19void List::append(Date *date, bool isShuo, bool isJieqi, bool isZhongqi,
20 int JieqiIndex, double time) {
21 point *p = new point(date->tm_year + 1900, date->tm_mon + 1, date->tm_mday,
22 date->tm_hour, date->tm_min, date->tm_sec, isShuo,
23 isJieqi, isZhongqi, JieqiIndex, 0, false, time);
24 tail->next = p;
25 tail = p;
26
27 if (isShuo) {
28 shuori_vec.push_back(p);
29 } else if (isJieqi && isZhongqi) {
30 zhongqi_vec.push_back(p);
31 }
32}
33
34/* 对链表进行“排序”
35 * 先按照时间先后排序
36 * 而后检查是否有中气和朔日在同一日但中气在前的情况
37 * 如果出现，将中气人为地放在朔日之后，否则会影响置闰的判断
38 * 这是今历沿用的顺治时宪历之规定，
39 * 即朔日只要和中气同日，不论先后，该中气即算作下一月
40 */
41void List::lunar_sort() {
42 point *cur, *pre;
43 bool flag;
44 do {
45 flag = false;
46 pre = head;
47 do {
48 cur = pre->next;
49 if (later(cur, cur->next)) {
50 flag = true;
51 pre->next = cur->next;
52 cur->next = cur->next->next;
53 pre->next->next = cur;
54 }
55 pre = pre->next;
56 } while (pre->next != tail);
57 } while (flag);
58
59 for (pre = head; pre->next != tail; pre = pre->next) {
60 cur = pre->next;
61 if (cur->isZhongqi && cur->next->isShuo &&
62 cur->next->year == cur->year && cur->next->mon == cur->mon &&
63 cur->next->day == cur->day) {
64 pre->next = cur->next;
65 cur->next = cur->next->next;
66 pre->next->next = cur;
67 }
68 }
69}
70
71bool List::later(point *a, point *b) {
72 // 比较a是否比b晚
73 // 比较年
74 if (a->year > b->year) {
75 return true;
76 } else if (a->year < b->year) {
77 return false;
78 }
79 // 比较月
80 else if (a->mon > b->mon) {
81 return true;
82 } else if (a->mon < b->mon) {
83 return false;
84 }
85 // 比较日
86 else if (a->day > b->day) {
87 return true;
88 } else if (a->day < b->day) {
89 return false;
90 }
91 // 比较时
92 else if (a->hour > b->hour) {
93 return true;
94 } else if (a->hour < b->hour) {
95 return false;
96 }
97 // 比较分
98 else if (a->min > b->min) {
99 return true;
100 } else if (a->min < b->min) {
101 return false;
102 }
103 // 比较秒
104 else if (a->sec > b->sec) {
105 return true;
106 } else if (a->sec < b->sec) {
107 return false;
108 } else {
109 return false;
110 }
111}
112
113//置闰月
114void List::Run() {
115 point *cur;
116 int i = 0, j = 0;
117 int mon_index = 11;
118 bool flag;
119 while (later(zhongqi_vec[12], shuori_vec[i])) {
120 i++;
121 }
122
123 if (i == 13) {
124 //冬至之间有13个朔日，置闰
125 for (j = 0; j < shuori_vec.size() - 1; j++) {
126 flag = true;
127 cur = shuori_vec[j]->next;
128 while (cur != shuori_vec[j + 1]) {
129 if (cur->isZhongqi) {
130 //本月有中气，不置闰
131 flag = false;
132 }
133 cur = cur->next;
134 }
135 if (flag) {
136 shuori_vec[j]->MonthIndex = mon_index;
137 shuori_vec[j]->RunYue = true;
138 j++;
139 break;
140 } else {
141 mon_index++;
142 if (mon_index == 13) {
143 mon_index = 1;
144 }
145 shuori_vec[j]->MonthIndex = mon_index;
146 shuori_vec[j]->RunYue = false;
147 }
148 }
149 }
150 //冬至之间有12个朔日，不置闰
151 //同时处理置闰完毕之后的朔日
152 for (; j < shuori_vec.size() - 1; j++) {
153 mon_index++;
154 if (mon_index == 13) {
155 mon_index = 1;
156 }
157 shuori_vec[j]->MonthIndex = mon_index;
158 shuori_vec[j]->RunYue = false;
159 }
160}
161
162//输出链表内容
163void List::output() {
164 for (point *i = head->next; i != NULL; i = i->next) {
165 printf("%d-%02d-%02d %02d:%02d:%02d ", i->year, i->mon, i->day,
166 i->hour, i->min, i->sec);
167 if (i->isJieqi) {
168 printf("%s\n", jieqi[i->JieqiIndex]);
169 } else if (i->isShuo) {
170 printf("朔日");
171 if (i->RunYue) {
172 printf(" 闰%d月", i->MonthIndex);
173 } else {
174 printf(" %d月", i->MonthIndex);
175 }
176 printf("\n");
177 }
178 }
179}
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 @@
2#include "calendar.h" 2#include "calendar.h"
3#include <cstdio> 3#include <cstdio>
4#include <ctime> 4#include <ctime>
5#include <iomanip>
5 6
6/*以下为"calendar.h规定之量*/ 7/*以下为"calendar.h规定之量*/
7double pi = 3.14159265358979323846;
8double day = 86400; 8double day = 86400;
9double delta = 1e-11; 9double delta = 1e-11;
10char jieqi[25][10] = {"冬至", "小寒", "大寒", "立春", "雨水", "惊蛰", "春分",
11 "清明", "谷雨", "立夏", "小满", "芒种", "夏至", "小暑",
12 "大暑", "立秋", "处暑", "白露", "秋分", "寒露", "霜降",
13 "立冬", "小雪", "大雪", "冬至"};
14 10
15/*计算需要用的类的实例化*/ 11/*计算需要用的类的实例化*/
16Julian julian; 12Julian julian;
17parameter p; 13parameter p;
18List list;
19 14
20Date *input(char *argv) { 15void radToDMS(const char *message, int lowerBound, int upperBound,
21 Date *date = new Date; 16 double radians) {
17 // Convert radians to degrees
18 double degrees = radians * 180.0 / M_PI;
22 19
23 // 输入年份，自当年元旦起算 20 while (degrees < lowerBound)
24 date->tm_year = atoi(argv) - 1900; 21 degrees += (upperBound - lowerBound);
25 date->tm_mon = 0; 22 while (degrees >= upperBound)
26 date->tm_mday = 1; 23 degrees -= (upperBound - lowerBound);
27 date->tm_hour = 0;
28 date->tm_min = 0;
29 date->tm_sec = 0;
30 24
31 return date; 25 // Calculate degrees, minutes, and seconds
32} 26 int deg = static_cast<int>(degrees);
33 27 double min_sec = (degrees - deg) * 60.0;
34/*计算当年二十四节气时刻，自去年冬至到本年冬至*/ 28 int min = static_cast<int>(min_sec);
35void JieQi(Date *date) { 29 double sec = (min_sec - min) * 60.0;
36 //计算起始时间为去年12月19日0时，开始寻找去年冬至
37 time_t time = mktime(date);
38 time -= 13 * day;
39
40 double t, t_end, t_middle; //采用二分法计算，t与t_end为间隔24h的两个时刻
41 double target_angle; //目标角度
42 Date *date_tmp;
43
44 //计算每个节气对应太阳地心黄经
45 for (int i = 0; i < 25; i++, time += 12 * day) {
46 target_angle = 15 * i + 270;
47 if (target_angle >= 360)
48 target_angle -= 360;
49
50 t = julian.getJulianKiloYear(time);
51 t_end = julian.getJulianKiloYear(time + day);
52
53 //定春分
54 if (target_angle == 0) {
55 /* 春分太阳地心视黄经为0，一般二分法不便于计算
56 * 因而限定范围为[355,360)∪(0,5]，即春分前后各5度
57 * 当t_end时刻太阳地心视黄经在[355,360)时，证明春分在后边，t和t_end前进24h
58 */
59 while (p.sun_longitude(t_end) < 360 &&
60 p.sun_longitude(t_end) > 355) {
61 time += day;
62 t = t_end;
63 t_end = julian.getJulianKiloYear(time + day);
64 }
65 /* 循环体time每次在已经计算出的上一节气基础上前进12日，理论上春分在t_end之后
66 * 当跳出上一while循环时理论上春分时刻正在t与t_end之间
67 * 为保险起见，添加这个循环
68 * 当t时刻太阳地心视黄经在[0,5]时，证明春分在前边，t和t_end后退24h
69 */
70 while (p.sun_longitude(t) < 5) {
71 time -= day;
72 t_end = t;
73 t = julian.getJulianKiloYear(time);
74 }
75
76 //二分法计算春分
77 while (t_end - t > delta) {
78 t_middle = (t + t_end) / 2;
79 if (p.sun_longitude(t_middle) < 360 &&
80 p.sun_longitude(t_middle) > 355) {
81 t = t_middle;
82 } else {
83 t_end = t_middle;
84 }
85 }
86 }
87 //其他情况
88 else {
89 //当t_end时刻太阳地心视黄经小于目标角度，节气在后，t和t_end前进24h
90 while (p.sun_longitude(t_end) < target_angle) {
91 time += day;
92 t = t_end;
93 t_end = julian.getJulianKiloYear(time + day);
94 }
95 30
96 //当t时刻太阳地心视黄经大于目标角度，节气在前，t和t_end后退24h 31 // Output the result
97 //该循环也是一个保险循环 32 printf("%s: %d°%d'%06.3f''\n", message, deg, min, sec);
98 while (p.sun_longitude(t) > target_angle) {
99 time -= day;
100 t_end = t;
101 t = julian.getJulianKiloYear(time);
102 }
103
104 while (t_end - t > delta) {
105 t_middle = (t + t_end) / 2;
106 if (p.sun_longitude(t_middle) < target_angle) {
107 t = t_middle;
108 } else {
109 t_end = t_middle;
110 }
111 }
112 }
113
114 //链表记录节气信息
115 time = julian.kiloYearToTime(t, i == 0 ? date->tm_year + 1899
116 : date->tm_year + 1900);
117 date_tmp = localtime(&time);
118 list.append(date_tmp, false, true, i % 2 == 0, i, t);
119 }
120}
121
122// 用以判断朔日的函数，返回true则太阳地心黄经大于月球地心黄经，即朔日在此后
123bool date_forward(double sun_longitude, double moon_longitude) {
124 /* 判断逻辑
125 * 当moon_longitude<sun_longitude时，由于时刻已经在朔日附近
126 * 因而可以判断朔日在该时刻之后
127 * 反之，当moon_longitude>sun_longitude时，
128 * 如果moon_longitude-sun_longitude<20，即度数相差不大，可知朔日在此前
129 * 否则，大于20时候，度数相差过大，理论上只可能为moon_longitude已接近360，sun_longitude略微大于0
130 * 此时实际上是月亮应当追赶太阳，因而朔日在此后
131 */
132 return moon_longitude < sun_longitude
133 ? true
134 : (moon_longitude - sun_longitude >= 20);
135}
136
137void ShuoRi(int year) {
138 double t, t_end, end,
139 t_middle; //采用二分法计算，t与t_end为间隔24h的两个时刻
140 time_t time;
141 Date *date_tmp = new Date;
142
143 //根据zhongqi_vec的最后一项（今年冬至）计算冬至下一日的0时刻
144 date_tmp->tm_year = list.zhongqi_vec[12]->year - 1900;
145 date_tmp->tm_mon = list.zhongqi_vec[12]->mon - 1;
146 date_tmp->tm_mday = list.zhongqi_vec[12]->day + 1;
147 date_tmp->tm_hour = 0;
148 date_tmp->tm_min = 0;
149 date_tmp->tm_sec = 0;
150 date_tmp->tm_isdst = -1;
151 time = mktime(date_tmp);
152 end = julian.getJulianKiloYear(time);
153
154 // 计算去年冬至时刻的time
155 date_tmp->tm_year = list.zhongqi_vec[0]->year - 1900;
156 date_tmp->tm_mon = list.zhongqi_vec[0]->mon - 1;
157 date_tmp->tm_mday = list.zhongqi_vec[0]->day;
158 date_tmp->tm_hour = list.zhongqi_vec[0]->hour;
159 date_tmp->tm_min = list.zhongqi_vec[0]->min;
160 date_tmp->tm_sec = list.zhongqi_vec[0]->sec;
161 date_tmp->tm_isdst = -1;
162 time = mktime(date_tmp);
163
164 delete date_tmp;
165
166 // 定朔日
167 for (t = julian.getJulianKiloYear(time); t <= end; time += 25 * day) {
168 // 每月长度为28-30天，因而每次循环前进25天
169 t = julian.getJulianKiloYear(time);
170 t_end = julian.getJulianKiloYear(time + day);
171
172 // 太阳在月亮之前，月亮追赶，故朔日在后，t和t_end前进24h
173 while (date_forward(p.sun_longitude(t_end), p.moon_longitude(t_end))) {
174 time += day;
175 t = t_end;
176 t_end = julian.getJulianKiloYear(time + day);
177 }
178 // 太阳在月亮之后，月亮落后，故朔日在前，t和t_end后退24h（用以保险）
179 while (!date_forward(p.sun_longitude(t), p.moon_longitude(t))) {
180 time -= day;
181 t_end = t;
182 t = julian.getJulianKiloYear(time);
183 }
184
185 while (t_end - t >= delta) {
186 t_middle = (t + t_end) / 2;
187 if (date_forward(p.sun_longitude(t_middle),
188 p.moon_longitude(t_middle))) {
189 t = t_middle;
190 } else {
191 t_end = t_middle;
192 }
193 }
194 // 输出朔日信息
195 time = julian.kiloYearToTime(t, year - 1);
196 date_tmp = localtime(&time);
197 list.append(date_tmp, true, false, false, -1, t);
198 }
199} 33}
200 34
201int main(int argc, char *argv[]) { 35int main(int argc, char *argv[]) {
202 36 Date date;
203 if (argc != 2) { 37 if (argc != 2) {
204 printf("Usage: %s year\\n\n", argv[0]); 38 printf("Input the time you want to calculate in <YYYY-MM-DD,HH:MM:SS> "
205 return 0; 39 "format: ");
40 scanf("%d-%d-%d,%d:%d:%d", &date.tm_year, &date.tm_mon, &date.tm_mday,
41 &date.tm_hour, &date.tm_min, &date.tm_sec);
42 } else {
43 sscanf(argv[1], "%d-%d-%d,%d:%d:%d", &date.tm_year, &date.tm_mon,
44 &date.tm_mday, &date.tm_hour, &date.tm_min, &date.tm_sec);
206 } 45 }
207 Date *date = input(argv[1]); 46 date.tm_year -= 1900;
208 int year = date->tm_year + 1900; 47 date.tm_mon -= 1;
209 double begin, end; 48 date.tm_isdst = -1;
210
211 JieQi(date);
212 ShuoRi(year);
213
214 list.lunar_sort();
215
216 // 奉旨置闰
217 list.Run();
218
219 // 输出
220 list.output();
221
222 // 释放内存，使用List的析构函数
223 49
224 delete date; 50 time_t time = mktime(&date);
51 double t = julian.getJulianKiloYear(time);
225 52
53 pair<double, double> sun = p.sun_longitude(t);
54 radToDMS("Sun longitude: ", 0, 360, sun.first);
55 radToDMS("Sun latitude: ", -90, 90, sun.second);
226 return 0; 56 return 0;
227} 57}
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 @@
9#include <vector> 9#include <vector>
10using namespace std; 10using namespace std;
11 11
12extern double pi;
13extern double delta; 12extern double delta;
14extern char jieqi[25][10]; 13extern char jieqi[25][10];
15 14
@@ -67,7 +66,7 @@ class parameter {
67 // 计算地球日心黄经 66 // 计算地球日心黄经
68 double get_longitude(vector<double> l, double t); 67 double get_longitude(vector<double> l, double t);
69 68
70 // 计算地球日心黄纬 69 // 计算阳心黄纬
71 double get_latitude(vector<double> b, double t); 70 double get_latitude(vector<double> b, double t);
72 71
73 double get_R(vector<double> r, double t); 72 double get_R(vector<double> r, double t);
@@ -88,36 +87,9 @@ class parameter {
88 vector<vector<double>> get_parameters(double t); 87 vector<vector<double>> get_parameters(double t);
89 88
90 public: 89 public:
91 double sun_longitude(double t); 90 pair<double, double> sun_longitude(double t);
92 91
93 double moon_longitude(double t); 92 double moon_longitude(double t);
94}; 93};
95 94
96class List {
97 private:
98 point *head, *tail; //链表头尾指针
99
100 bool later(point *a, point *b);
101
102 public:
103 vector<point *> shuori_vec; //记录各朔日的地址
104 vector<point *> zhongqi_vec; //记录各中气的地址
105
106 List();
107 ~List();
108
109 // 链表尾部插入
110 void append(Date *date, bool isShuo, bool isJieqi, bool isZhongqi,
111 int JieqiIndex, double time);
112
113 //对链表进行“排序”
114 void lunar_sort();
115
116 //置闰月
117 void Run();
118
119 //输出链表内容
120 void output();
121};
122
123#endif 95#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 @@
1#include "parameters.h"
2#include "calendar.h" 1#include "calendar.h"
2#include "parameters.h"
3#include <utility>
3using namespace std; 4using namespace std;
4 5
5//凡无特殊说明的函数，全部返回弧度制 6//凡无特殊说明的函数，全部返回弧度制
@@ -11,7 +12,7 @@ double parameter::get_longitude(vector<double> l, double t) {
11 L += l[i] * pow(t, i); 12 L += l[i] * pow(t, i);
12 } 13 }
13 //返回弧度制，太阳地心黄经=地球日心黄经+180° 14 //返回弧度制，太阳地心黄经=地球日心黄经+180°
14 return L + pi; 15 return L + M_PI;
15} 16}
16 17
17// 计算地球日心黄纬 18// 计算地球日心黄纬
@@ -36,19 +37,19 @@ double parameter::get_R(vector<double> r, double t) {
36// 转换FK5误差 37// 转换FK5误差
37double parameter::delta_FK5(double L, double B, double T) { 38double parameter::delta_FK5(double L, double B, double T) {
38 // L转角度制 39 // L转角度制
39 L *= 180 / pi; 40 L *= 180 / M_PI;
40 //计算L'，角度制 41 //计算L'，角度制
41 double L_p = L - 1.397 * T - 0.00031 * T * T; 42 double L_p = L - 1.397 * T - 0.00031 * T * T;
42 //计算L'，弧度制 43 //计算L'，弧度制
43 L_p *= pi / 180; 44 L_p *= M_PI / 180;
44 while (L_p < 0) 45 while (L_p < 0)
45 L_p += 2 * pi; 46 L_p += 2 * M_PI;
46 while (L_p >= 2 * pi) 47 while (L_p >= 2 * M_PI)
47 L_p -= 2 * pi; 48 L_p -= 2 * M_PI;
48 //计算delta_L，单位为角秒 49 //计算delta_L，单位为角秒
49 double delta_L = -0.09033 + 0.03916 * (cos(L_p) + sin(L_p)) * tan(B); 50 double delta_L = -0.09033 + 0.03916 * (cos(L_p) + sin(L_p)) * tan(B);
50 //转换为弧度制 51 //转换为弧度制
51 delta_L *= pi / (180 * 3600); 52 delta_L *= M_PI / (180 * 3600);
52 return delta_L; 53 return delta_L;
53} 54}
54 55
@@ -88,14 +89,14 @@ double parameter::nutation(double T) {
88 nutation_parameters[i][2] * ang[2] + 89 nutation_parameters[i][2] * ang[2] +
89 nutation_parameters[i][3] * ang[3] + 90 nutation_parameters[i][3] * ang[3] +
90 nutation_parameters[i][4] * ang[4]; 91 nutation_parameters[i][4] * ang[4];
91 sita *= pi / 180; 92 sita *= M_PI / 180;
92 s += (nutation_parameters[i][5] + nutation_parameters[i][6] * T) * 93 s += (nutation_parameters[i][5] + nutation_parameters[i][6] * T) *
93 sin(sita); 94 sin(sita);
94 // c+=(nutation_parameters[i][7]+nutation_parameters[i][8]*T)*cos(sita); 95 // c+=(nutation_parameters[i][7]+nutation_parameters[i][8]*T)*cos(sita);
95 } 96 }
96 // 计算得到的s和c单位0.0001秒，转换为弧度制 97 // 计算得到的s和c单位0.0001秒，转换为弧度制
97 s *= 0.0001; 98 s *= 0.0001;
98 s *= pi / (180 * 3600); 99 s *= M_PI / (180 * 3600);
99 return s; 100 return s;
100} 101}
101 102
@@ -104,7 +105,7 @@ double parameter::aberration(double R) {
104 //公式，单位0.0001角秒 105 //公式，单位0.0001角秒
105 double delta = -20.4898 / R; 106 double delta = -20.4898 / R;
106 //转换为弧度制 107 //转换为弧度制
107 delta *= pi / (180 * 3600); 108 delta *= M_PI / (180 * 3600);
108 return delta; 109 return delta;
109} 110}
110 111
@@ -157,7 +158,7 @@ vector<vector<double>> parameter::get_parameters(double t) {
157} 158}
158 159
159// 计算太阳地心黄经，外部调用，返回角度制 160// 计算太阳地心黄经，外部调用，返回角度制
160double parameter::sun_longitude(double t) { 161pair<double, double> parameter::sun_longitude(double t) {
161 //从文件中计算地日运行参数 162 //从文件中计算地日运行参数
162 vector<vector<double>> p = get_parameters(t); 163 vector<vector<double>> p = get_parameters(t);
163 164
@@ -177,13 +178,8 @@ double parameter::sun_longitude(double t) {
177 //光行差修正，返回弧度制 178 //光行差修正，返回弧度制
178 L += aberration(R); 179 L += aberration(R);
179 180
180 // L转角度制 181 // 返回太阳地心黄经和黄纬
181 L *= 180 / pi; 182 return make_pair(L, -B);
182 while (L < 0)
183 L += 360;
184 while (L >= 360)
185 L -= 360;
186 return L;
187} 183}
188 184
189// 计算月球地心黄经，外部调用，返回角度制 185// 计算月球地心黄经，外部调用，返回角度制
@@ -202,7 +198,7 @@ double parameter::moon_longitude(double t) {
202 moon_longitude_parameters[i][1] * angles[1] + 198 moon_longitude_parameters[i][1] * angles[1] +
203 moon_longitude_parameters[i][2] * angles[2] + 199 moon_longitude_parameters[i][2] * angles[2] +
204 moon_longitude_parameters[i][3] * angles[3]; 200 moon_longitude_parameters[i][3] * angles[3];
205 sita *= pi / 180; 201 sita *= M_PI / 180;
206 EI += moon_longitude_parameters[i][4] * sin(sita) * 202 EI += moon_longitude_parameters[i][4] * sin(sita) *
207 pow(E, fabs(moon_longitude_parameters[i][1])); 203 pow(E, fabs(moon_longitude_parameters[i][1]));
208 } 204 }
@@ -212,12 +208,12 @@ double parameter::moon_longitude(double t) {
212 a1 = 119.75 + 131.849 * T; 208 a1 = 119.75 + 131.849 * T;
213 a2 = 53.09 + 479264.290 * T; 209 a2 = 53.09 + 479264.290 * T;
214 210
215 EI += 3958.0 * sin(a1 * pi / 180); 211 EI += 3958.0 * sin(a1 * M_PI / 180);
216 EI += 1962.0 * sin((L - angles[3]) * pi / 180); 212 EI += 1962.0 * sin((L - angles[3]) * M_PI / 180);
217 EI += 318.0 * sin(a2 * pi / 180); 213 EI += 318.0 * sin(a2 * M_PI / 180);
218 214
219 L += EI / 1000000; 215 L += EI / 1000000;
220 L += 180 / pi * nutation(T); //地球章动修正 216 L += 180 / M_PI * nutation(T); //地球章动修正
221 217
222 while (L < 0) 218 while (L < 0)
223 L += 360; 219 L += 360;
generated by cgit v1.2.3-70-g09d2 (git 2.46.0) at 2025-04-08 14:12:09 +0800