楕円アルゴリズム

MITのHakmemには意外な話題がたくさんある. プロッタで「円を描くアルゴリズム」もその1つだ.

newX=oldX-epsiron*oldY
newY=oldY+epsiron*newX

下の式で, newXを使うところがなんとなく怪しげなところ. さっそくPostscriptでやってみる.

200 200 translate
/x 100 def /y 0 def /eps 0.01 def
x y moveto
1 1 628 {pop
/x x y eps mul sub def  %x=x-eps*y
/y y x eps mul add def  %y=y+eps*x
x y lineto
} for stroke

とやるとちゃんと円が描ける.

一方, yの計算に昔のxの値を使うには, Postscriptのスタックを活用し

/x x y eps mul sub /y y x eps mul add def def

とすると, おかしいことになる. やはり新しいxが正解である.

epsの値を徐々に大きくすると, 円ではなく, 45度傾いた楕円を描くプログラムだったことが分かる. そこで楕円の長軸と短軸の長さを計算し重ねてみた.

200 200 translate
/x 100 def /y 0 def /eps 0.5 def
x y moveto
1 1 6.28 eps div {pop
/x x y eps mul sub def
/y y x eps mul add def
x y lineto
} for stroke
gsave
1 0 0 setrgbcolor       %赤にする
45 rotate               %45度傾ける
1.154700 0.894427 scale %長軸と短軸の長さにscaleする
100 0 moveto 0 0 100 0 360 arc stroke
grestore
0 0 1 setrgbcolor
-100 0 moveto 100 0 lineto 
0 -100 moveto 0 100 lineto stroke

(100,0)から出発したときのepsに対する長軸と短軸の長さと離心率を数値計算で求めると

eps   長軸   短軸  離心率
0.5  115.47 89.44  0.6325
0.1  102.60 97.59  0.3086
0.01 100.25 99.75  0.0997

のようになる. 要するに「円を描くアルゴリズム」ではなく, 「楕円を描くアルゴリズム」であった.