Particles

A good chunk of my weekend went towards making a simple particle simulation. I like the amount of emergent behavior it demonstrates despite its very simple particle interaction rules.

Each particle is defined by its position, velocity, and color. Rather than use a library, I just wrote my own simple Vector2 module¹ to handle operations on 2d vectors (constructed as Vector2 { x :: Float, y :: Float }), and started using it for positions and velocities. We use Float instead of Double because it plays more nicely with OpenGL.

data Particle = Particle
  { pos :: Vector2
  , vel :: Vector2
  , color :: Color
  }
  deriving (Eq, Show)

Each particle color is one of a small set that we define ourselves:

data Color = Red | Green | Blue | Yellow
  deriving (Eq, Show, Enum, Bounded)

In order to draw a particle simulation, you first need to be able to draw a single particle. I found some example Haskell GLFW code that provided window-handling code and an OpenGL context.

In OpenGL, we can set a background color by setting the clear color, and then clearing the screen.

main :: IO ()
main = do

  ... -- setup

  glClearColor 0 0 0.1 1 -- dark blue
  glClear GL_COLOR_BUFFER_BIT

We’ll draw particles as little circles (actually, as little n-gons).

GL_TRIANGLE_FAN uses the first vertex as a central hub point, so we’ll place it at the center of our particle with glVertex2f p.pos.x p.pos.y. Then, triangles will fan out from there to points along the circumference of the circle.

For each vertex, we find its offset angle from the start, and convert into a position via cos and sin. Then, we place a vertex there. The triangle fan mode will fill in our shape for us as we go.

This is not the most efficient way to draw circles in OpenGL, but it worked well enough for the relatively small numbers of particles in this simulation. Processing time was overwhelmed by the application of particle force updates anyways.

drawParticle :: Particle -> IO ()
drawParticle p = do
  let vertices = 12

  glBegin GL_TRIANGLE_FAN
  setGlColor p.color
  glVertex2f p.pos.x p.pos.y

  forM_ [0..vertices] $ \i -> do
    let angle = i * 2 * pi / vertices
    glVertex2f
      (p.pos.x + particleRadius * cos angle)
      (p.pos.y + particleRadius * sin angle)

  glEnd

On each frame, we update all the particles, print out the frames per second² as of the last frame, draw all the particles, and then swap the newly drawn buffer into the existing window from GLFW.

draw :: Window -> [Particle] -> DiffTime -> IO [Particle]
draw window oldParticles deltaTime = do
  let
    particles = (flip map) oldParticles $ \particle ->
      updateParticle oldParticles particle deltaTime

  putStrLn $ "fps = " <> show (truncate $ 1 / deltaTime :: Int)

  glClear GL_COLOR_BUFFER_BIT
  forM_ particles drawParticle
  GLFW.swapBuffers window

  pure particles

The draw function takes a DiffTime representing how long it’s been since the last frame. We can thread time through³ to each draw call by getting the initial time, and then on each loop getting an updated time, sending the difference to draw as deltaTime, and having loop call itself with the new time.

Each call to loop reads window events, updates the time, draws particles given a TimeStep value, and then calls itself again to continue the loop.

main = do
  ... -- existing code
  time <- utctDayTime <$> getCurrentTime
  particles <- initialParticles
  loop particles time window

loop :: [Particle] -> DiffTime -> Window -> IO ()
loop oldParticles oldTime window = do
  GLFW.pollEvents
  time <- utctDayTime <$> getCurrentTime
  particles <- draw window oldParticles $ case timeStep of
    Fixed t -> t
    RealTime -> time - oldTime
  loop particles time window

Because this is a simulation, it’s nice to have the option of a fixed timestep. This helps avoid some amount of instability in the simulation, at the cost of no longer being able to update visually in real time.

The RealTime option works well on fast simulations where visuals matter more than accuracy, but if the simulation seems erratic it’s better to switch to Fixed.

data TimeStep = Fixed DiffTime | RealTime

We can initialize a simulation with a collection of random particles.

initialParticles :: IO [Particle]
initialParticles =
  sequence $ replicate numParticles randomParticle
  where
    randomParticle = do
      randomX :: Float <- randomRIO (-1, 1)
      randomY :: Float <- randomRIO (-1, 1)
      randomColor :: Color <- toEnum <$> randomRIO
        ( fromEnum (minBound :: Color)
        , fromEnum (maxBound :: Color)
        )
      pure $ Particle
        { pos = Vector2 randomX randomY
        , vel = zeroVec
        , color = randomColor
        }

Here the bounds are set between (-1, 1) to correspond with OpenGL coordinates. We can use these bounds as collision boundaries by clamping particle positions to stay within them, and inverting velocities upon collision for a “bounce” effect. Including a simple damping factor makes it seem like the particles are losing energy when colliding with the walls.

resolveWallCollisions :: Particle -> Particle
resolveWallCollisions p = p
  { pos = p.pos
    { x = clamp p.pos.x
    , y = clamp p.pos.y
    }
  , vel = p.vel
    { x = adjustVelocity p.pos.x p.vel.x
    , y = adjustVelocity p.pos.y p.vel.y
    }
  }
  where
    upperBound = 1
    lowerBound = negate upperBound

    realLowerBound = lowerBound + particleRadius
    realUpperBound = upperBound - particleRadius

    adjustVelocity position velocity =
      if isCollision position
      then -velocity * collisionDamping
      else velocity

    isCollision n =
      n < realLowerBound ||
      n > realUpperBound

    clamp n =
      if
        | n < realLowerBound -> realLowerBound
        | n > realUpperBound -> realUpperBound
        | otherwise -> n

The final piece of our simulation is the inclusion of attraction and repulsion forces between particles. There are two radii associated with every particle: a smaller one, within which other particles are repulsed, and a big one within which other particles are attracted.

I spent a while messing with attraction and repulsion radii, and implementing these forces as inverse square laws. Then I found the video The code behind Particle Life, and decided to use its math instead.

We find the accelerative force between each pair of particles, and then sum up those influences to get the overall force acting on a single particle. Once we do this for all particles, we’ll have calculated how each particle will move thanks to every other.

Unfortunately, this algorithm is in \(O(n^2)\). Faster methods exist! For example, spatial partitioning—instead of looking at every particle, only look at those in nearby partitions of the overall space. I implemented it in a separate space-partitioning branch. However, if you just want your simulation to run, the naive algorithm should be fine.⁴

attractionAcceleration :: [Particle] -> Particle -> Vector2
attractionAcceleration particles p =
  sumVec $ (flip map) otherParticles $ \other ->
    let dist = distance p.pos other.pos
        dir = scale (other.pos `minus` p.pos) (1 / dist)
    in
      if
        | dist < repulsionRadius ->
            scale dir $ dist / repulsionRadius - 1
        | dist < attractionRadius ->
            scale dir $ attractionFactor p.color other.color *(1 - (abs $ 2 * dist - 1 - repulsionRadius) / (1 - repulsionRadius))
        | otherwise -> zeroVec

  where
    attractionRadius = 1
    otherParticles = filter (/=p) particles

We can associate each pair of colors with an attractionFactor between them. For example, here red particles are attracted to other red particles with a strength of 0.6.

An efficient way to represent these pairings is with a triangular matrix. I considered using a Bimap for this, but instead went with a solution that lets us still use GHC’s totality checking on Haskell enums. This is nice because when we add a new Color, the compiler will tell us which pairs of colors we’re missing here.

attractionFactor :: Color -> Color -> Float
attractionFactor a b = case (a, b) of
  (Red   , Red)    ->  0.6
  (Red   , Green)  -> -0.3
  (Red   , Blue)   ->  0.2
  (Red   , Yellow) -> -0.4
  (Green , Green)  ->  2.0
  (Green , Blue)   ->  0.3
  (Green , Yellow) -> -0.2
  (Blue  , Blue)   -> -0.4
  (Blue  , Yellow) ->  0.5
  (Yellow, Yellow) -> -0.2
  -- repeats of above, but want to enumerate
  -- all cases for totality checking
  (Green , Red)   -> swap
  (Blue  , Red)   -> swap
  (Blue  , Green) -> swap
  (Yellow, Red)   -> swap
  (Yellow, Green) -> swap
  (Yellow, Blue)  -> swap
  where
    swap = attractionFactor b a

The final update of each particle involves updating the velocity by applying gravity (a simple downwards acceleration), applying the attractive/repulsive forces between this particle and every other particle, scaling down by a slowdownConstant factor, updating the position from this new velocity, and finally checking for and resolving collisions with the walls.

As an aside, it might be interesting to do something like Verlet integration to implement collisions between particles as well, but I didn’t do that here.

updateParticle :: [Particle] -> Particle -> DiffTime -> Particle
updateParticle allParticles p deltaTime =
  resolveWallCollisions $ p
    { pos = p.pos `plus` (newVel `scale` realToFrac deltaTime)
    , vel = newVel
    }
  where
    newVel =
      sumVec
        [ p.vel
        , scale downVec $ gravity * realToFrac deltaTime -- gravity
        , attractionAcceleration allParticles p `scale` realToFrac deltaTime -- particle attractions
        ]
        `scale` slowdownConstant

There are a large variety of controllable parameters that result in different-looking simulations. One effect I enjoyed was accomplished by turning off the screen clearing between each frame. This led to particles being drawn as trails of points over time rather than single instantaneous points.

Again, here’s the code for this project.