点線で楕円を描く・まとめ - cocoatomo衝動日記〈移行後版〉

【追記・2011/01/22】最初に書いた内容に大幅に誤りがあったので, 再度編集しています. おっちょこちょいですみません.

【追記・2011/02/03】修正完了. ソースコードを色々変えました.

昨日、一昨日の記事には数式に間違いがあったので訂正したまとめ記事を書きます。

【修正・2013/03/12】数式に対するコメントの通り, 修正しました.

楕円の弧長

まず、 $x = a \cos \varphi,\quad y = b \sin \varphi$ とパラメータ表示して、点 (a, 0) から左回りに弧長 $u(\varphi)$ を測っていくことにします。

$\begin{align*} u &= \int_0^{\varphi} \sqrt{dx^2 + dy ^2} \\ &= \int_0^{\varphi} \sqrt{a^2 \sin^2 \varphi + b^2 \cos^2 \varphi} d\varphi \end{align*}$

となり、この逆関数 $\varphi (u)$ の微分を考えます。(ここの時点で式の形に間違いがありましたね.)

$\frac{d \varphi}{du} = \frac{1}{\frac{du}{d \varphi}} = \frac{1}{\sqrt{a^2 \sin^2 \varphi + b^2 \cos^2 \varphi}}$

これを差分形にして、

$\begin{align*} \Delta \varphi &= \frac{\Delta u}{\sqrt{a^2 \sin^2 \varphi + b^2 \cos^2 \varphi}} \\ \varphi_0 &= 0 \\ \varphi_n &= \varphi_{n - 1} + \frac{\Delta u}{\sqrt{a^2 \sin^2 \varphi_{n-1} + b^2 \cos^2 \varphi_{n-1}}} \end{align*}$

$\sin$ と $\cos$ を間違えていましたね。

ここで $a$ と $b$ の大小については何も言っていないことに注意してください。

これを実装します。

Python ソースコード

開発をローカルから Github に移したので, 最新のソースはこちらを見てください.

import sys
import math

DIVISION = 1000.0
CYCLE = 10

def angles(du, a, b):
    phi = 0
    while phi <= 2 * math.pi:
        yield phi
        phi += du / math.sqrt((a * math.sin(phi))**2 + (b * math.cos(phi))**2)

def coordinate(du, a, b):
    for angle in angles(du, a, b):
        yield (a * math.cos(angle), b * math.sin(angle))

def circumference(a, b):
    """Return a length of a circumference of an ellipse.

    @param a, b length of two semi-axes

    reference: http://en.wikipedia.org/wiki/Ellipse#Circumference
    """
    expr = ((a - b)/(a + b))**2
    return math.pi * (a + b) * (1 + (3 * expr)/(10 + math.sqrt(4 - 3 * expr)))

def draw(argv=None):
    if not argv:
        argv = sys.argv
    filename = argv[1]
    a = float(argv[2])
    b = float(argv[3])
    du = circumference(a, b) / DIVISION

    with open(filename, 'w') as out_file:
        print >> out_file, 'plot "-" w l'
        print >> out_file, '#x\ty'
        for i, coord in enumerate(coordinate(du, a, b)):
            i %= CYCLE * 2
            if i < CYCLE:
                print >> out_file, '{0}\t{1}'.format(*coord)
            else:
                print >> out_file, ''
        print >> out_file, 'end'
        print >> out_file, 'pause -1'


if __name__ == '__main__':
    draw(sys.argv)

コマンド

$ python ellipse.py ellipse.dat 1.0 0.5
gnuplot> plot "ellipse.dat" w l

出力結果

以下のようにちゃんと点線が等長、等間隔で描けています。(実際にプログラムで点線の長さと間隔を測ったところ、1 % 程度の誤差はありましたがだいたい同じでした。)

$a$ と $b$ には大小の制限は無いので縦長の楕円も描けます。

あとがき

昨日一昨日の記事を書いてから、色々な方からコメントやツッコミをいただきました。ありがとうございます。

今まで何かを作ったり、記事を書いたりして公開したことはあったけれど、今回もらった反応は今までで一番多いものでした。
そして、何かを発表して反応をもらう喜びを強く感じました。

これは病み付きになりそうです。

もしかして、オープンソースで活動してる人の中には、この喜びを原動力にやってる人もいるのかなぁ? と思ったり。

それでは、お休みなさい。