Welcome, guest | Sign In | My Account | Store | Cart

Hofstadter Butterfly Fractal

Python, 77 lines
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
# Hofstadter Butterfly Fractal
# http://en.wikipedia.org/wiki/Hofstadter%27s_butterfly
# Wolfgang Kinzel/Georg Reents,"Physics by Computer" Springer Press (1998)
# FB36 - 20130922
import math
from PIL import Image
imgSize = 200
image = Image.new("RGB", (imgSize, imgSize))
pixels = image.load() 

def gcd(a, b): # Greatest Common Divisor
    if b == 0: return a
    return gcd(b, a % b)

pi2 = math.pi * 2.0
MAXX = imgSize + 1
MAXY = imgSize + 1
qmax = imgSize
for q in range(4, qmax, 2):
    print str(100 * q / qmax).zfill(2) + "%"
    for p in range(1, q, 2):
        if gcd(p, q) <= 1:
            sigma = pi2 * p / q
            nold = 0
            ie = 0
            for ie in range(0, MAXY + 2):
                e = 8.0 * ie / MAXY - 4.0 - 4.0 / MAXY
                n = 0
                polyold = 1.0
                poly = 2.0 * math.cos(sigma) - e
                if polyold * poly < 0.0: n += 1

                for m in range(2, q / 2):
                    polynew = (2.0 * math.cos(sigma * m) - e) * poly - polyold
                    if poly * polynew < 0.0: n += 1
                    polyold = poly
                    poly = polynew

                polyold = 1.0
                poly = 2.0 - e
                if polyold * poly < 0.0: n += 1
                polynew = (2.0 * math.cos(sigma) - e) * poly - 2.0 * polyold
                if poly * polynew < 0.0: n += 1
                polyold = poly
                poly = polynew

                for m in range(2, q / 2):
                    polynew = (2.0 * math.cos(sigma * m) - e) * poly - polyold
                    if poly * polynew < 0.0: n += 1
                    polyold = poly
                    poly = polynew

                polynew = (2.0 * math.cos(sigma * q / 2.0) - e) * poly - 2.0 * polyold
                if poly * polynew < 0.0: n += 1

                polyold = 1.0
                poly = 2.0 - e
                if polyold * poly < 0.0: n += 1
                polynew = (2.0 * math.cos(sigma) - e) * poly - 2.0 * polyold
                if poly * polynew < 0.0: n += 1
                polyold = poly
                poly = polynew

                for m in range(2, q / 2):
                    polynew = (2.0 * math.cos(sigma * m) - e) * poly - polyold
                    if poly * polynew < 0.0: n += 1
                    polyold = poly
                    poly = polynew

                polynew = (2.0 * math.cos(sigma * q / 2.0) - e) * poly - 2.0 * polyold
                if poly * polynew < 0.0: n += 1
                if n > nold:
                    pixels[int(MAXY - ie), int(MAXX * p / q)] = (255, 255, 255)
                    pixels[int(MAXX * p / q), int(MAXY - ie)] = (255, 255, 255)
                nold = n

image.save("HofstadterButterflyFractal.png", "PNG")
Created by FB36 on Sun, 22 Sep 2013 (MIT)
Python recipes (4591)
FB36's recipes (148)

Required Modules

Other Information and Tasks