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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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