1 /*-
2 * SPDX-License-Identifier: BSD-3-Clause
3 *
4 * Copyright (c) 2019-2020 The DragonFly Project. All rights reserved.
5 *
6 * This code is derived from software contributed to The DragonFly Project
7 * by Aaron LI <aly@aaronly.me>
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions
11 * are met:
12 *
13 * 1. Redistributions of source code must retain the above copyright
14 * notice, this list of conditions and the following disclaimer.
15 * 2. Redistributions in binary form must reproduce the above copyright
16 * notice, this list of conditions and the following disclaimer in
17 * the documentation and/or other materials provided with the
18 * distribution.
19 * 3. Neither the name of The DragonFly Project nor the names of its
20 * contributors may be used to endorse or promote products derived
21 * from this software without specific, prior written permission.
22 *
23 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
26 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
27 * COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
28 * INCIDENTAL, SPECIAL, EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING,
29 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
31 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
32 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
33 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
34 * SUCH DAMAGE.
35 *
36 * Reference:
37 * Calendrical Calculations, The Ultimate Edition (4th Edition)
38 * by Edward M. Reingold and Nachum Dershowitz
39 * 2018, Cambridge University Press
40 */
41
42 /*
43 * Implement the Chinese calendar of the 1645 version, established in the
44 * second year of the Qīng dynasty.
45 *
46 * The winter solstice (dōngzhì; 冬至) always occurs during the eleventh
47 * month of the year. The winter-solstice-to-winter-solstice period is
48 * called a suì (岁).
49 *
50 * The leap month of a 13-month winter-solstice-to-winter-solstice period
51 * is the first month that does not contain a major solar term --- that is,
52 * the first lunar month that is wholly within a solar month.
53 */
54
55 #include <math.h>
56 #include <stdbool.h>
57 #include <stdio.h>
58
59 #include "calendar.h"
60 #include "basics.h"
61 #include "chinese.h"
62 #include "dates.h"
63 #include "gregorian.h"
64 #include "moon.h"
65 #include "sun.h"
66 #include "utils.h"
67
68 /*
69 * Fixed date of the start of the Chinese calendar.
70 * Ref: Sec.(19.3), Eq.(19.15)
71 */
72 static const int epoch = -963099; /* Gregorian: -2636, February, 15 */
73
74 /*
75 * Timezone (in fraction of days) of Beijing adopted in Chinese calendar
76 * calculations.
77 * Ref: Sec.(19.1), Eq.(19.2)
78 */
79 static double
chinese_zone(int rd)80 chinese_zone(int rd)
81 {
82 if (gregorian_year_from_fixed(rd) < 1929)
83 return (1397.0 / 180.0 / 24.0);
84 else
85 return (8.0 / 24.0);
86 }
87
88 /*
89 * The universal time (UT) of (clock) midnight at the start of the fixed date
90 * $rd in China.
91 * Ref: Sec.(19.1), Eq.(19.7)
92 */
93 static double
midnight_in_china(int rd)94 midnight_in_china(int rd)
95 {
96 return (double)rd - chinese_zone(rd);
97 }
98
99 /*
100 * Calculate the last Chinese major solar term (zhōngqì) in range of [1, 12]
101 * before the fixed date $rd.
102 * Ref: Sec.(19.1), Eq.(19.1)
103 */
104 static int
current_major_solar_term(int rd)105 current_major_solar_term(int rd)
106 {
107 double t_u = midnight_in_china(rd);
108 double lon = solar_longitude(t_u);
109 return mod1(2 + div_floor((int)lon, 30), 12);
110 }
111
112 /*
113 * Calculate the fixed date (in China) of winter solstice on or before the
114 * given fixed date $rd.
115 * Ref: Sec.(19.1), Eq.(19.8)
116 */
117 static int
chinese_winter_solstice_onbefore(int rd)118 chinese_winter_solstice_onbefore(int rd)
119 {
120 /* longitude of Sun at winter solstice */
121 const double winter = 270.0;
122
123 /* approximate the time of winter solstice */
124 double t = midnight_in_china(rd + 1);
125 double approx = estimate_prior_solar_longitude(winter, t);
126
127 /* search for the date of winter solstice */
128 int day = (int)floor(approx) - 1;
129 while (winter >= solar_longitude(midnight_in_china(day+1)))
130 day++;
131
132 return day;
133 }
134
135 /*
136 * Calculate the fixed date (in China) of the first new moon on or after the
137 * given fixed date $rd.
138 * Ref: Sec.(19.2), Eq.(19.9)
139 */
140 static int
chinese_new_moon_onafter(int rd)141 chinese_new_moon_onafter(int rd)
142 {
143 double t = new_moon_atafter(midnight_in_china(rd));
144 double t_s = t + chinese_zone(rd); /* in standard time */
145
146 return (int)floor(t_s);
147 }
148
149 /*
150 * Calculate the fixed date (in China) of the first new moon before the
151 * given fixed date $rd.
152 * Ref: Sec.(19.2), Eq.(19.10)
153 */
154 static int
chinese_new_moon_before(int rd)155 chinese_new_moon_before(int rd)
156 {
157 double t = new_moon_before(midnight_in_china(rd));
158 double t_s = t + chinese_zone(rd); /* in standard time */
159
160 return (int)floor(t_s);
161 }
162
163 /*
164 * Determine whether the Chinese lunar month starting on the given fixed
165 * date $rd has no major solar term.
166 * Ref: Sec.(19.2), Eq.(19.11)
167 */
168 static bool
chinese_no_major_solar_term(int rd)169 chinese_no_major_solar_term(int rd)
170 {
171 int rd2 = chinese_new_moon_onafter(rd + 1);
172 return (current_major_solar_term(rd) ==
173 current_major_solar_term(rd2));
174 }
175
176 /*
177 * Recursively determine whether there is a Chinese leap month on or
178 * after the lunar month starting on fixed date $m1 and at or before
179 * the lunar month starting on fixed date $m2.
180 * Ref: Sec.(19.2), Eq.(19.12)
181 */
182 static bool
chinese_prior_leap_month(int m1,int m2)183 chinese_prior_leap_month(int m1, int m2)
184 {
185 int m2_prev = chinese_new_moon_before(m2);
186 return ((m2 >= m1) &&
187 (chinese_no_major_solar_term(m2) ||
188 chinese_prior_leap_month(m1, m2_prev)));
189 }
190
191 /*
192 * Calculate the fixed date of Chinese New Year in suì containing the
193 * given date $rd.
194 * Ref: Sec.(19.2), Eq.(19.13)
195 */
196 static int
chinese_new_year_in_sui(int rd)197 chinese_new_year_in_sui(int rd)
198 {
199 /* prior and following winter solstice */
200 int s1 = chinese_winter_solstice_onbefore(rd);
201 int s2 = chinese_winter_solstice_onbefore(s1 + 370);
202
203 /* month after the 11th month: either 12 or leap 11 */
204 int m12 = chinese_new_moon_onafter(s1 + 1);
205 /* month after m12: either 12 (or leap 12) or 1 */
206 int m13 = chinese_new_moon_onafter(m12 + 1);
207 /* the next 11th month */
208 int m11_next = chinese_new_moon_before(s2 + 1);
209
210 bool leap_year = (lround((m11_next - m12) / mean_synodic_month) == 12);
211 if (leap_year && (chinese_no_major_solar_term(m12) ||
212 chinese_no_major_solar_term(m13))) {
213 /* either m12 or m13 is a leap month */
214 return chinese_new_moon_onafter(m13 + 1);
215 } else {
216 return m13;
217 }
218 }
219
220 /*
221 * Calculate the fixed date of Chinese New Year on or before the given
222 * fixed date $rd.
223 * Ref: Sec.(19.2), Eq.(19.14)
224 */
225 static int
chinese_new_year_onbefore(int rd)226 chinese_new_year_onbefore(int rd)
227 {
228 int newyear = chinese_new_year_in_sui(rd);
229 if (rd >= newyear)
230 return newyear;
231 else
232 return chinese_new_year_in_sui(rd - 180);
233 }
234
235 /*
236 * Calculate the fixed date of Chinese New Year in Gregorian year $year.
237 * Ref: Sec.(19.6), Eq.(19.26)
238 */
239 int
chinese_new_year(int year)240 chinese_new_year(int year)
241 {
242 struct date date = { year, 7, 1 };
243 int july1 = fixed_from_gregorian(&date);
244 return chinese_new_year_onbefore(july1);
245 }
246
247 /*
248 * Calculate the Chinese date (cycle, year, month, leap, day) corresponding
249 * to the fixed date $rd.
250 * Ref: Sec.(19.3), Eq.(19.16)
251 */
252 void
chinese_from_fixed(int rd,struct chinese_date * date)253 chinese_from_fixed(int rd, struct chinese_date *date)
254 {
255 /* prior and following winter solstice */
256 int s1 = chinese_winter_solstice_onbefore(rd);
257 int s2 = chinese_winter_solstice_onbefore(s1 + 370);
258
259 /* start of month containing the given date */
260 int m = chinese_new_moon_before(rd + 1);
261 /* start of the previous month */
262 int m_prev = chinese_new_moon_before(m);
263 /* month after the last 11th month: either 12 or leap 11 */
264 int m12 = chinese_new_moon_onafter(s1 + 1);
265 /* the next 11th month */
266 int m11_next = chinese_new_moon_before(s2 + 1);
267
268 bool leap_year = (lround((m11_next - m12) / mean_synodic_month) == 12);
269 int month = (int)lround((m - m12) / mean_synodic_month);
270 if (leap_year && chinese_prior_leap_month(m12, m))
271 month--;
272 month = mod1(month, 12);
273 bool leap_month = (leap_year &&
274 chinese_no_major_solar_term(m) &&
275 !chinese_prior_leap_month(m12, m_prev));
276 int elapsed_years = (int)floor(1.5 - month/12.0 +
277 (rd - epoch) / mean_tropical_year);
278
279 date->cycle = div_floor(elapsed_years - 1, 60) + 1;
280 date->year = mod1(elapsed_years, 60);
281 date->month = month;
282 date->leap = leap_month;
283 date->day = rd - m + 1;
284 }
285
286 /*
287 * Calculate the fixed date corresponding to the given Chinese date $date
288 * (cycle, year, month, leap, day).
289 * Ref: Sec.(19.3), Eq.(19.17)
290 */
291 int
fixed_from_chinese(const struct chinese_date * date)292 fixed_from_chinese(const struct chinese_date *date)
293 {
294 int midyear = (int)floor(epoch + mean_tropical_year *
295 ((date->cycle - 1) * 60 + date->year - 0.5));
296 int newyear = chinese_new_year_onbefore(midyear);
297
298 /* new moon before the given date */
299 int newmoon = chinese_new_moon_onafter(
300 newyear + (date->month - 1) * 29);
301 struct chinese_date date2;
302 chinese_from_fixed(newmoon, &date2);
303 if (date->month != date2.month || date->leap != date2.leap) {
304 /* there was a prior leap month, so get the next month */
305 newmoon = chinese_new_moon_onafter(newmoon + 1);
306 }
307
308 return newmoon + date->day - 1;
309 }
310
311 /*
312 * Calculate the fixed date of Qīngmíng (清明) in Gregorian year $g_year.
313 * Ref: Sec.(19.6), Eq.(19.28)
314 */
315 int
chinese_qingming(int g_year)316 chinese_qingming(int g_year)
317 {
318 double lambda = 15.0; /* Solar longitude of Qīngmíng */
319 struct date date = { g_year, 4, 1 }; /* Qīngmíng is around April 5 */
320 int rd = fixed_from_gregorian(&date);
321 double zone = chinese_zone(rd);
322 double t = solar_longitude_atafter(lambda, rd) + zone;
323 return (int)floor(t);
324 }
325
326 /**************************************************************************/
327
328 /* celestial stem, tiāngān (天干) */
329 static const struct stem {
330 const char *name;
331 const char *zhname;
332 } STEMS[] = {
333 { "Jiǎ", "甲" },
334 { "Yǐ", "乙" },
335 { "Bǐng", "丙" },
336 { "Dīng", "丁" },
337 { "Wù", "戊" },
338 { "Jǐ", "己" },
339 { "Gēng", "庚" },
340 { "Xīn", "辛" },
341 { "Rén", "壬" },
342 { "Guǐ", "癸" },
343 };
344
345 /* terrestrial branch, dìzhī (地支) */
346 static const struct branch {
347 const char *name;
348 const char *zhname;
349 const char *zodiac;
350 const char *zhzodiac;
351 } BRANCHES[] = {
352 { "Zǐ", "子", "Rat", "鼠" },
353 { "Chǒu", "丑", "Ox", "牛" },
354 { "Yín", "寅", "Tiger", "虎" },
355 { "Mǎo", "卯", "Rabbit", "兔" },
356 { "Chén", "辰", "Dragon", "龙" },
357 { "Sì", "巳", "Snake", "蛇" },
358 { "Wǔ", "午", "Horse", "马" },
359 { "Wèi", "未", "Goat", "羊" },
360 { "Shēn", "申", "Monkey", "猴" },
361 { "Yǒu", "酉", "Rooster", "鸡" },
362 { "Xū", "戌", "Dog", "狗" },
363 { "Hài", "亥", "Pig", "猪" },
364 };
365
366 /* Chinese names of months and days */
367 static const char *months[] = {
368 "正", "二", "三", "四", "五", "六",
369 "七", "八", "九", "十", "冬", "腊",
370 };
371 static const char *mdays[] = {
372 "初一", "初二", "初三", "初四", "初五", "初六",
373 "初七", "初八", "初九", "初十", "十一", "十二",
374 "十三", "十四", "十五", "十六", "十七", "十八",
375 "十九", "二十", "廿一", "廿二", "廿三", "廿四",
376 "廿五", "廿六", "廿七", "廿八", "廿九", "三十",
377 };
378
379 /* 24 major and minor solar terms (节气) */
380 static const struct chinese_jieqi jieqis[] = {
381 { "Lìchūn", "立春", false, 315 },
382 { "Yǔshuǐ", "雨水", true, 330 },
383 { "Jīngzhé", "惊蛰", false, 345 },
384 { "Chūnfēn", "春分", true, 0 },
385 { "Qīngmíng", "清明", false, 15 },
386 { "Gǔyǔ", "谷雨", true, 30 },
387 { "Lìxià", "立夏", false, 45 },
388 { "Xiǎomǎn", "小满", true, 60 },
389 { "Mángzhòng", "芒种", false, 75 },
390 { "Xiàzhì", "夏至", true, 90 },
391 { "Xiǎoshǔ", "小暑", false, 105 },
392 { "Dàshǔ", "大暑", true, 120 },
393 { "Lìqiū", "立秋", false, 135 },
394 { "Chǔshǔ", "处暑", true, 150 },
395 { "Báilù", "白露", false, 165 },
396 { "Qiūfēn", "秋分", true, 180 },
397 { "Hánlù", "寒露", false, 195 },
398 { "Shuāngjiàng", "霜降", true, 210 },
399 { "Lìdōng", "立冬", false, 225 },
400 { "Xiǎoxuě", "小雪", true, 240 },
401 { "Dàxuě", "大雪", false, 255 },
402 { "Dōngzhì", "冬至", true, 270 },
403 { "Xiǎohán", "小寒", false, 285 },
404 { "Dàhán", "大寒", true, 300 },
405 };
406
407 /*
408 * Calculate the fixed date (in China) of the following jiéqì on or after
409 * the given fixed date $rd, with the calculated jiéqì information stored
410 * in $jieqi.
411 * If $type is C_JIEQI_ALL, find either major or minor jiéqì;
412 * If $type is C_JIEQI_MAJOR, find only major jiéqì;
413 * If $type is C_JIEQI_MINOR, find only minor jiéqì.
414 */
415 int
chinese_jieqi_onafter(int rd,int type,const struct chinese_jieqi ** jieqi)416 chinese_jieqi_onafter(int rd, int type, const struct chinese_jieqi **jieqi)
417 {
418 const struct chinese_jieqi *jq1, *jq2;
419 const double zone = chinese_zone(rd);
420 double t_u = midnight_in_china(rd);
421 double lon = solar_longitude(t_u);
422 double lon1, lon2;
423 size_t n = nitems(jieqis);
424
425 for (size_t i = 0; i < n; i++) {
426 jq1 = &jieqis[i];
427 jq2 = &jieqis[(i+1) % n];
428 if ((type == C_JIEQI_MINOR && jq2->is_major) ||
429 (type == C_JIEQI_MAJOR && !jq2->is_major)) {
430 jq2 = &jieqis[(i+2) % n];
431 }
432
433 lon1 = jq1->longitude;
434 lon2 = jq2->longitude;
435 if (lon2 < lon1)
436 lon2 += 360;
437
438 if ((lon >= lon1 && lon < lon2) ||
439 (lon + 360 >= lon1 && lon + 360 < lon2)) {
440 break;
441 }
442 }
443
444 t_u = solar_longitude_atafter(jq2->longitude, floor(t_u));
445 *jieqi = jq2;
446 return (int)floor(t_u + zone);
447 }
448
449
450 /*
451 * Format the Chinese date of the given fixed date $rd in $buf.
452 * Return the formatted string length.
453 */
454 int
chinese_format_date(char * buf,size_t size,int rd)455 chinese_format_date(char *buf, size_t size, int rd)
456 {
457 struct chinese_date cdate;
458 chinese_from_fixed(rd, &cdate);
459 return snprintf(buf, size, "%s%s月%s",
460 cdate.leap ? "闰" : "", months[cdate.month - 1],
461 mdays[cdate.day - 1]);
462 }
463
464
465 /*
466 * Find days of the specified Chinese month ($month) and day ($day).
467 * If month $month < 0, then month is ignored (i.e., all months).
468 * (NOTE: The year $year is ignored.)
469 */
470 int
chinese_find_days_ymd(int year __unused,int month,int day,struct cal_day ** dayp,char ** edp __unused)471 chinese_find_days_ymd(int year __unused, int month, int day,
472 struct cal_day **dayp, char **edp __unused)
473 {
474 struct cal_day *dp;
475 struct chinese_date cdate;
476 int rd, rd_begin, rd_end;
477 int count = 0;
478
479 rd_begin = chinese_new_moon_before(Options.day_begin);
480 rd_end = chinese_new_moon_onafter(Options.day_end);
481 rd = rd_begin;
482 while (rd <= rd_end) {
483 chinese_from_fixed(rd, &cdate);
484 if (cdate.month == month || month < 0) {
485 cdate.day = day;
486 rd = fixed_from_chinese(&cdate);
487 if ((dp = find_rd(rd, 0)) != NULL) {
488 if (count >= CAL_MAX_REPEAT) {
489 warnx("%s: too many repeats",
490 __func__);
491 return count;
492 }
493 dayp[count++] = dp;
494 }
495 }
496
497 /* go to next month */
498 rd = chinese_new_moon_onafter(rd + (day == 0 ? 2 : 1));
499 }
500
501 return count;
502 }
503
504 /*
505 * Find days of the specified Chinese day of month ($dom) of all months.
506 */
507 int
chinese_find_days_dom(int dom,struct cal_day ** dayp,char ** edp __unused)508 chinese_find_days_dom(int dom, struct cal_day **dayp, char **edp __unused)
509 {
510 return chinese_find_days_ymd(-1, -1, dom, dayp, NULL);
511 }
512
513
514 /*
515 * Print the Chinese calendar of the given date $rd and events of the year.
516 */
517 void
show_chinese_calendar(int rd)518 show_chinese_calendar(int rd)
519 {
520 struct chinese_date date;
521 chinese_from_fixed(rd, &date);
522 /* Ref: Sec.(19.3), Eq.(19.18) */
523 struct stem stem = STEMS[mod1(date.year, 10) - 1];
524 struct branch branch = BRANCHES[mod1(date.year, 12) - 1];
525
526 struct date gdate;
527 gregorian_from_fixed(rd, &gdate);
528 int g_year = gdate.year;
529 printf("公历 (Gregorian): %d-%02d-%02d\n",
530 gdate.year, gdate.month, gdate.day);
531
532 printf("农历: %s%s年 [%s], %s%s月%s\n",
533 stem.zhname, branch.zhname, branch.zhzodiac,
534 date.leap ? "闰" : "", months[date.month - 1],
535 mdays[date.day - 1]);
536 printf("Chinese calendar: year %s%s [%s], %smonth %d, day %d\n",
537 stem.name, branch.name, branch.zodiac,
538 date.leap ? "leap " : "", date.month, date.day);
539
540 /* the following Chinese New Year */
541 int newyear = chinese_new_year_onbefore(
542 chinese_new_year_onbefore(rd) + 30*13);
543 gregorian_from_fixed(newyear, &gdate);
544 printf("新年 (Chinese New Year): %d-%02d-%02d\n",
545 gdate.year, gdate.month, gdate.day);
546
547 const double zone = chinese_zone(rd);
548 const struct chinese_jieqi *jq;
549 double t_jq;
550 char buf_time[32], buf_zone[32];
551 int lambda;
552
553 /* 1st solar term (Lìchūn) is generally around February 4 */
554 date_set(&gdate, g_year, 2, 1);
555 t_jq = (double)fixed_from_gregorian(&gdate);
556
557 format_zone(buf_zone, sizeof(buf_zone), zone);
558
559 printf("\n二十四节气 (solar terms):\n");
560 for (size_t i = 0; i < nitems(jieqis); i++) {
561 jq = &jieqis[i];
562 lambda = jq->longitude;
563 t_jq = solar_longitude_atafter(lambda, floor(t_jq)) + zone;
564 gregorian_from_fixed(floor(t_jq), &gdate);
565 format_time(buf_time, sizeof(buf_time), t_jq);
566
567 printf("%s (%-13s): %3d°, %d-%02d-%02d %s %s\n",
568 jq->zhname, jq->name, lambda,
569 gdate.year, gdate.month, gdate.day,
570 buf_time, buf_zone);
571 }
572 }
573