第七章 儒略日 Julian Day

1. 儒略历和格里高利历

1.1 定义

儒略历，是格里历的前身，由罗马共和国独裁官儒略·凯撒采纳埃及亚历山大的希腊数学家兼天文学家索西琴尼计算的历法，在公元前45年1月1日起执行，取代旧罗马历历法的历法。一年设12个月，大小月交替，四年一闰，平年365日，闰年于二月底增加一闰日，年平均长度为365.25日。由于累积误差随着时间越来越大，1582年后由教皇格里高利十三世改良，变为格里历，即沿用至今的公历。但大英帝国、北美十三州等直到1752年才从儒略历改用格里历。现今儒略历只有苏格兰昔德兰群岛之富拉岛、阿索斯神权共和国和一些北非的柏柏尔人使用。

格里历（拉丁语：Calendarium Gregorianum，又译、国瑞历、额我略历、格列高利历、格里高利历、葛瑞格里历、格列高历，也称基督历[1]），是由意大利医生兼哲学家阿洛伊修斯·里利乌斯改革儒略历制定的历法，由罗马大公教会教宗格列高利十三世在1582年颁行。公历是阳历的一种，于1912年在中国引进采用，因农历等中国传统历法是阴阳历，故公历在中文中又称阳历、西历、新历、公历。格里历与儒略历一样，格里历也是每四年在2月底置一闰日，但格里历特别规定，除非能被400整除，所有的世纪年（能被100整除）都不设闰日；如此，每四百年，格里历仅有97个闰年，比儒略历减少3个闰年 [注 1]。 格里历的历年平均长度为365.2425日，接近平均回归年的365.242199074日，即约每3300年误差一日，也更接近春分点回归年的365.24237日，即约每8000年误差一日；而儒略历的历年为365.25日，约每128年就误差一日[注 1]。到1582年时，儒略历的春分日（3月21日）与地球公转到春分点的实际时间已相差10天。因此，格里历开始实行时，将儒略历1582年10月4日星期四的次日，为格里历1582年10月15日星期五，即有10天被删除，但原星期的周期保持不变。格里历的纪年沿用儒略历，自传统的耶稣诞生年开始，称为“公元”，亦称“西元”。

1.2 区别

　　由定义可见，儒略历和格里高利历都属于阳历，后者是以太阳回归年为基础，对前者置润误差进行调整后的历法。所以两者之间最大的不同就是置润规则。

1.2.1 置润规则

• 儒略历

  能被4整除的年份为闰年（产生多润的原因，从而使儒略历在日期上落后于格里历）

• 格里历

  能被4整除但不能被100整除的非世纪年 + 能被400整除的世纪年，即要在儒略历闰年中扣除那些不能被400整除的世纪年（如1600，1700等）

// LeapYearJulian returns true if year y in the Julian calendar is a leap year.
// 儒略历闰年判断
func LeapYearJulian(y int) bool {
	return y%4 == 0
}

// LeapYearGregorian returns true if year y in the Gregorian calendar is a leap year.
// 格里历闰年判断
func LeapYearGregorian(y int) bool {
	return (y%4 == 0 && y%100 != 0) || y%400 == 0
}

1.2.2 转换

　　知道了这个不同，我们就能推导出两种历法之间实际的相差天数，如下：

• 1582年：
  格里历10月15日，合儒略历10月5日，或之后的日期：格里历日期减10日等于儒略历日期。
• 1583年——1699年：
  格里历日期减10日等于儒略历日期。
• 1700年（格里历没有闰日，但儒略历有）：
  格里历2月28日，合儒略历2月18日，或之前的日期：格里历日期减10日等于儒略历日期。
  格里历3月1日，合儒略历2月19日，或之后的日期：格里历日期减11日等于儒略历日期。
• 1701年——1799年：
  格里历日期减11日等于儒略历日期。
• 1800年（格里历没有闰日，但儒略历有）：
  格里历2月28日，合儒略历2月17日，或之前的日期：格里历日期减11日等于儒略历日期。
  格里历3月1日，合儒略历2月18日，或之后的日期：格里历日期减12日等于儒略历日期。
• 1801年——1899年：
  格里历日期减12日等于儒略历日期。
• 1900年（格里历没有闰日，但儒略历有）：
  格里历2月28日，合儒略历2月16日，或之前的日期：格里历日期减12日等于儒略历日期。
  格里历3月1日，合儒略历2月17日，或之后的日期：格里历日期减13日等于儒略历日期。
• 1901年——2099年：
  格里历日期减13日等于儒略历日期。
• 2100年（格里历没有闰日，但儒略历有）：
  格里历2月28日，合儒略历2月15日，或之前的日期：格里历日期减13日等于儒略历日期。
  格里历3月1日，合儒略历2月16日，或之后的日期：格里历日期减14日等于儒略历日期。

2. 儒略日 Julian Day

• 儒略日（Julian Day）
  儒略日是在儒略周期内以连续的日数计算时间的计时法，主要是天文学家在使用。

• 儒略日数（Julian Day Number，JDN）
  儒略日数的计算是从格林威治标准时间的中午开始，包含一个整天的时间，起点的时间（0日）回溯至儒略历的公元前4713年1月1日中午12点（在格里历是公元前4714年11月24日），这个日期是三种多年周期的共同起点，且是历史上最接近现代的一个起点。例如，2000年1月1日的UT12:00是儒略日2,451,545。

• 儒略日期（Julian date，JD）
  儒略日期是以格林威治标准时中午12:00的儒略日加上那一天的瞬时时间的分数。儒略日期是儒略日添加小数部分所表示的儒略日数。例如，2013年1月1日00:30:00（UT）是儒略日期2,456,293.520833。

• 儒略周期（Julian Period）
  儒略周期是开始于公元前4713年，长达7980年的纪年法，被用于历史上各种不同历法的日期转换。公元2018年是儒略周期的6731年，下一个儒略周期将开始于公元3268年。

• 儒略日起点
  儒略日的起点订在西元前4713年（天文学上记为 -4712年）1月1日格林威治时间平午（世界时12:00），即JD 0指定为UT时间B.C.4713年1月1日12:00到UC时间B.C.4713年1月2日12:00的24小时。每一天赋予了一个唯一的数字，顺数而下，如：1996年1月1日12:00:00的儒略日是2450084。这个日期是考虑了太阳、月亮的轨道运行周期，以及当时收税的间隔而订出来的。Joseph Scaliger定义儒略周期为7980年，是因28、19、15的最小公倍数为28×19×15=7980。其中：

  28年为一太阳周期（solar cycle），经过一太阳周期，则星期的日序与月的日序会重复。

  19年为一太阴周期，或称默冬章（Metonic cycle），因235朔望月=19回归年，经过一太阴周期则阴历月年的日序重复。

  15年为一小纪（indiction cycle），此为罗马皇帝君士坦丁一世所颁，每15年评定财产价值以供课税，成为古罗马用的一个纪元单位，

  故以7980年为一儒略周期，而所选的起点公元前4713年，则是这三个循环周期同时开始的最近年份。

• 简化儒略日（MJD）
  由于儒略日数字位数太多，国际天文学联合会于1973年采用简化儒略日（MJD），其定义为MJD = JD - 2400000.5。MJD相应的起点是1858年11月17日世界时0时。

　　算法中主要用的是 儒略日期（Julian date，JD），后续简称为儒略日

3. 儒略日和阳历日期（格里历、儒略历日期）的转换

3.1 阳历日期$(Y,M,D)$ $\Rightarrow$ 儒略日$(JD)$

1. 设$Y$为给定年份，$M$为月份，$D$为该月日期（可以带小数）。

2. 若 $M > 2$，$Y$和$M$不变，

3. 若 $M =1或2$，以$Y–1$代$Y$，以$M+12$代$M$，换句话说，如果日期在1月或2月，则被看作是在前一年的13月或14月。

4. 对格里历有 ：$A = \lfloor \frac {Y}{100}\rfloor$ $B = 2 - A + \lfloor \frac{A}{4}\rfloor$

5. 对儒略历，取 $B = 0$

6. 要求的儒略日即为：$$JD = \lfloor 365.25(Y+4716)\rfloor+\lfloor30.6001(M+1)\rfloor+D+B-1524.5$$

此处作对30.6001取值的解释是为了确保小数在运算时的准确性，如30.6 * 5 应该等于153，但是很多计算机得出的结果是152.999 9998。

另一种方式就是用306代替，最后再除以10取整

// FloorDiv returns the integer floor of the fractional value (x / y).
//
// It uses integer math only, so is more efficient than using floating point
// intermediate values.  This function can be used in many places where INT()
// appears in AA.  As with built in integer division, it panics with y == 0.
func FloorDiv(x, y int) (q int) {
	q = x / y
	if (x < 0) != (y < 0) && x%y != 0 {
		q--
	}
	return
}

// FloorDiv64 returns the integer floor of the fractional value (x / y).
//
// It uses integer math only, so is more efficient than using floating point
// intermediate values.  This function can be used in many places where INT()
// appears in AA.  As with built in integer division, it panics with y == 0.
func FloorDiv64(x, y int64) (q int64) {
	q = x / y
	if (x < 0) != (y < 0) && x%y != 0 {
		q--
	}
	return
}

// CalendarGregorianToJD converts a Gregorian year, month, and day of month
// to Julian day.
//
// Negative years are valid, back to JD 0.  The result is not valid for
// dates before JD 0.
// 格里历日期转儒略日
func CalendarGregorianToJD(y, m int, d float64) float64 {
	switch m {
	case 1, 2:
		y--
		m += 12
	}
	a := base.FloorDiv(y, 100)
	b := 2 - a + base.FloorDiv(a, 4)
	// (7.1) p. 61
	return float64(base.FloorDiv64(36525*(int64(y+4716)), 100)) +
		float64(base.FloorDiv(306*(m+1), 10)+b) + d - 1524.5
}

// CalendarJulianToJD converts a Julian year, month, and day of month to Julian day.
//
// Negative years are valid, back to JD 0.  The result is not valid for
// dates before JD 0.
// 儒略历日期转儒略日
func CalendarJulianToJD(y, m int, d float64) float64 {
	switch m {
	case 1, 2:
		y--
		m += 12
	}
	return float64(base.FloorDiv64(36525*(int64(y+4716)), 100)) +
		float64(base.FloorDiv(306*(m+1), 10)) + d - 1524.5
}

3.2 儒略日$(JD)$ $\Rightarrow$ 阳历日期$(y,m,d)$

1. 将 $JD$ 加上 0.5，令$Z$为其整数部分，$F$ 为尾数（小数）部分。
2. 若 $Z$ < 2299161，取 $A$ = $Z$;
   若 $Z$ 大于等于 2299161，计算 \begin{align} α &=\lfloor \frac {Z-1867216.25}{36524.25}\rfloor\\[2ex] A &=Z+1+α-\lfloor \frac α4\rfloor \end{align}
3. 然后计算 \begin{align} B &= A+1524\\[2ex] C &= \lfloor \frac {B-122.1}{365.25}\rfloor\\[2ex] D &= \lfloor 365.25C\rfloor\\[2ex] E &= \lfloor \frac {B-D}{30.6001}\rfloor\\[2ex] \end{align}
4. 该月日期（带小数部分）则为：$$d = B - D - \lfloor 30.6001E\rfloor + F$$
5. 月份 m 为： \begin{cases} m = E – 1 ;& \text {if $E$ < 14}\\[2ex] m = E – 13 ;& \text {if $E$ $\geq$ 14} \end{cases}
6. 年份为 y： \begin{cases} y = C – 4716 ;& \text {if $m$ > 2}\\[2ex] y = C – 4715 ;& \text {if $m$ $\leq$ 2} \end{cases}

// JDToCalendar returns the calendar date for the given jd.
//
// Note that this function returns a date in either the Julian or Gregorian
// Calendar, as appropriate.
// 儒略日转公历日期
// 如果儒略日对应的格里历时间点在1582-10-15 12点之前，则转为儒略历日期
// 如果儒略日对应的格里历时间点在1582-10-15 12点之后，则转为格里历日期
func JDToCalendar(jd float64) (year, month int, day float64) {
	zf, f := math.Modf(jd + .5)
	z := int64(zf)
	a := z
	// if z >= 2299151 { // typo 应该是2299161对应于现实格里历起始日1582-10-15中午12点
	if z >= 2299161 { // typo 应该是2299161对应于现实格里历起始日1582-10-15中午12点
		α := base.FloorDiv64(z*100-186721625, 3652425)
		a = z + 1 + α - base.FloorDiv64(α, 4)
	}
	b := a + 1524
	c := base.FloorDiv64(b*100-12210, 36525)
	d := base.FloorDiv64(36525*c, 100)
	e := int(base.FloorDiv64((b-d)*1e4, 306001))
	// compute return values
	day = float64(int(b-d)-base.FloorDiv(306001*e, 1e4)) + f
	switch e {
	default:
		month = e - 1
	case 14, 15:
		month = e - 13
	}
	switch month {
	default:
		year = int(c) - 4716
	case 1, 2:
		year = int(c) - 4715
	}
	return
}

// jdToCalendarGregorian returns the Gregorian calendar date for the given jd.
//
// Note that it returns a Gregorian date even for dates before the start of
// the Gregorian calendar.  The function is useful when working with Go
// time.Time values because they are always based on the Gregorian calendar.
// 始终转为格里历，忽略儒略历转格里历被扣除的那10天，即把1582-10-15 12点之前的日期也当成格里历算
func jdToCalendarGregorian(jd float64) (year, month int, day float64) {
	zf, f := math.Modf(jd + .5)
	z := int64(zf)
	α := base.FloorDiv64(z*100-186721625, 3652425)
	a := z + 1 + α - base.FloorDiv64(α, 4)
	b := a + 1524
	c := base.FloorDiv64(b*100-12210, 36525)
	d := base.FloorDiv64(36525*c, 100)
	e := int(base.FloorDiv64((b-d)*1e4, 306001))
	// compute return values
	day = float64(int(b-d)-base.FloorDiv(306001*e, 1e4)) + f
	switch e {
	default:
		month = e - 1
	case 14, 15:
		month = e - 13
	}
	switch month {
	default:
		year = int(c) - 4716
	case 1, 2:
		year = int(c) - 4715
	}
	return
}

// JDToTime takes a JD and returns a Go time.Time value.
// 儒略历转为格里历日期对应的 Go Time 类型
func JDToTime(jd float64) time.Time {
	// time.Time is always Gregorian
	y, m, d := jdToCalendarGregorian(jd)
	t := time.Date(y, time.Month(m), 0, 0, 0, 0, 0, time.UTC)
	return t.Add(time.Duration(d * 24 * float64(time.Hour)))
}

// TimeToJD takes a Go time.Time and returns a JD as float64.
//
// Any time zone offset in the time.Time is ignored and the time is
// treated as UTC.
// Go Time 类型转为儒略日（默认日期为格里历）
func TimeToJD(t time.Time) float64 {
	ut := t.UTC()
	y, m, _ := ut.Date()
	d := ut.Sub(time.Date(y, m, 0, 0, 0, 0, 0, time.UTC))
	// time.Time is always Gregorian
	return CalendarGregorianToJD(y, int(m), float64(d)/float64(24*time.Hour))
}

3.3 儒略日$(JD)$ $\Rightarrow$ 星期几

　　因为儒略日0.0(-4712-1-1 12:00) 是星期一，往后7天一循环就能得出星期几了
　　计算该日0时的儒略日，加上1.5，再除以7，所得余数将指示出星期几：
若余数为0，则为星期日，1为星期一，2为星期二，3为星期三，4为星期四，5为星期五，6为星期六。
　　儒略历到格里高利历的换算并不影响星期。 因而，在1582年10月4日星期四接下来的一天便是10月15日星期五。

// DayOfWeek determines the day of the week for a given JD.
//
// The value returned is an integer in the range 0 to 6, where 0 represents
// Sunday.  This is the same convention followed in the time package of the
// Go standard library.
// 儒略日0.0为-4712-1-1 12:00 星期一，往后7天一循环就能得出星期几了
func DayOfWeek(jd float64) int {
	return int(jd+1.5) % 7
}

3.4 日阳历日期$(Y,M,D)$ $\Rightarrow$ 年内序数日$(N)$

　　年内的序数日$N$可由以下公式得出：$$N = \lfloor \frac {275M}{9}\rfloor - K\lfloor \frac {M+9}{12}\rfloor + D - 30 $$ 　　此处$M$为月份，$D$为该月日期，闰年$K = 1$，平年$K = 2$
　　$N$ 取整数，自1月1日开始取值1，直至12月31日取值365（或闰年取值366）。

　　如何推导出这个公式？
　　乍一看挺没头绪的，我们先来看一个数列$$ a_m = \lfloor \frac {5(m+1)}{9}\rfloor, (m=0,1,2,\dots,12)$$

\begin{array}{c|ccl} m & a_m & a_m - a_{m-1} & 当月天数 \\
\hline 0 & 0 & - & - \\
1 & 1 & 1 & 31 \\
2 & 1 & 0 & 30(看做30) \\
3 & 2 & 1 & 31 \\
4 & 2 & 0 & 30 \\
5 & 3 & 1 & 31 \\
6 & 3 & 0 & 30 \\
7 & 4 & 1 & 31 \\
8 & 5 & 1 & 31 \\
9 & 5 & 0 & 30 \\
10 & 6 & 1 & 31 \\
11 & 6 & 0 & 30 \\
12 & 7 & 1 & 31 \end{array}

　　如果把$m$看成是月份，并且暂时把二月看成是30天，那么第二列恰好是$当月天数-30$，所以我们就可以用以下公式来表示每个月的天数。 $$D_m = 30 + a_m - a_{m-1}$$

　　所以从1月到m月的总天数： \begin{align} S_m = & D_1+D_2+,\dots,+D_m\\[2ex] = & 30 + a_m - a_{m-1} + \\[2ex] & 30 + a_{m-1} - a_{m-2} + \\[2ex] & \cdots \\[2ex] & 30 + a_2 - a_1 + \\[2ex] & 30 + a_1 - a_0 \\[2ex] = & 30m + a_m - a_0\\[2ex] = & 30m + \lfloor \frac {5(m+1)}{9}\rfloor \end{align}

　　由于上述表达式只有在2月等于30天的情况下才成立，所以当$m\geq 2$时还要根据平年、闰年分别减去相应的天数2，1，得到1月到$m$月的总天数和为： \begin{align} S_m=30m + \lfloor \frac {5(m+1)}{9}\rfloor - \lfloor \frac {m+10}{12}\rfloor k,\ 其中 \end{align} \begin{cases} k = 1 , \ 闰年\\[2ex] k = 2 , \ 平年 \end{cases}

　　所以$M$月$D$日的年内序数$N$为： \begin{align} N &= S_{M-1} + D\\[2ex] &= 30(M-1) + \lfloor \frac {5M}{9}\rfloor - \lfloor \frac {M+9}{12}\rfloor k + D\\[2ex] &= \lfloor \frac {275M}{9}\rfloor - \lfloor \frac {M+9}{12}\rfloor k -30 + D\\[2ex] \end{align}

// DayOfYearGregorian computes the day number within the year of the Gregorian
// calendar.
// 格里历年内序数
func DayOfYearGregorian(y, m, d int) int {
	return DayOfYear(y, m, d, LeapYearGregorian(y))
}

// DayOfYearJulian computes the day number within the year of the Julian
// calendar.
// 儒略历年内序数
func DayOfYearJulian(y, m, d int) int {
	return DayOfYear(y, m, d, LeapYearJulian(y))
}

// DayOfYear computes the day number within the year.
//
// This form of the function is not specific to the Julian or Gregorian
// calendar, but you must tell it whether the year is a leap year.
// 输入闰年标识，计算年内序数
func DayOfYear(y, m, d int, leap bool) int {
	k := 2
	if leap {
		k--
	}
	return wholeMonths(m, k) + d
}

// m月之前的所有月份天数之和
func wholeMonths(m, k int) int {
	return 275*m/9 - k*((m+9)/12) - 30
}

3.5 年内序数日$(N)$ $\Rightarrow$ 阳历日期$(M,D)$

// DayOfYearToCalendar returns the calendar month and day for a given
// day of year and leap year status.
// 年内序数求对应的日期
func DayOfYearToCalendar(n int, leap bool) (m, d int) {
	k := 2
	if leap {
		k--
	}
	if n < 32 {
		m = 1
	} else {
		m = (900*(k+n) + 98*275) / 27500
	}
	return m, n - wholeMonths(m, k)
}