Gravner-Griffeath snowflake simulation.
For better results increase the image size and number of growth steps and prepare to wait!
Python, 170 lines
Download
Copy to clipboard1 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")
|
