点線で楕円を描く・その2

[追記] この記事には誤りが含まれています。訂正した記事はこちらです。

昨日の記事の続きです。

昨日は計算式までは出したのですが、実際に計算せずにはいられなかったのでやってみました。

Numpy とかあるのは知ってるのですが、標準の Python だけでどこまでできるかやってみました。
描画は gnuplot 使ってやってます。

ソースコード

ちゃちゃっと書いたので、コメントとかは省いてます。

import sys
import math

def delta(phi, du, k2):
    return phi + du / math.sqrt(1 - k2 * math.sin(phi)**2)

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

    k2 = 1 - (b/a)**2
    phi = 0

    with open(filename, 'w') as out_file:
        print >> out_file, '#x\ty'
        i = 0
        while phi <= 2 * math.pi:
            i += 1
            i %= 10
            if i < 5:
                print >> out_file, '{0}\t{1}'.format(a * math.cos(phi), b * math.sin(phi))
            else:
                print >> out_file, ''
            phi = delta(phi, du, k2)

                

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

コマンド

$ python ellipse.py ellipse.dat 1.0 0.5 0.01

gnuplot> plot "ellipse.dat" w l

plot 結果

あとがき

思ったより誤差が発散しなかったので楽でした。

今日ばっかりは俺ってすげーって思いました :P

それでは。