# Keep track of what is ice and what's not: A : 2d Array from -alot to +alot A[0,0] = 1 while True: place a particle at (x:Int,y:Int) somewhere random but far away. while True: move the particle a small distance randomly, updating x and y if A[x,y]: A[previous x, previous y] = 1 break
This gives us a fractal shape like this:
This has some of the branching that real ice crystals in snowflakes have, but it is far more random than real life. There are lots of reasons for this, but one part of it is that simulating individual particles is just too random. In real life, these crystals may have millions of trillions of water molecules in them, which is more than we could simulate individually.
So how do we tackle this? We can model the particle density as a 2- or 3-Dimensional array, and model diffusion on that array.
Diffusion can be modelled with a convolution between the particle density array and a Green's function - a gaussian kernel. Basically, if you know the particle density, ρ(x,y,t), at time t, then the particle density at time t + 𝛿t is given by the following:
ρ(x,y,t + 𝛿t) = ρ(x,y,t) ∗ G,
where G(x,y) = exp (-(x^2 + y^2) / s^2), and ∗ is the convolution operator.
To simulate the rest of the process, we also need to set the density to zero anywhere where there is an ice crystal --- to model the idea that any particles hitting the crystal will stick. We also need to allow the crystal to grow into adjacent cells according to adjacent particle density.
There is a fair bit more to explain here, but It's probably best for the interested reader to get a jupyter notebook up and running, and just copy the following code in and play around with it:
Note that the convolutions are done using Fourier transforms, which is efficient, and benefits from being able to take one of the Fourier transforms out of the inner loop.
import numpy as np import math from PIL import Image from tqdm import tqdm n = 800 # This def getSmoother(n): smoother = np.zeros([n, n]) smoother[0, 0] = 1 smoother[0, 1] = 1 smoother[1, 0] = 1 smoother[-1, 0] = 1 smoother[0, -1] = 1 smoother[1, -1] = 1 smoother[-1, 1] = 1 return smoother def getG(n): G = np.zeros([n, n]) for i in range(n): for j in range(n): i2 = i j2 = j if i2 > n/2: i2 -= n if j2 > n/2: j2 -= n # we're on a hexagonal grid, so this # is the x displacement for two tiles dx = i2 + 0.5 * j2 # and the y displacement dy = j2 * math.sqrt(3)/2 r = math.sqrt(dx*dx+dy*dy) G[i, j] = math.exp(-r*r/8) G = G / np.sum(np.sum(G)) return G frame = 0 def save(A): global frame B = A * 0.2 n = len(A) B[B < 0.01] = 0 # This handles the hexagonal grid, somewhat approximately: B = np.array( [B[i, range((n // 2 - i // 2), (n - i // 2), 1)] for i in range(n // 4, (3 * n) // 4, 1) if (i % 5) != 0]) def red(x): return np.uint(math.exp(-x * 4.0) * 255) def green(x): return np.uint(math.exp(-x) * 255) def blue(x): return np.uint(math.exp(-x / 4.0) * 255) col = [[[red(x), green(x), blue(x), 255] for x in y] for y in B] im = Image.fromarray(np.uint8(col)) im.save(f"frames/im{10000 + frame}.png") frame += 1 def main(): n = 1200 A = np.zeros([n, n]) A[n//2, n//2] = 1 G = getG(n) smoother = getSmoother(n) rho = np.ones([n, n]) sf = np.fft.fft2(smoother) gf = np.fft.fft2(G) for i in tqdm(range(5000)): B = A.copy() B[B < 1] = 0 B[B > 1] = 1 A2 = np.fft.ifft2(np.multiply(sf, np.fft.fft2(B))).real A2[A2 > 0.9] = 1 A2[A2 < 0.9] = 0 A += A2 * (rho * 0.5) A[A > 5] = 5 rho[A > 0.5] = 0 rho = np.fft.ifft2(np.multiply(gf, np.fft.fft2(rho))).real rho = rho / np.mean(np.mean(rho)) + np.random.randn(n, n) * 0.02 if (i%5)==0: save(A) save(A) if __name__ == '__main__': main()
If anyone wants a simple python code that demonstrates using Fourier transforms to grow an ice crystal, here it is.
Enjoy
AtmosphereA simple demo of a simulation of an atmosphere. It looks quite cool, but there's not a lot you can do with it yet, and the physics isn't yet all that accurate. |
How to simulate fractal ice crystal growth in PythonThis presents python code to draw snowflakes, simulating a diffusion process with Fourier transforms. |
The double slit and observersA look at the double slit experiment. The first half is meant to be a clear explanation, using simulations. The second half discusses some of the philosophy / interpretations of quantum mecahnics. |
© Hugo2015. Session @sessionNumber