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

Gravner-Griffeath snowflake simulation.

For better results increase the image size and number of growth steps and prepare to wait!

Python, 170 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 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170``` ```# Gravner-Griffeath Snowflake Simulation # http://psoup.math.wisc.edu/Snowfakes.htm # "A typical simulated crystal reaches a final diameter of 400-600 cells # over 10,000-100,000 updates (growth steps)..." # FB36 - 20130526 import math; import random from PIL import Image imgx = 300; imgy = 300 # image size imgx1 = imgx - 1; imgy1 = imgy - 1 image = Image.new("RGB", (imgx, imgy)) pixels = image.load() ver = "deterministic" # random.choice(["deterministic", "randomized"]) print "Version: " + ver maxIt = 3000 # 10000 - 100000 # of growth steps p = 0.58 # random.random() * 0.6 + 0.3 # rho: homogeneous vapor density k = random.random() * 0.05 # kappa b = 2.0 # random.random() * 1.95 + 1.05 # beta a = random.random() * 0.3 # alpha t = random.random() * 0.5595 + 0.02 # theta m = random.random() * 0.01 # mu g = 0.0000515 # gamma s = random.random() * -0.5 # sigma print "Parameters:" print "rho = ", p print "kappa = ", k print "beta = ", b print "alpha = ", a print "theta = ", t print "mu = ", m print "gamma = ", g print "sigma = ", s print mx = imgx; my = imgy # width and height of 2DCA dx = [-1, 0, -1, 1, 0, 1]; dy = [-1, -1, 0, 0, 1, 1] # 6 directions to grow # set initial state # ice cells belong to crystal at = [[[0 for x in range(mx)] for y in range(my)] for z in range(2)] # boundary masses of cells bt = [[[0.0 for x in range(mx)] for y in range(my)] for z in range(2)] # crystal masses of cells ct = [[[0.0 for x in range(mx)] for y in range(my)] for z in range(2)] # diffusive masses of cells dt = [[[p for x in range(mx)] for y in range(my)] for z in range(2)] # set ice seed ox = (mx - 1) / 2; oy = (my - 1) / 2 at[0][oy][ox] = 1 ct[0][oy][ox] = 1.0 dt[0][oy][ox] = 0.0 def isBoundary(x, y): global dx, dy, at, za flag = False if at[za][y][x] == 0: for j in range(6): # neighbors jx = x + dx[j]; jy = y + dy[j] if jx >= 0 and jx < mx and jy >= 0 and jy < my: if at[za][jy][jx] == 1: flag = True; break return flag def numIce(x, y): # of ice neighbors of boundary cell global dx, dy, at, za ni = 0 if at[za][y][x] == 0: for j in range(6): # neighbors jx = x + dx[j]; jy = y + dy[j] if jx >= 0 and jx < mx and jy >= 0 and jy < my: if at[za][jy][jx] == 1: ni += 1 return ni # total diffusive mass def difMass(x, y): global dx, dy, dt, zd wsum = dt[zd][y][x] for j in range(6): # neighbors jx = x + dx[j]; jy = y + dy[j] if jx >= 0 and jx < mx and jy >= 0 and jy < my: wsum += dt[zd][jy][jx] return wsum za = 0; wa = 1 zb = 0; wb = 1 zc = 0; wc = 1 zd = 0; wd = 1 for i in range(maxIt): # growth steps print "Growth Step: " + str(i + 1) + " of " + str(maxIt) # step 1: diffusion for iy in range(my): for ix in range(mx): if not(at[za][iy][ix] == 1 or isBoundary(ix, iy)): dt[wd][iy][ix] = difMass(ix, iy) / 7.0 elif isBoundary(ix, iy): wsum = dt[zd][iy][ix] for j in range(6): # neighbors jx = ix + dx[j]; jy = iy + dy[j] if jx >= 0 and jx < mx and jy >= 0 and jy < my: if at[za][jy][jx] == 1: wsum += dt[zd][iy][ix] else: wsum += dt[zd][jy][jx] dt[wd][iy][ix] = wsum / 7.0 zd = 1 - zd; wd = 1 - wd # switch planes # step 2: freezing for iy in range(my): for ix in range(mx): if isBoundary(ix, iy): bt[wb][iy][ix] = bt[zb][iy][ix] + (1.0 - k) * dt[zd][iy][ix] ct[wc][iy][ix] = ct[zc][iy][ix] + k * dt[zd][iy][ix] dt[wd][iy][ix] = 0.0 zb = 1 - zb; wb = 1 - wb # switch planes zc = 1 - zc; wc = 1 - wc # switch planes zd = 1 - zd; wd = 1 - wd # switch planes # step 3: attachment for iy in range(my): for ix in range(mx): nIce = numIce(ix, iy) if nIce > 0: if (nIce == 1 or nIce == 2) and bt[zb][iy][ix] >= b: at[wa][iy][ix] = 1 if nIce == 3: if bt[zb][iy][ix] >= 1.0: at[wa][iy][ix] = 1 elif bt[zb][iy][ix] >= a: if difMass(ix, iy) < t: at[wa][iy][ix] = 1 if nIce >= 4: at[wa][iy][ix] = 1 if at[wa][iy][ix] == 1: ct[wc][iy][ix] = bt[zb][iy][ix] + ct[zc][iy][ix]; bt[zb][iy][ix] = 0.0 za = 1 - za; wa = 1 - wa # switch planes zb = 1 - zb; wb = 1 - wb # switch planes zc = 1 - zc; wc = 1 - wc # switch planes # step 4: melting for iy in range(my): for ix in range(mx): if isBoundary(ix, iy): bt[wb][iy][ix] = bt[zb][iy][ix] * (1.0 - m) ct[wc][iy][ix] = ct[zc][iy][ix] * (1.0 - g) dt[wd][iy][ix] = dt[zd][iy][ix] + bt[zb][iy][ix] * m + ct[zc][iy][ix] * g zb = 1 - zb; wb = 1 - wb # switch planes zc = 1 - zc; wc = 1 - wc # switch planes zd = 1 - zd; wd = 1 - wd # switch planes # step 5: noise if ver == "randomized": for iy in range(my): for ix in range(mx): if at[wa][iy][ix] == 0: # if not ice dt[wd][iy][ix] = dt[zd][iy][ix] * (1.0 + s * random.choice([1, -1])) zd = 1 - zd; wd = 1 - wd # switch planes # paint final state of the snowflake an45 = - math.pi / 4.0 sn45 = math.sin(an45); cs45 = math.cos(an45) scale = math.sqrt(3.0); ox = imgx1 / 2.0; oy = imgy1 / 2.0 for ky in range(imgy): for kx in range(imgx): # apply geometric transformation (scaling and rotation) tx = kx - ox; ty = (ky - oy) * scale tx0 = tx * cs45 - ty * sn45 + ox ty = tx * sn45 + ty * cs45 + oy; tx = tx0 if tx >= 0 and tx <= imgx1 and ty >= 0 and ty <= imgy1: c = at[wa][int((my - 1) * ty / imgy1)][int((mx - 1) * tx / imgx1)] pixels[kx, ky] = (c * 255, c * 255, c * 255) image.save("Gravner-Griffeath_Snowfake_Simulation.png", "PNG") ```
 Created by FB36 on Sun, 26 May 2013 (MIT)