# 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
Atomistic Simulation of MetalsThis presents an interactive simulation of atoms making up a nanoscopic particle of metal. |
Mars ReentryA ThreeJS simulation of mars reentry in a spaceship. |
BodyWorks: Kinematics ToolA page with a javascript application where you can set body positions, calculate joint angles and animate human motion. |
© Hugo2015. Session @sessionNumber