私は暦が好きなので, Shakespeareの劇を読んでいても「暦」の文字があると立ち止ってしまう.
例えば, 誰でも思い出すであろう「夏の夜の夢」3幕1場で, 6人の職人が森の中で芝居の稽古をする場面:
スナウト おれたちが芝居をやる晩に月は出るのか?
ボトム 暦だ, 暦だ! 今年一年の暦を見て, 月が出るか調べるんだ, 月は出るか.
クインス うん, その晩月は出る.
(ここでの引用は, 私と同時代人である小田島さんの訳による.)
旧暦の国の人なら, 日付けで月のありなしは分る.
1日の前後は月出は夜明, 月入は夕方
8日の前後は月出は夜中, 月入は昼間
15日の前後は月出は夕方, 月入は夜明
22日の前後は月出は昼間, 月入は夜中
である. だから明治5年12月2日, 旧暦の最後の日は, 地球照で鎌のように細く光る三日月が
西の空にあったに違いない.
英国は, というかヨーロッパの殆んどは太陽暦だったから, 月の様子を知るには
暦が必要だったらしい.
しかしまだ疑問がある. 今我々の手元にある普通の暦は, 日付と曜日の対応を示
すのが主目的で, 次に祝日, 二十四節気, 雑節, さらに詳しいのには大安や友引の六曜の日付
も記載がある程度. 日出, 日入, 月出, 月入の時刻や月相は殆んど書いて
ない. 天体の情報があるのは天文暦almanacだが, 英文「夏の夜の夢」の
当該部分はcalendarになっている.
東京の毎日の日出入, 月出入は, 理科年表の暦部や
暦象年表あるいは新聞の隅をよく見なければならない.
ところが1幕1場では:
シーシュース 楽しい日々をあと四日すごせば
新月の宵となる.
ヒポリタ そうすれば新月が, 天高く引きしぼられた
銀の弓のようにかかり, 私たちの婚儀の夜を
見まもってくれるでしょう.
だから, いつが新月か知る人は知っていた. 婚儀の日は新月なのに
「うん, その晩月は出る」というShakespeareの無頓着ぶりも楽しい.
ユリウス暦を採用した張本人「ジュリアス・シーザー」2幕1場では:
ブルータス そうだ, 明日は三月十五日であったかな?
ルーシャス さあどうでしたか.
ブルータス 暦を調べて, すぐ知らせてくれぬか.
ルーシャス かしこまりました.
ルーシャス 旦那様, 三月は十四日すんだところです.
正確に管理された日めくりカレンダーでなければ, 今日が何日かを知るの
も普通には難しい. 最近ならスマホが手元にあれば分るが... さらにローマの暦の
日付は非常にややこしいから, 何日経過したかを彼等はどう計算したのか
興味がある.
ローマ暦のややこしい日付けに関しては
この解説を見て欲しい.
その3幕1場で:
シーザー 三月十五日がきたぞ.
占い師 きましたが, シーザー, まだすぎさってはいません.
この「まだすぎさってはいません」がずっしり重い.
Shakespeareの史劇に登場する薔薇戦争時代で最も有名で残忍で話題も多いのは
「リチャード三世」であろう. 最近でも遺骨が発見されて改葬されたり, エドワード四世
の2人の王子はロンドン塔で殺害されたのではないというニュースがあったりした.
「リチャード三世」3幕8場:
リチャード おい, 暦をとってくれ.
だれか今日太陽を見たものはおるか?
これは日出の時刻を知ろうとしたものらしい.
「冬物語」の4幕4場, ごろつきのオートリカスが道化に, お前は酷い目に会うぞと驚かす場面
では:
オートリカス 一年中で一番暑い日と暦に出ている日に煉瓦塀にもたせかけられ, 太陽の
カッと燃える目に見つめられながら, ...
7月22日に大暑と書いてあるようなことをオートリカスはいったのだろう
か. 理科年表の「気温の最高および最低記録」のようなものを指すのか.
ついでだが, 昔の人は時刻を何で知ったか? ヨーロッパの町にいると, 教会の
鐘の音で「ヨーロッパに来たなぁ」と思う. でも教会から遠い戦場では, 太陽の
方向以外に時刻を知るのは難しかったであろう.
ロミオとジュリエット 1幕1場:
ベンヴォーリオ おはよう, ロミオ.
ロミオ 早いのか, まだ?
ベンヴォーリオ 九時だよ.
近くの教会の鐘が鳴ったか.
再び暦に戻って:
「マクベス」 4幕1場, マクベスといい合った魔女や幻影が消えた後:
マクベス どうした? 行ってしまったか? 畜生, この忌ましい時め, 暦の
つづくかぎり呪われるがいい!
自然の「天地の続く限り」ではなく, 人工の「暦の続く限り」というのが
意外な発想である. Shakespeareの暦についての考えが少しは分ったように思う.
2025年3月14日金曜日
2025年1月26日日曜日
月齢カレンダー
私がこの数年愛用している予定表の話はこのブログに何回か書いた. ところで, 今年の予定表の表紙にある週日定数も月齢定数も共に0である. つまり下のようなのだ.
そうするとこの次にこうなるのは何年後かも気になる.
y年の週日定数, week(y)と月齢定数, moon(y)は次の式で計算する.
この表を見た序でに, 今年の復活祭を調べよう. まず3月21日の月齢は21. 次の満月は月齢45として4月14日. その日の曜日は1で月曜だから, 日曜はその6日後で4月20日になる.
もうひとつ序でに. Winning WaysのDoomsdayは3月0日の曜日なので, 今年は5, つまり金曜日である.
そうするとこの次にこうなるのは何年後かも気になる.
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, つまり金曜日である.
2024年8月4日日曜日
Heronの式
前回のブログ「四面体の体積」に名称だけ出てきたHeronの式というのがある.
三角形の3辺の長さをそれぞれa, b, cとし, (a+b+c)/2をSとすると, その三角形の
面積は√S(S-a)(S-b)(S-c)というのである.
たとえば, 辺の長さが3,4,5のPythagoras三角形ではS=6だから, 平方根号の中は 6×3×2×1=36で, 面積は6だ. 右の図, 1辺が2の正三角形では, S=3で, 根号内は3×1×1×1=3で, 面積は√3である.
Heronの式を覚えた時, これが正しいという証明を見たかどうかは分らない. 勿論, Sの次元は長さ, S-a,S-b,S-cも長さなので, 根号内の次元は長さの4乗; 開平すると 2乗になり, つまり面積だ. そういえば, 前回のブログにあった, 四面体の体積を6本の 辺の長さから求める式も体積の次元になっている. あたりまえだ.
でも定量的に合っているかを調べるのが今回の頭の体操である. 方針としては, 根号の あるのがウンシャンだから, 面積の二乗の値を考えることにする. 例題の三角形は こういう頂点名と辺名を持つ. 計算の途中で時々数値的にチェックするのにもこれを 使う. この三角形は∠Aが直角なので, 面積は√8×√18/2=6だ.
式を使わざるを得ないから, Texの出力のpdfにも助けてもらおう.
式(0)はHeronの式で, 2乗すると(1). Sを定義で戻し, 展開したのが式(2)である. この程度 の計算はWolframAlphaにやって貰う.
三角形ABCの各頂点の対辺をa,b,cとすると, その各々の2乗もpdfにある通り. これを 式(2)の分子に入れて展開するわけだが, これは中々面倒である. WolframAlphaも 時間がかかるとか文句をいうし, 我が家のRapsberry PiのMathematicaは今は ネットに繋っていないので, 使い勝手が悪い. 結局は多項式の乗算のプログラムを Schemeでちょこっと作り, それで展開した. 時々例題の値を入れて験算する.
そして展開した結果は項数21のpdfにある式である. これも目がくらくらするから, xixjykylについて, ijと klの対の項の係数を表にしてみた. さらに例題の座標で験算すると576. 16で 割ると36だ. しめしめ.
さてこの式は本当に三角形の面積であろうか. 頂点の座標(x0y0), (x1y1),(x2y2)の三角形の 面積(頂点の回り方によってはその負数)の式は, 幸いにも私は高校生の頃から知っている. 次のpdfにある式だ. 行列式を計算式にしてその分子の2乗を展開すると, またもや項数21の多項式が得られた. 今回もその係数の表を作るとブラボー! 先の表の1/4の値である. という訳けで Heronの式は面積になるのであった.
だが待てよ. この行列式はなぜ正しいのか. こちらも気になったので, 計算してみた.
図の三角形の, 座標軸に平行な線での外接長方形を作る. それから長方形の辺と 三角形の辺で作る3個の三角形の面積を引く. そういう式を作ってみる.
WoflramAlpha様に
Expand[(x2-x1)(y0-y2)-(x0-x1)(y0-y1)/2-(x2-x1)(y1-y2)/2-(x2-x0)(y0-y2)/2]
と伺いを立てたら
(x0y1)/2-(x0y2)/2-(x1y0)/2+(x1y2)/2+(x2 y0)/2-(x2y1)/2
のお告げが出た. 整理すると
(x0y1+x1y2+x2y0-x0y2-x1y0-x2y1)/2
となり, 行列式の値と一致している. もう一度ブラボー! ちゃんとした証明でもないが, 私にとっては十分納得出来る頭の体操であった.
たとえば, 辺の長さが3,4,5のPythagoras三角形ではS=6だから, 平方根号の中は 6×3×2×1=36で, 面積は6だ. 右の図, 1辺が2の正三角形では, S=3で, 根号内は3×1×1×1=3で, 面積は√3である.
Heronの式を覚えた時, これが正しいという証明を見たかどうかは分らない. 勿論, Sの次元は長さ, S-a,S-b,S-cも長さなので, 根号内の次元は長さの4乗; 開平すると 2乗になり, つまり面積だ. そういえば, 前回のブログにあった, 四面体の体積を6本の 辺の長さから求める式も体積の次元になっている. あたりまえだ.
でも定量的に合っているかを調べるのが今回の頭の体操である. 方針としては, 根号の あるのがウンシャンだから, 面積の二乗の値を考えることにする. 例題の三角形は こういう頂点名と辺名を持つ. 計算の途中で時々数値的にチェックするのにもこれを 使う. この三角形は∠Aが直角なので, 面積は√8×√18/2=6だ.
式を使わざるを得ないから, Texの出力のpdfにも助けてもらおう.
式(0)はHeronの式で, 2乗すると(1). Sを定義で戻し, 展開したのが式(2)である. この程度 の計算はWolframAlphaにやって貰う.
三角形ABCの各頂点の対辺をa,b,cとすると, その各々の2乗もpdfにある通り. これを 式(2)の分子に入れて展開するわけだが, これは中々面倒である. WolframAlphaも 時間がかかるとか文句をいうし, 我が家のRapsberry PiのMathematicaは今は ネットに繋っていないので, 使い勝手が悪い. 結局は多項式の乗算のプログラムを Schemeでちょこっと作り, それで展開した. 時々例題の値を入れて験算する.
そして展開した結果は項数21のpdfにある式である. これも目がくらくらするから, xixjykylについて, ijと klの対の項の係数を表にしてみた. さらに例題の座標で験算すると576. 16で 割ると36だ. しめしめ.
さてこの式は本当に三角形の面積であろうか. 頂点の座標(x0y0), (x1y1),(x2y2)の三角形の 面積(頂点の回り方によってはその負数)の式は, 幸いにも私は高校生の頃から知っている. 次のpdfにある式だ. 行列式を計算式にしてその分子の2乗を展開すると, またもや項数21の多項式が得られた. 今回もその係数の表を作るとブラボー! 先の表の1/4の値である. という訳けで Heronの式は面積になるのであった.
だが待てよ. この行列式はなぜ正しいのか. こちらも気になったので, 計算してみた.
図の三角形の, 座標軸に平行な線での外接長方形を作る. それから長方形の辺と 三角形の辺で作る3個の三角形の面積を引く. そういう式を作ってみる.
WoflramAlpha様に
Expand[(x2-x1)(y0-y2)-(x0-x1)(y0-y1)/2-(x2-x1)(y1-y2)/2-(x2-x0)(y0-y2)/2]
と伺いを立てたら
(x0y1)/2-(x0y2)/2-(x1y0)/2+(x1y2)/2+(x2 y0)/2-(x2y1)/2
のお告げが出た. 整理すると
(x0y1+x1y2+x2y0-x0y2-x1y0-x2y1)/2
となり, 行列式の値と一致している. もう一度ブラボー! ちゃんとした証明でもないが, 私にとっては十分納得出来る頭の体操であった.
2024年8月1日木曜日
四面体の体積
筑波大の三谷さんのツィッターに四面体の体積を計算する話があった.
角封筒の上の辺の中央と下の辺の両端を結ぶ折目をつけ, それを山折りにすると, 合同な 二等辺三角形4枚による四面体が出来る. 元の封筒の横と縦が1と√2の時の 体積を計算せよという問題である. ちょっとやってみたところ, 存外簡単であった.
下の図で, 赤い線は元の封筒を展開したもの. 下の辺と上の辺は繋っている. この斜めの 線を谷折りにすると, ABCを底面, ACD'を手前の面, ABDを向うの面, 辺DEとD'E'は 繋っていて斜面だが上の面になる.
錐体の体積を計算するには, 底面積と高さが必要だが, 底面積は簡単だ. 辺BCが1, 辺OAは辺ABが√2だから, √7/2である.
高さについては, 三角形ABDは, 辺ABを折目として頂点Dが起き上がるから, Dに対応する 平面上の点はABと直交する破線で示す線に沿って下がり, AO上のGまで来ると下から上って来たD’と 出会って頂点になる. GAの距離はADの距離にtan αを掛けたものだが, tan αはOB/OAだしADは1だから, GA=1/√7である. よって 高さは√6/√7.
求める体積は (√7/2)*(√6)/((√7)*2*3)=1/2√6
三谷さんのツィッターにはヒントと称して, 四面体の6辺の長さからその体積を 求める式があった. 平面三角形のHeronの公式の三次元版である. そこにある式の 添字の付け方を, 私流に修正したのが次の図である.
このように添字を付けると, 四面体の体積Vの2乗は の一番上の式で得られる. しかしこれも目がくらくらするので, 始めの3項と終りの4項を それぞれ6変数のtfと3変数のtgに分け, 2行目以下のように定義する.
最初のV2を計算する式を書いたLatexは次の様だ. 計算式と同様に分解してある.
ところで, 最初の封筒の四面体の体積は, 形状が特殊であったせいか, 簡単であった. 底面の3点の座標と, その各点から頂点までの距離が与えられている時にも体積は 得られるだろうか. 以下はその計算である. 式が沢山あるので, Latexでpdfに したもので示す.
底辺の3点, A, B, Cの座標をそれぞれ(x0,y0,0),
(x1,y1,0),(x2,y2,0);
頂点Dの座標を(x,y,z)とする. A,B,Cから頂点までの距離をそれぞれ
l0,l1,l2とする.
今回もpythonで計算すると
角封筒の上の辺の中央と下の辺の両端を結ぶ折目をつけ, それを山折りにすると, 合同な 二等辺三角形4枚による四面体が出来る. 元の封筒の横と縦が1と√2の時の 体積を計算せよという問題である. ちょっとやってみたところ, 存外簡単であった.
下の図で, 赤い線は元の封筒を展開したもの. 下の辺と上の辺は繋っている. この斜めの 線を谷折りにすると, ABCを底面, ACD'を手前の面, ABDを向うの面, 辺DEとD'E'は 繋っていて斜面だが上の面になる.
錐体の体積を計算するには, 底面積と高さが必要だが, 底面積は簡単だ. 辺BCが1, 辺OAは辺ABが√2だから, √7/2である.
高さについては, 三角形ABDは, 辺ABを折目として頂点Dが起き上がるから, Dに対応する 平面上の点はABと直交する破線で示す線に沿って下がり, AO上のGまで来ると下から上って来たD’と 出会って頂点になる. GAの距離はADの距離にtan αを掛けたものだが, tan αはOB/OAだしADは1だから, GA=1/√7である. よって 高さは√6/√7.
求める体積は (√7/2)*(√6)/((√7)*2*3)=1/2√6
三谷さんのツィッターにはヒントと称して, 四面体の6辺の長さからその体積を 求める式があった. 平面三角形のHeronの公式の三次元版である. そこにある式の 添字の付け方を, 私流に修正したのが次の図である.
このように添字を付けると, 四面体の体積Vの2乗は の一番上の式で得られる. しかしこれも目がくらくらするので, 始めの3項と終りの4項を それぞれ6変数のtfと3変数のtgに分け, 2行目以下のように定義する.
\newcommand{\sq}[1]{l_#1^2}
\newcommand{\tf}[6]{\sq#1\sq#2(\sq#3+
\sq#4+\sq#5+\sq#6-\sq#1-\sq#2)}
\newcommand{\tg}[3]{-\sq#1\sq#2\sq#3}
$V^2=\frac{1}{144}[\tf012345+\tf234501\\
+\tf450123\tg025\tg034\tg124\tg135]$
これを使ってpythonで計算したのが次だ.
import math
def sq(x):
return x*x
def tf(a,b,c,d,e,f):
return sq(a)*sq(b)
*(sq(c)+sq(d)+sq(e)+sq(f)-sq(a)-sq(b))
def tg(a,b,c):
return -sq(a)*sq(b)*sq(c)
def v2(l0,l1,l2,l3,l4,l5):
return (tf(l0,l1,l2,l3,l4,l5)+tf(l2,l3,l4,l5,l0,l1)
+tf(l4,l5,l0,l1,l2,l3)
+tg(l0,l2,l5)+tg(l0,l3,l4)+tg(l1,l2,l4)
+tg(l1,l3,l5))/144
l0=1;l1=1;l2=math.sqrt(2);l3=math.sqrt(2)
l4=math.sqrt(2);l5=math.sqrt(2)
print(math.sqrt(v2(l0,l1,l2,l3,l4,l5)))
=> 0.20412414523193143
この値は確かに1/2√6 である. ところで, 最初の封筒の四面体の体積は, 形状が特殊であったせいか, 簡単であった. 底面の3点の座標と, その各点から頂点までの距離が与えられている時にも体積は 得られるだろうか. 以下はその計算である. 式が沢山あるので, Latexでpdfに したもので示す.
x0=0;y0=-1/2;x1=math.sqrt(7)/2;y1=0;x2=0;y2=1/2 l0=math.sqrt(2);l1=1;l2=math.sqrt(2) a0=2*(x0-x1); b0=2*(y0-y1) c0=sq(l1)-sq(l0)+sq(x0)-sq(x1)+sq(y0)-sq(y1) a1=2*(x1-x2); b1=2*(y1-y2) c1=sq(l2)-sq(l1)+sq(x1)-sq(x2)+sq(y1)-sq(y2) x=(c0*b1-c1*b0)/(a0*b1-a1*b0) y=(a0*c1-a1*c0)/(a0*b1-a1*b0) z=math.sqrt(sq(l0)-sq(x0-x)-sq(y0-y)) print(x,y,z) #=> 0.9449111825230684 0.0 0.9258200997725514期待通りの値が得られた. でもこんな実数の値より, 1/2√6 の方が嬉しい.
2024年7月2日火曜日
Fixed Day Number
前回のブログのアルゴリズムに, 年初から前月末までの総日数を計算する式:
のような図になる. どちらも2月は30日あるとする.
まず左の図に注目する. 1月から2月は下がり, 2月から3月は上がり, 4月へは下がり, ... 8月と9月へは下がり, ... つまり大の月の後は下がり, 小の月の後は上がるらしい.
計算して見ると下がる量は0.4167, 上がる量は0.5833である.
もう一度式を見ると, monthが1増えると, 367 / 12 つまり 30.5833 増える. これと30及び31の差であった. これで納得できるが, Calendrical Calculation にはこういう説明があった.
第1項はm = 1, 2, ..., 12について,
そこで上の第1項と第2項を足すと, (367 * m - 362) // 12になるのであった.
2月を30日にしたという発想がよかったのだ. 一方, 年末から計算する方の式はこう導く.
欲しい値は1月, 2月, ..., 12月で, 2月も30日とするから
7-(7 * m - 2) // 12
でこれを(13 - m) * 30 に足す. これを一纏めにしたい.
まず下の計算の1行目のように, 床関数を天井関数にする. そして1を足して床関数に戻す. しかしmの値によっては, 整数になることがあり, 1を足すと失敗する. 従って, 2行目のように, 1の代りに 11/12を足してみる. 3行目はそれに30日の倍数を足したところ. 最後に分数の外に ある部分を分子に足す. すると(4777-367m)//12が得られる.
どうだろうか.
(367 * month - 362) // 12があった. monthを1から12まで変えて値を計算すると.
[0, 31, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336]が得られる. 演算子「//」で小数点以下を落しているが, その落された値を調べると, ( (4777 - 367 * month) // 12 の図も右に示す.)
のような図になる. どちらも2月は30日あるとする.
まず左の図に注目する. 1月から2月は下がり, 2月から3月は上がり, 4月へは下がり, ... 8月と9月へは下がり, ... つまり大の月の後は下がり, 小の月の後は上がるらしい.
計算して見ると下がる量は0.4167, 上がる量は0.5833である.
もう一度式を見ると, monthが1増えると, 367 / 12 つまり 30.5833 増える. これと30及び31の差であった. これで納得できるが, Calendrical Calculation にはこういう説明があった.
(7 * m - 2) // 12 + 30 * (m - 1)この第2項 30 * (m - 1) はすべて小の月として前月末までの総日数.
第1項はm = 1, 2, ..., 12について,
0, 1, 1, 2, 2, 3, 3, 4, 5, 5, 6, 6 31 30 31 30 31 30 31 31 30 31 30である(上の行), 下の行は1月から11月までの日数で, こうして見ると, 大の月の 上で上の行の値は1増える. つまり小の月とした時の誤差を得ている.
そこで上の第1項と第2項を足すと, (367 * m - 362) // 12になるのであった.
2月を30日にしたという発想がよかったのだ. 一方, 年末から計算する方の式はこう導く.
欲しい値は1月, 2月, ..., 12月で, 2月も30日とするから
[367, 336, 306, 275, 245, 214, 184, 153, 122, 92, 61, 31]である. 30の倍数の列
[360, 330, 300, 270, 240, 210, 180, 150, 120, 90, 60, 30]に誤差の列
[7, 6, 6, 5, 5, 4, 4, 3, 2, 2, 1, 1]を足したものだ. これは7から年初からの計算の誤差を引いたのもだ. つまり,
7-(7 * m - 2) // 12
でこれを(13 - m) * 30 に足す. これを一纏めにしたい.
まず下の計算の1行目のように, 床関数を天井関数にする. そして1を足して床関数に戻す. しかしmの値によっては, 整数になることがあり, 1を足すと失敗する. 従って, 2行目のように, 1の代りに 11/12を足してみる. 3行目はそれに30日の倍数を足したところ. 最後に分数の外に ある部分を分子に足す. すると(4777-367m)//12が得られる.
どうだろうか.
2024年7月1日月曜日
Fixed Day Number
私の好きな本に, Edward M. ReingoldとNachum DershowitzのCalendrical Caluculation
がある. 随分昔に購入したもので, 手元にあるのは2004年のThird Printing.
要するに暦に関するアルゴリズムが満載で面白い. この本では暦の変換の基準は, Fixed Day Numberといい, Gregorian暦が過去に遡って使われたとして, その1年1月1日(の正子)を 第1日とするものである. いちいちfixed day numberといわず, R.D.と いうらしい. (Rata Die)
似たようのものに, -4712年1月1日正午から起算するユリウス日があるが, これは途中でJulian暦 からGregolian暦になるので変換が面倒である.
Calendrical Caluculationにある, Gregorian暦の年月日からそのR.Dを計算する アルゴリズムを見てみよう. 本書のアルゴリズムはCommon Lispで記述されているが, 今日はPythonに翻訳する.
次は前月末までの日数を計算する. 目標としては, 1月なら0. 2月なら31, 3月は 平年なら59, 閏年なら60, ... 閏年の補正は後でやることにし, 本書にあるアルゴリズム は(367 * month - 362) // 12というものだ. monthを 1から12まで変えてこの値を計算すると:
これで前年末までと前月末までの総日数が分ったから, 後はdayを足せばよい.
閏年の補正に本書ではこういう式を使う.
各々の項の値は682915, 467, -18, 4, 275, -1, 14で, R.D.は683656 である.
ところで私は何度も暦の計算のプログラムを書いたが, その年の終までの日数を計算し, その月の前までの日数を引いて日を足すのが好きだ. 上のアルゴリズムに year - 1が4回もあるのが気に食わないのである.
それで普通は年末から今月0日までの日数を引くことになるが, 大体はそこは定数の 表を利用する. しかし, Calendricalの本にあるような式が使えないかと 考えた. Calendricalの式と同様, 1,2月は特別扱いにしてもよいとして, 定数など を変えてテストした. 欲しい値は:
ここで使う定数4777は非常にクリティカルで, 4776や 4778に変えてみると, 前者は7月が小さく, 後者は2月が大きい.
年末から月始めまでの日数を今回のように式で計算すると, 実引数に0月とか13月が与え られても計算してしまう恐れがある. 表なら範囲外としてエラーに出来る. それを 心配したとしても, 日の値に範囲外が与えられる可能性はあるのだから, まぁ我慢することに しよう.
要するに暦に関するアルゴリズムが満載で面白い. この本では暦の変換の基準は, Fixed Day Numberといい, Gregorian暦が過去に遡って使われたとして, その1年1月1日(の正子)を 第1日とするものである. いちいちfixed day numberといわず, R.D.と いうらしい. (Rata Die)
似たようのものに, -4712年1月1日正午から起算するユリウス日があるが, これは途中でJulian暦 からGregolian暦になるので変換が面倒である.
Calendrical Caluculationにある, Gregorian暦の年月日からそのR.Dを計算する アルゴリズムを見てみよう. 本書のアルゴリズムはCommon Lispで記述されているが, 今日はPythonに翻訳する.
def fixeddaynumber(year, month, day):
return (365 * (year - 1) + (year - 1) // 4
- ((year - 1) // 100) + (year - 1) // 400
+ (367 * month - 362) // 12
+ (0 if (month <= 2) else
(-1 if leap_year (year) else -2))
+ day)
最初の365 * (year - 1)は前年までがすべて平年であったとしての, 前年の
最後の日までの総日数である. しかし, 4年に一度閏年があるから, その分を足す.
+ (year - 1) // 4. Gregorian暦では, 100で整除される年は平年に
するのであったからそれを引く. - ((year - 1) // 100). しかし400で
整除されるなら閏年にするからそれを足す. + (year - 1) // 400.
前年末までの総日数はこれで計算出来た. 次は前月末までの日数を計算する. 目標としては, 1月なら0. 2月なら31, 3月は 平年なら59, 閏年なら60, ... 閏年の補正は後でやることにし, 本書にあるアルゴリズム は(367 * month - 362) // 12というものだ. monthを 1から12まで変えてこの値を計算すると:
print(list(map (lambda m:(367*m-362)//12, range(1,13))))
[0, 31, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336]
3月からは, 平年なら2, 閏年なら1多いことがわかる. 従ってこの関数の
month番目を使い, monthが1, 2ならこの値そのまま, 3以降はこの値から
平年は2, 閏年は1を引いて使う.これで前年末までと前月末までの総日数が分ったから, 後はdayを足せばよい.
閏年の補正に本書ではこういう式を使う.
def leap_year (year):
return (year % 4) == 0 & ((year % 400)
in [100, 200, 300])
これで1年1月1日を計算すると:
print (fixeddaynumber (1, 1, 1)) 1なるほど. また, おととし, 開業150年と話題になった鉄道開通の日, 1872年10月14日(明治5年には日本は まだ旧暦を使っていたが, これは太陽暦に換算してある)のR.D.を計算してみる.
各々の項の値は682915, 467, -18, 4, 275, -1, 14で, R.D.は683656 である.
ところで私は何度も暦の計算のプログラムを書いたが, その年の終までの日数を計算し, その月の前までの日数を引いて日を足すのが好きだ. 上のアルゴリズムに year - 1が4回もあるのが気に食わないのである.
それで普通は年末から今月0日までの日数を引くことになるが, 大体はそこは定数の 表を利用する. しかし, Calendricalの本にあるような式が使えないかと 考えた. Calendricalの式と同様, 1,2月は特別扱いにしてもよいとして, 定数など を変えてテストした. 欲しい値は:
[[365, 366], [334, 335], 306, 275, 245, 214,
184, 153, 122, 92, 61, 31]
である. 左端が1月と2月で, かっこ内は左が平年, 右が閏年. 右端が12月.
その結果, (4777 - 367 * m) // 12がよいらしかった. つまり
print(list(map (lambda m:(4777 - 367 * m) // 12,
range(1,13))))
[367, 336, 306, 275, 245, 214, 184, 153, 122, 92, 61, 31]
1月, 2月は平年は2を引き, 閏年は1を引く. ここで使う定数4777は非常にクリティカルで, 4776や 4778に変えてみると, 前者は7月が小さく, 後者は2月が大きい.
print(list(map (lambda m:(4776 - 367 * m) // 12,
range(1,13))))
[367, 336, 306, 275, 245, 214, 183, 153, 122, 92, 61, 31]
print(list(map (lambda m:(4778 - 367 * m) // 12,
range(1,13))))
[367, 337, 306, 275, 245, 214, 184, 153, 122, 92, 61, 31]
これを利用して書いたfixed day number関数は次の通り:
def fixeddaynumber2(year, month, day):
return (365 * year + year // 4 - (year // 100)
+ year // 400 - (4777 - 367 * month) // 12
+ (0 if (month > 2) else
(-1 (year % 400 == 0)
if (year % 100 == 0) else
(year % 4 == 0) -2))
+ day)
こちらで鉄道記念日を計算したのと, 本年6月30日の値の7で割った剰余とを見ると:
print (fixeddaynumber2 (1872, 10, 14)) 683656 print ((fixeddaynumber2 (2024, 6, 30)) % 7) 0だから, R.D.を7で割った剰余が曜日になるわけだ. つまり, R.D.1の日, 1年1月1日は月曜であった.
年末から月始めまでの日数を今回のように式で計算すると, 実引数に0月とか13月が与え られても計算してしまう恐れがある. 表なら範囲外としてエラーに出来る. それを 心配したとしても, 日の値に範囲外が与えられる可能性はあるのだから, まぁ我慢することに しよう.
2022年9月20日火曜日
満月の十五夜
今年は9月10日が中秋の名月, つまりお月見であった. 幸いよく晴れて, 満月が堪能出来た.
ところで新聞が「今年の中秋の名月は満月」と報じた. 私はもちろん旧暦の15日が満月にならない方が普通で, たまには満月になることも知っていたから, この新聞記事をみて何とも思わなかったが, どの程度すれすれに満月なのかと天文年鑑を見ると満月は18時59分であった.
旧暦の15日が満月にならないのは, 新月の時点のある日を旧暦の1日とするからである. 仮に1朔望月を29.5とし, 月齢14.75を満月とすれば, ある日の朝0.25日までに新月があればその日が1日で, 15日の晩には月齢が14.75に達する. そうでないと, 月齢が14.75になるのは, 旧暦16日になる.
ところが月の公転には遅速があり, 満月は太陽と月の黄経の差が180度の時のことだから, 月齢14.75からプラスマイナス1日くらいずれ得る.
そこで新月から満月, 満月から新月の経過時間の変化を見てみたいと思った. 2013年から2022年の天文年鑑から, 新月と満月の日時を書き出し, パソコンに入力して新月と次の満月, 満月と次の新月の経過時間を計算した.
10年間にある朔望月は, 19年7閏法を考えると, 120+3.5回だから, 日時データは247, 間隔は246あった. 最大値は15.601, 最小値は13.907, 平均値は14.767.
最後の10個は
... 14.32 15.278 14.071 15.497 13.958
15.579 14.009 15.497 14.216 15.256
である.
それを絵にしたのが下だ. 赤は新月から満月, 青は満月から新月までの経過時間である. 縦軸の単位は日.
図を見ると, 時々赤と青が同じになり, つまり新月から満月までと満月から新月までが同時間になり, その中間は一方は増えて減り, 他方が減って増える.
その理由は多分こうであろう. 月の公転軌道も楕円であり, 長軸の一方が近地点, 他方が遠地点である. 近地点では公転速度が速く, 遠地点では遲いのは常識だ.
近地点と遠地点の近くに新月と満月の場所があれば, 新月から満月も満月から新月もほぼ同時間になる. しかし, 新月と満月がその中間くらいにあると, 近地点側を通る方は経過時間が短かく, 遠地点側は長い. それがこの図の謎解きであろう.
2022年の新月と満月の日の月の地心距離を天文年鑑から調べ, 新月から満月までやその逆の経過時間がどの地心距離からどの地心距離であったかの図を描いてみた. それを下に示す. 横軸は経過時間(単位は日)で, 縦軸は地心距離(単位は万km).
これを見ると, 上の推論が正しいことが分る. 左上から右下へ来る矢印は, 遠地点から近地点までの経過で14.8日くらい掛ることで, 15.7日くらいの水平の矢印は, 中間から中間への遠地点側を通る経過, 14.0日くらいの水平の矢印は, 近地点側を通る経過である.
十五夜が満月になるというのには, こういう仕掛けがあったわけだ. その後, 国立天文台暦計算室のページを見付けた.
ところで新聞が「今年の中秋の名月は満月」と報じた. 私はもちろん旧暦の15日が満月にならない方が普通で, たまには満月になることも知っていたから, この新聞記事をみて何とも思わなかったが, どの程度すれすれに満月なのかと天文年鑑を見ると満月は18時59分であった.
旧暦の15日が満月にならないのは, 新月の時点のある日を旧暦の1日とするからである. 仮に1朔望月を29.5とし, 月齢14.75を満月とすれば, ある日の朝0.25日までに新月があればその日が1日で, 15日の晩には月齢が14.75に達する. そうでないと, 月齢が14.75になるのは, 旧暦16日になる.
ところが月の公転には遅速があり, 満月は太陽と月の黄経の差が180度の時のことだから, 月齢14.75からプラスマイナス1日くらいずれ得る.
そこで新月から満月, 満月から新月の経過時間の変化を見てみたいと思った. 2013年から2022年の天文年鑑から, 新月と満月の日時を書き出し, パソコンに入力して新月と次の満月, 満月と次の新月の経過時間を計算した.
10年間にある朔望月は, 19年7閏法を考えると, 120+3.5回だから, 日時データは247, 間隔は246あった. 最大値は15.601, 最小値は13.907, 平均値は14.767.
最後の10個は
... 14.32 15.278 14.071 15.497 13.958
15.579 14.009 15.497 14.216 15.256
である.
それを絵にしたのが下だ. 赤は新月から満月, 青は満月から新月までの経過時間である. 縦軸の単位は日.
図を見ると, 時々赤と青が同じになり, つまり新月から満月までと満月から新月までが同時間になり, その中間は一方は増えて減り, 他方が減って増える.
その理由は多分こうであろう. 月の公転軌道も楕円であり, 長軸の一方が近地点, 他方が遠地点である. 近地点では公転速度が速く, 遠地点では遲いのは常識だ.
近地点と遠地点の近くに新月と満月の場所があれば, 新月から満月も満月から新月もほぼ同時間になる. しかし, 新月と満月がその中間くらいにあると, 近地点側を通る方は経過時間が短かく, 遠地点側は長い. それがこの図の謎解きであろう.
2022年の新月と満月の日の月の地心距離を天文年鑑から調べ, 新月から満月までやその逆の経過時間がどの地心距離からどの地心距離であったかの図を描いてみた. それを下に示す. 横軸は経過時間(単位は日)で, 縦軸は地心距離(単位は万km).
これを見ると, 上の推論が正しいことが分る. 左上から右下へ来る矢印は, 遠地点から近地点までの経過で14.8日くらい掛ることで, 15.7日くらいの水平の矢印は, 中間から中間への遠地点側を通る経過, 14.0日くらいの水平の矢印は, 近地点側を通る経過である.
十五夜が満月になるというのには, こういう仕掛けがあったわけだ. その後, 国立天文台暦計算室のページを見付けた.
登録:
投稿 (Atom)















