[pukiwiki]
歯車にエンピツを差してグルグルと回すと不思議な図形が描ける不思議な定規　「スピログラフ」をPythonで再現してみました。

”＜＜追記＞＞”テストしたら、数式が間違ってました
が、この模様もきれいなので、こちらも残しておきますー

-[[ググる:スピログラフ]]

Pythonネタは久しぶりです。

で、定規だと外側の歯車と内側の歯車、二つの比で、何種類もの図形が描けるのですけれども、、、、せっかくコンピュ－タでプログラムを作るのだし、　ギアをさらに何段も加えた図形を描けるようにしてみました。　ホンモノの歯車で作るのは　自分　不器用なのでムリですけれども、ＣＧなのでラクラクですー　
[/pukiwiki]

サンプル画像

コード自体は、かなりヘボいですが、それでも　結構　綺麗な模様が描けます。
SVGなど、ベクターデータを出力するようにすると楽しいかも。（宿題。　Blenderあたりを使うのが近道かな？）

スピログラフで描く図形は　サイクロイド曲線とかトロコイド曲線と呼ばれるもの。
スピログラフで描けるのは、内サイクロイドのみですが、両方描けるようにしました。

動作にはPILが必要です。
動作確認はCPython2.5で行ってるので、2.6以降だと、一部書き換えが必要かも？

# -*- coding: utf8 -*-
from __future__ import division
import Image,ImageDraw
from math import cos,sin,pi

W,H=2000,2000
im=Image.new("RGB",(W,H))

draw=ImageDraw.ImageDraw(im)
#rlist ギア比  マイナスで内トロコイド　プラスで外トロコイド
#scale  画像のスケーリング
#holelist　ギアのどのあたりにエンピツを差すか

#ヘッダの図形のデータ
rlist1=[1,-3,-5,7]; scale1=50; holelist1=[1,1,1,1]
rlist2=[3,-7,-9]; scale2=100; holelist2=[.9,.4,.3]

def nround(m,n):
    u"""線が閉じるには何周必要か計算する
      具体的には、mとnの最小公倍数でいいはずなのですが、、、
     スピログラフと同じような動作にしてあります """

    m,n=abs(m),abs(n)
    if m==n or m%n==0:
        return m//n
    if m<n :
        m,n=n,m
    b=0
    
    for k in xrange(1,100) :
        b+=n
        if b>0 :
          b-=m
        elif b==0:
            break 
    return k


def spiro(rlist,scale,holelist):
    draw.rectangle((0,0,W,H),fill=(255,255,255))

    ox,oy=W/2,H/2

    last=200
    k=13/19


    rtxt="_R_"+str(rlist)[1:-1]+"_hole_"+",".join(["%1.2f"%x for x in holelist])

    roundlist=[]
    ratio=1
    ratiolist=[ratio,]
    last=rlist[0]

    for r in rlist[1:]:
        rnd=nround(last,r)
        roundlist.append(rnd)
        ratio=ratio*last/r
        ratiolist.append(ratio)
        last=r

    rlist=map(lambda x: x*scale,rlist)

    imax0=300
    imax=imax0*reduce(lambda x,y:x*y,roundlist)

    maxdepth=len(rlist)

    def f(t,depth=0,x=ox,y=oy,maxdepth=maxdepth,ratio=1.0,lastr=rlist[0]):
        u"""この関数を、必要なだけ再帰呼び出しします　（遅いですが）
        スピログラフが何重にも入れ子になっている、と思ってください
        """
        if depth<maxdepth:
            r=rlist[depth]
            if not r : return x,y
            ratio=ratiolist[depth]
            hole=holelist[depth]
            t2=t*ratio
            r2=r*hole
            x2=x+r2*cos(t2)
            y2=y+r2*sin(t2)
            return f(t,depth+1,x2,y2,ratio=ratio,lastr=r)
        else:
            return x,y
        
    lastx,lasty=f(0)

    for i in xrange(1,imax+1):
        t=2.0*i*pi/imax0
        x,y=f(t)
        draw.line((lastx,lasty,x,y), fill=(0,0,0))
        lastx,lasty=x,y 
    im.save("spiro%s.jpg"%str(rtxt))

spiro(rlist1,scale1, holelist1)
spiro(rlist2,scale2, holelist2)