# 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
Experimental Flying GameFly around in a plane. Some physics, but mainly just playing with websockets. If you can get a friend to play at the same time, you should be able to shoot each other down. |
Maths Exam GeneratorAuto-generated maths exams, with and without answers. Set at A-level / end of high school / beginning of university. |
How SpaceX land first stage boostersThe algorithms that SpaceX (probably) use to control their first stage boosters. Several animations. |
© Hugo2015. Session @sessionNumber