2012年4月22日日曜日

再帰曲線

東京大学新聞に2012年度後期日程入試問題が掲載されていた. その総合科目IIにフラクタルの絵があった. 問題はフラクタル次元を計算するものだが, 私としてはその絵を描くプログラムに関心があった.

まず図はこういうものだ.


左からF1, F2, F3


左からG1, G2, G3

問題に記述はこうだ.

はじめに, 1辺の長さ1の正方形をF0とする. F0を, 1辺の長さ1/3の9個の小正方形に分割し, 中央の小正方形を取り去ったものをF1とする.  また, それぞれの小正方形を1辺の長さ1/3のユニットと呼ぶ. 次に, F1 を構成する8個のユニットのそれぞれについて, これを9個の小正方形に分割し, そこから中央の小正方形を取り去ったものをF2とする.  また, この分割で得られたそれぞれの小正方形を1辺の長さ1/9のユニットと呼ぶ. 以下, 同様の操作を繰り返して得られる図形をF3, F4, ...とする.

PostScriptのこのプログラムは次のようだ.
/n 5 def
/l 540 def
/xs [0 1 2 0 2 0 1 2] def
/ys [0 0 0 1 1 2 2 2] def
/a {n 0 eq {0 0 moveto l 0 rlineto 0 l rlineto l neg 0 rlineto
 closepath fill}
 {/n n 1 sub def
 gsave 1 3 div dup scale
 0 1 7 {/i exch def
 xs i get l mul ys i get l mul gsave translate a grestore} for
 grestore
 /n n 1 add def} ifelse} def
 40 40 translate a

まず全体の次数nを決める. この例では5. 1辺の長さlも決める. 次に8個の正方形を描くので, その順に左下のxとy座標のリストxs, ysを書く. そしてユニットを描く手続きaの定義だ.



n=0なら, 長さlの正方形を書く. そうでないならnを1減らし, スケールを1/3にする. dupはx座標とy座標の両方を1/3にするためのもの. gsaveとそれに見合うgrestoreは, スケールや原点を変えるとき, 以前の環境をスタックするものである. xsとysから各小正方形の原点の座標をとり, 図を各手続きaを呼ぶ.

こうしてn=5を描いたのがこれだ.



次のは少し手ごわい.

はじめに, 1辺の長さ1の正方形G0を, 1辺の長さ1/2の正方形ユニット4個に分け, 右上のユニットについては, 左下の1辺1/4の正方形を取り去る. こうして得られた図形をG1とする. 次に, 残った3個の, 欠陥のないユニットのそれぞれについて, 同様の操作を行なう. すなわち, それぞれを4個の1辺の長さ1/4のユニットに分け, 右上にユニットについては左下の1辺の長さ1/8の正方形を取り去る. 得られた図形G2とする. 続いて, G2に含まれる1辺の長さ1/4の, 欠陥のない各正方形ユニットについて, 同様の操作を行なう. 得られた図形をG3とする. これを繰り返して得られる図形の列をG1, G2, G3, ...とする.

目で見ると簡単だが, 欠陥のない正方形をどう判定するか.

私の考え方はこうだ.

各正方形を4区画に分け, 左下, 右下, 左上, 右上の順に下請けに渡すとする. G0を描く手続きをaとすると, 1/2のサイズで, a,a,aと3回呼び, 最後に左下を取り除いた正方形を描く手続きbを呼ぶ.

aはnを1減らして, a,a,a,bと呼べば良い. bはnを1減らし, 左下を除いて, a,a,aと呼べば良い. そう考えて書いたのが, このプログラムだ. 上のプログラムが理解出来れば, こちらも分かるだろう.

/n 7 def
/l 540 def /l2 l 2 div def
/a {n 0 eq{0 0 moveto l 0 rlineto 0 l rlineto l neg 0 rlineto
 closepath fill}
 {/n n 1 sub def
 gsave
 0.5 dup scale
 gsave 0 0 translate a grestore
 gsave l 0 translate a grestore
 gsave 0 l translate a grestore
 gsave l l translate b grestore
 grestore
 /n n 1 add def} ifelse} def
/b {n 0 eq{
 l2 0 moveto l2 0 rlineto 0 l rlineto l neg 0 rlineto
 0 l2 neg rlineto l2 0 rlineto closepath fill}
 {/n n 1 sub def
 gsave 0.5 dup scale
 gsave l 0 translate a grestore
 gsave 0 l translate a grestore
 gsave l l translate a grestore
 grestore
 /n n 1 add def} ifelse} def
 40 40 translate a


この結果, n=7で描いたのが, 下の図である.



PostScripでこういう図を描くのは思ったより易しい.

2012年4月11日水曜日

太陽太陰暦

英語ではlunisolar calendarというらしいから, 日本語とは順が反対だが, 太陰太陽暦でもいいようだ.

要するに月齢に合わせて月を決めるが, ときどき閏月をおいて, 季節を合わせる. その規則はどこにでも見つけることができる.

• 新月の起きる日を, 新しい月の朔日(ついたち, 1日)とする.

• 24節気のうち, 太陽の黄経が30度の倍数になる中気に各月を対応させる. 中気の対応しない月を閏月とする.

昔の本を読むと, 春一月, 夏四月, 秋七月, 冬十月と書いてあるが, 旧暦では1, 2, 3月が春, 4, 5, 6月が夏, 7, 8, 9月が秋, 10, 11, 12月が冬である. その各月に中気が対応し, 中気のない月が閏と書いてあるから, 手元の天文年鑑を見ながら計算してみた.

新月と節気の日付と時刻を順に書いてみると,

新月 2月22日 7h35m --- a月1日(2月)
啓蟄 3月 5日 13h21m 345度
春分 3月20日 14h14m 0度
新月 3月22日 23h37m --- b月1日
清明 4月 4日 18h 6m 15度
穀雨 4月20日 1h12m 30度
新月 4月21日 16h18m --- c月1日
立夏 5月 5日 11h20m 45度
小満 5月21日 0h16m 60度
新月 5月21日 8h47m --- d月1日
芒種 6月 5日 15h26m 75度
新月 6月20日 0h 2m --- e月1日(5月)
夏至 6月21日 8h 9m 90度

aの月は春分があるから2月である. eの月は夏至があるから5月だ.




b,c,dが3月と4月とその閏月になる. bの月には黄経30度の穀雨があるから3月になる. cの月は小満があるから4月で, dの月には30度の倍数になる節気がないから, 閏4月だと思った.

ところが旧暦に関するウェブページを見てみると違うのである. 閏は4月21日からの3月の方であって, dの月が4月であった.

なぜかというと, 小満と5月の新月は, 小満の方が7時間半ほど早いが, 同じ日である. 新月が何時でも, その日の午前0時から5月になるのである. 従って小満はd月1日にあることになっているらしい.

天体の運行に基づくカレンダーはこういうところが微妙だ. 平朔とかカレンダー用のアルゴリズムをCalendrical Calculationsでもっと勉強する必要がある.