usr.bin/calendar/chinese.c

   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 }