2016年2月24日水曜日

MathematicaのCal

いろいろなプログラム言語でCalのプログラムを書いてみた. 今回はMathematicaである.

Mathematicaというと, 一発処理のような印象である. 例えば Prime[25] と入力すると, 25番目の素数 97が出力される. しかしなんとかプログラムを書いてみようと思い, 例によってCalを書いてみた. 結局いろいろなノウハウを覚えつつ試行錯誤の結果完成したのが次のプログラムである. その見どころを説明しよう.



0行目. 仮引数はyr_のように後に下線をつける. 本体で使うときはその下線はいらない. システムの関数手続きは大文字で始めるが, ユーザの関数は小文字で始める. 本体の定義は:=の後に書く.

1行目. Module[{ の次にローカル変数を書く. ローカルな関数名も書く.

2,3行目. printd[d_]は, その月の最後の日付をzとして, 1≤d≤zの時は出力し, それ以外は3文字の空白にする. さらにdが1桁なら前に2文字の空白, 2桁なら1文字の空白をつける.

4,5行目. prline[d_]は, dからd+6までをprintdで文字列にしたものを, StringJoin(<>を使う)し, Printで出力する.

6,7行目. 年と月を出力. また曜日名を出力.

8,9行目. うるう年なら, ローカル変数lを1に, そうでなければ0にする.

10行目. 各月の朔日の曜日を決めるため, 次の年の初めまでの日数を計算.

11行目. 1,2月の日数と, 月初めから年末までの日数(mod 7)を設定.

12,13行目. 3月以降の月に上のような値を計算する.

14行目. 朔日の曜日wを確定する.

15行目. dを第1週の日曜に相当する日(1-w)から月末日まで, 7ずつ増やしてprlineする.

こうして出来たのが次だ.



なれない言語でプログラムを書くのはしんどい.

2016年2月15日月曜日

月齢カレンダー

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

それではその年の月齢定数をCとすると, その年のmd日の近似的な月齢を (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の年が回って来る.

2016年2月14日日曜日

菱形六十面体

2010年5月17日のこのブログの続きである. そこにはWolfram MathWorldのRhombic Hexecontahedronの項に 20 golden rhombohedra can be combined to form a solid rhombic hexecontahedron. と書いてあるのが気になっていて, そのrhombohedronの頂点の立体角を計算し, たしかに20個あると全立体角, 4πステラジアンになると書いた.

golden rhombohedra は菱形六面体とでも訳そうか, 6つの面が黄金菱形になっているものである.

Wolfram MathWorldのGolden Rhombohedronの項は, これには鋭角のものと鈍角のものがあるというが, ここで使うのは鋭角の方である. この項には1辺の長さaの鋭角黄金六面体の 体積は1/5√(10-2√5) a3とあり, また菱形六十面体の項にはその体積が4√(2(5+√5))a3とあるから, 体積的にも確かにそうである.

ざっと考えてみると, 菱形六十面体の母体である正十二面体には頂点が20あるから, 多分その頂点の数と菱形六面体が対応しているであろうと推察できるが, あまり自信はない.

昨今, 3Dプリンタが使えるようになったので, ちょっと作ってみようという気になった.



上の写真の上は20個の菱形六面体が出来たところ, 下はその5個を使い, 菱形六十面体の内, 正十二面体の一面に相当する面が完成した所である.

この写真では分らないが, 3Dプリンタで作ったものは, 存外 外面がざらざらで, 接着するのが大変であった. 次回はもっと時間がかかっても, もう少し精密な面にしてみたい.

さて, 次の5個で同様に正十二面体の反対側の面に相当する面も作る. この両方の5個組を北極と南極とすれば, 残りの10個で赤道に相当する部分を構成することになる.

それには10個の内2個ずつを接着して2個組を5組作り, それを一周するように繋げるのである. これは仲々面倒で, どの面とどの面が合さるのか, いろいろこねくり回さなければならない. まぁそうこうする内に固定出来た. 下の写真は接着中の写真で, 小型のクランプやしゃこ万力で押えているところである. 左に転がっている2個は北極と南極の5個組である.



出来上がった赤道部分のこのようである.



次の図はPostScriptで描いた, 菱形六十面体と上でいう北極, 赤道, 南極の図である. 手前右上から光が当っているように影を付けてみたが, 効果はいまいちであった.

全体
北極
赤道
南極

相対的な位置が分るように, 外接する正十二面体と, 中心を通るx,y,z軸も描いてある.

これで見ると20個すべての菱形六面体が中心に接しているのが分る. つまり赤道の組は中央が厚さ0である.

3Dプリンタで作ったモデルにはこのような精度は残念ながらない. でもこんな実験が出来るのも3Dプリンタのお蔭である.

2016年2月3日水曜日

菱形六十面体

菱形六十面体についてこのブログで何度も取り上げた. その最初, 2010年5月16日のに「角D'B'E'は黄金菱形なので, 黄金比φ=1.61に対して2cot-1(φ)=63.435°である. その一端を持ち上げると, 角DBEは少しずつ広がり, 72°になるわけだ. その時の傾き角を計算すると, 驚いたことに31.7175°であった. 」と書いている.


長い間 気になっていたが, 先頃 証明した. 図を見てほしい.

図の下は破線のAD'B'E'が縦1, 横φの黄金比の菱形である. 実線のADBEは角Aが72度で縦が1の菱形である. この横の対角線の長さを以下eとする. 黄金比の方の角Aの半分は31.7175°である.

図の上の直角三角形は底辺が下の実線の菱形の横の対角線, 斜線が黄金比の対角線で, これを描いてみると角Aがやはり31.7175°になるのである. この角を以下θとする.

証明したいのは, 上の31.7175°, つまりcos θ=e/φが下の31.7175°, つまりtan 2/φに等しいということである. 一連の計算は下の通り.


φの値はよく知られている. tan 36°はcos 36°がφ/2 から計算できる.

(0)上の図のθが欲しいので, atanの式を書く. 以下atanの引数を計算する.

(1) φとeの値を代入する. (2)根号の中, 第1項を展開し, 第2項の分母を有理化する. (3)根号の中を通分し, 1/eの値を書く.

(4)根号の中の分子を計算する.

(5)根号の中, 分母分子を2倍すると, 分母分子とも開平出来た.

一方 1/φも計算すると(6)(7)同じ値になった. めでたしめでたし.