Welcome, guest | Sign In | My Account | Store | Cart
```# Hofstadter Butterfly Fractal
# 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))

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

```

#### Diff to Previous Revision

```--- revision 1 2013-09-22 22:25:07
+++ revision 2 2013-09-23 03:18:02
@@ -12,11 +12,9 @@
if b == 0: return a
return gcd(b, a % b)

-pix = [[0 for x in range(imgSize)] for y in range(imgSize)]
+pi2 = math.pi * 2.0
MAXX = imgSize + 1
MAXY = imgSize + 1
-pi2 = math.pi * 2.0
-nold = 0
qmax = imgSize
for q in range(4, qmax, 2):
print str(100 * q / qmax).zfill(2) + "%"
@@ -72,12 +70,8 @@
polynew = (2.0 * math.cos(sigma * q / 2.0) - e) * poly - 2.0 * polyold
if poly * polynew < 0.0: n += 1
if n > nold:
-                    pix[int(MAXX * p / q)][int(MAXY - ie)] = 1
-                    pix[int(MAXY - ie)][int(MAXX * p / q)] = 1
+                    pixels[int(MAXY - ie), int(MAXX * p / q)] = (255, 255, 255)
+                    pixels[int(MAXX * p / q), int(MAXY - ie)] = (255, 255, 255)
nold = n

-# paint
-for ky in range(imgSize):
-    for kx in range(imgSize):
-        pixels[kx, ky] = (pix[ky][kx] * 255, pix[ky][kx] * 255, pix[ky][kx] * 255)