月齢カレンダー

私がこの数年愛用している予定表の話はこのブログに何回か書いた. ところで, 今年の予定表の表紙にある週日定数も月齢定数も共に0である. つまり下のようなのだ.

そうするとこの次にこうなるのは何年後かも気になる.

y年の週日定数, week(y)と月齢定数, moon(y)は次の式で計算する.

def week(y): return (y-2000+(y-2000)//4 
    +4)%7
def moon(y): return (((y-11)%19)*11)%29

そこで次に(0,0)になる年を計算してみたら, なんと2272年; 247サロス周期の後であった. 今年は昭和100年, 戦後80年, 阪神淡路大震災30年といわれるが, 暦屋にとっても誠に貴重な年なのである.

この表を見た序でに, 今年の復活祭を調べよう. まず3月21日の月齢は21. 次の満月は月齢45として4月14日. その日の曜日は1で月曜だから, 日曜はその6日後で4月20日になる.

もうひとつ序でに. Winning WaysのDoomsdayは3月0日の曜日なので, 今年は5, つまり金曜日である.

月齢カレンダー

このブログの2008年12月に月齢計算のことを書いた. それは計算機を使って月齢を計算する自己流アルゴリズムであったが, その昔, 私は高校生の頃, 先輩に教えられて毎年月齢定数というのを使って月齢を概算していた.

それではその年の月齢定数をCとすると, その年のm月d日の近似的な月齢を (C + m + d) mod 30 で得るものであった.

それで時々年の初めころ, 今年の月齢定数はいくつかなとほどほどに計算して記憶していたこともあった.

ところでインターネットであちこち見ていたら堀さんの式というのに出会った. http://www.asj.or.jp/geppou/archive_open/1968/pdf/19680704.pdf

詳しくは元記事を見ていただくことにして, 大体は月齢定数と同様の趣旨であるが, mの代わりに f(m)を足すのである. この値が1月から順に0,2,0,2,2,4,5,6,7,8,9,10なのである. 6月から後はf(m) = m - 2なので, その辺は昔の月齢定数と同様だが, 年初がすこし違っていて, この式の方がより正確になるらしい. 元記事のタイトル「鬼鬼西」は1月から6月までの値の覚え方である.

Cの方は, (((y - 11) mod 19) × 11) mod 30 だそうで, 本年2016年なら 2005 mod 19 = 10なので, 110 mod 30 = 20が今年の月齢定数になる. 元記事の掲載された1968年はこの定数がちょうど0だったのでタイムリーな発表であった.

f(2)=2だから (20+2+8) mod 30 = 0 となり, 2月8日が新月つまり春節であった. 3月の満月をxとすると, 20 + 0 + x = 15 (mod 30) だからx = 25となって, 今年の復活祭は3月27日である. また3月9日には部分日蝕があるが, その月齢は 20 + 0 + 9 = 29 だ.

なおこの定数は 2005 mod 19 = 10 から後8年はこの値が1ずつ増え, 従って定数も毎年11ずつ増えるから初めから計算しなおすことはなく, 来年2017年は1, 2018年は12, 2019年は23,...と分かる. こうして..., 7, 18と来たら, その次は29ではなく0とする(なる?)のだ. 次回0なのは2025年である. 1968年, 1987年, 2006年, 2025年と19年ごとに定数が0の年が回って来る.

月齢カレンダー

直前のブログ「月齢カレンダー」の月齢計算で, 「ある定数」の決め方がなんとなく杜撰に思われるので, もう少し説得力ある計算法を考えた.

ある日の正午の月齢は, 直前の新月からの時間である. NAOJのウェブページを見ると, 2008年12月27日21時22分が新月だそうだ. すると12月28日正午までの経過時間は14時間38分. 24時間に対しては0.6097...となり, つまりそれが28日の月齢である.

従って望む日の月齢は, 1朔望月を29日499/940とすると, 望む日のRDから2008年12月28日のRDを引き, その29日499/940の剰余プラス0.6097となる.

(define (rd y m d) ;rdの計算
  (let* ((p (lambda (n) (= (modulo y n) 0)))
    (q (lambda (n) (quotient y n)))
    (s `(0 31 ,(if (if (p 100) (p 400) (p 4)) 29 28)
      31 30 31 30 31 31 30 31 30 31)))
    (- (+ d (* y 365) (q 4) (q -100) (q 400)) 
      (apply + (list-tail s m)))))

(define (ma y m d) ;月齢の計算
  (+ (/ (modulo (* (- (rd y m d) (rd 2008 12 28)) 940)
    27759) 940.) 0.6097))

(do ((i 1 (+ i 1))) ((> i 12)) ;1日の月齢表の計算
  (display (/ (round (* (ma 2009 i 1) 10)) 10)) (newline))

4.6
6.1
4.5
6.
6.5
8.
8.4
9.9
11.4
11.8
13.3
13.8

一方19年7閏法と有理化した朔望月の求め方は次の通り.

1太陽年=365.242198日
1朔望月=29.530589日
1太陽年から12朔望月を引くと10.875166日

これが1朔望月の何分の何に相当するか知りたいので, この比を連分数展開する.

小数xの連分数をn項とるのに

(define (cf x n)
 (if (> n 0)
  (let* ((y (/ 1 x)) (z (inexact->exact (truncate y))))
   (cons z (cf (- y z) (- n 1)))) '()))
として
(cf (/ 10.875144 29.530589) 10)
を計算すると
(2 1 2 1 1 17 2 2 2 2)
が得られる. 17分の1辺りで棄てることにし,
(define (cf2r l)
 (if (null? l) 0
  (/ 1 (+ (car l) (cf2r (cdr l))))))
で有理数を計算すると
(cf2r '(2 1 2 1 1))→ 7/19

となり, 余りは1朔望月の7/19であることが分かる. 従って19年に7回閏月を置くのである.

次に19太陽年は 365 1/4日の19倍となる. この時間に朔望月は19×12+7回(つまり235回)あるので,

(/ (* (+ 365 (/ 1 4)) 19) (+ (* 12 19) 7))→ 27759/940
(quotient 27759 940)→29
(modulo 27759 940)→499

すなわち29日499/940となる.

追記 上のような月齢の式を書いたが, 最後に0.6097を足すのはこれまたウンシャンである. その回避法は2つあり, まずは2008年12月28日のRDではなく, 2008年12月27日21時22分のRDを基準にすること. つぎは正午ころに新月のある日を見つけることである. たまたま手元の天文年鑑を見ていたら, 2005年8月5日の新月は12時5分であった. これはしたり! 1日は86400秒であり, 5分は300秒なので, 誤差は無視出来よう. 従って

(define (ma y m d) (/ (modulo                                        
(* (- (rd y m d) (rd 2005 8 5)) 940) 27759) 940.))  

はそれほどウンシャンではない式だと思っている.

月齢カレンダー

これまでカレンダーのプログラムは数えきれないほど書いている. 2007年の夏のプロシンのお題もカレンダーであった. Unixのコマンドcalのように, 年と月を入力し, 例えば

   December 2008
 S  M Tu  W Th  F  S
    1  2  3  4  5  6
 7  8  9 10 11 12 13
14 15 16 17 18 19 20
21 22 23 24 25 26 27
28 29 30 31

を出力するものだ. プロシンのときは, mss, scheme, haskell, teco, tex, postscriptで書いた.

最近は不況のせいか, いや私にカレンダーを送ってくれる同僚が殆んど退職したのであろうが, 富士通の「世界の車窓」以外は, カレンダーもあまり貰わなくなった. それでこの数日カレンダーを自作しようかと考えていた.

以前はpcalというプログラムがあったが, まぁそういうものは, 自分で書く方が早い. ところで, そのpcalのカレンダーに, 月齢の絵がついているのを見つけた. lcalともいうらしい. それはどうすれば書けるか, 思いがけなく熱中した.

月齢カレンダーを書くには, 2つの問題を解決しなければならない. まずy年m月d日の月齢の計算法. 次に月相の表示法である.

とりあえずサーチすると月齢の計算法は見つかるが, どの式を見てもウンシャンである. いろいろ考えた末, 私が採用したのは次のものである.

カレンダーに図示するのだから, それほど正確である必要はない. それよりエレガントに計算したい. まずy年m月d日のある起点からの日数がいる. それにはReingold, DershowitzのCalendrical CalculationsにあるRD(Rata Die (fixed date))を使うことにする. これはGregorian暦が昔まで使われたとして, 西暦1年1月1日(月曜)からの通日である.

一方, 朔望月は29.53...日であるが, これはTAOCPの復活祭公式にあるように, 有理数29 449/940にする. つまりRDにある定数を足し, 有理朔望月で除した剰余を月齢とするのである.

ある定数が問題であるが, たとえば自宅にあった2004年と2008年の理科年表の太陽, 月の表の正午の月齢からキャリブレートすれば得られる. こうして計算したのが, 以下の表で, 各年は1月1日から毎月1日の月齢で, 左が私の計算, 右が2004年と2008年は理科年表, 2009年と2010年はどこかのウェブページによるものである. まぁよかろう.

   2004        2008        2009       2010
 8.5  8.7 | 22.5 22.4 |  4.6  4.6 | 15.2 15.6
10.0 10.2 | 24.0 23.6 |  6.1  5.8 | 16.7 16.8
 9.5  9.7 | 23.5 23.0 |  4.6  4.1 | 15.2 15.0
10.9 11.2 | 24.9 24.4 |  6.0  5.5 | 16.7 16.3
11.4 11.6 | 25.4 25.0 |  6.5  6.0 | 17.1 16.6
12.9 12.9 | 26.9 26.6 |  8.0  7.6 | 18.6 18.1
13.3 13.3 | 27.3 27.3 |  8.4  8.3 | 19.1 18.7
14.8 14.7 | 28.8 29.0 |  9.9 10.0 | 20.5 20.3
16.3 16.1 |  0.7  1.3 | 11.4 11.7 | 22.0 22.0
16.8 16.5 |  1.2  1.8 | 11.8 12.3 | 22.5 22.7
18.2 18.0 |  2.7  3.2 | 13.3 13.9 | 23.9 24.3
18.7 18.5 |  3.1  3.4 | 13.8 14.3 | 24.4 24.9

次に月相の図であるが, これは要するに図学の話題で, 図学は教養学部の学生であった頃の得意科目であり, 簡単の単であった.


2つの円があり, 上は月を北の極から見たもの, 下は月を地球から見たものである. 上の図の真上に太陽があるのが新月で, 図はφだけ太陽が右によった, あるいは月が東に進んだことを陰で示す. この陰の境界を地球から見た図に書けばよいわけで, h-h'の水平面で月を輪切りにし, 各断面での陰の位置を下の図へ移せば良い.

これで難問は片付いたので, これも似たようなのがインターネットにあったが, 2009年の各月の奇数日の月相を書いてみた.



7月21日と23日に新月が見える. この中間の22日に日食がある.