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

 ``` 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") ```