# Generic Include Statements

```
component particlesim
import List.{...}
export Executable
```

# Objects and Overloaded Operators

```
object Vect3D(var x : RR64, var y : RR64, var z : RR64)
end
object Particle(var pos:Vect3D,var vel:Vect3D,var radius:RR64,var mass:RR64)
end
opr juxtaposition(v:Vect3D, n:Number): Vect3D = Vect3D(v.x n, v.y n, v.z n)
opr +(a:Vect3D,b:Vect3D): Vect3D = Vect3D(a.x + b.x, a.y + b.y, a.z + b.z)
opr DOT(a:Vect3D,b:Vect3D):RR64 = a.x b.x + a.y b.y + a.z b.z
magnitude(a:Vect3D):RR64 = SQRT (a.x a.x + a.y a.y + a.z a.z)
```

# Collision Code

```
collisions(pop,size,timestep) = do
var i : ZZ32 = 0
while i < size do
for j <- 0:(size-1) do
var distVect: Vect3D = pop[j].pos + (pop[i].pos (-1.0))
mag = magnitude(distVect)
if i NE j AND mag < (pop[i].radius + pop[j].radius) then
atomic do
vcm = ((pop[i].vel pop[i].mass) + (pop[j].vel pop[j].mass)) (1 / (pop[i].mass + pop[j].mass))
pop[i].vel := pop[i].vel + (vcm (-1.0))
pop[j].vel := pop[j].vel + (vcm (-1.0))
distVect := distVect (1 / mag)
if (distVect DOT pop[i].vel GE 0) then
var vnorm: Vect3D = distVect (distVect DOT pop[i].vel)
var e : RR64 := magnitude(vnorm) 0.34
e := e^(-0.234)
if (e GT 1) then
e := 1
end
pop[i].vel := pop[i].vel + (vnorm ((-1.0) (1+e)))
vnorm := distVect (distVect DOT pop[j].vel)
pop[j].vel := pop[j].vel + (vnorm ((-1.0) (1+e)))
end
pop[i].vel := pop[i].vel + vcm
pop[j].vel := pop[j].vel + vcm
end
end
end
i += 1
end
end
```

# Gravity

This code was a little trickier to write. Because of the implicit parallelism of for loops, the outer iteration is implemented using a sequential while loop in order to ensure thread-safe performance. There isn't much else of note in this code except for the clear usage of adjacency for multiplication and extensive use of overloaded operators.

```
gravity(pop,size,timestep) = do
var i : ZZ32 = 0
while i < size do
for j <- (i):(size-1) do
if i NE j then
var distVect: Vect3D = pop[i].pos + (pop[j].pos (-1.0))
dist = magnitude(distVect)
softDist = dist + 0.0000001
mag = timestep / (softDist softDist dist)
pop[i].vel := pop[i].vel + (distVect (-1.0 pop[j].mass) mag)
pop[j].vel := pop[j].vel + (distVect pop[i].mass mag)
end
end
i += 1
end
end
```

# Stepper

This code is called every timestep in order to actually perform the simulation given the particle system.

```
step(pop,size,timestep) = do
gravity(pop,size,timestep)
collisions(pop,size,timestep)
for i <- 0:(size-1) do
var stepVector : Vect3D = pop[i].vel timestep
pop[i].pos := pop[i].pos + stepVector
end
end
```

The step function takes an array of particles that is referred to as a 'population', the total size of said array, and a timestep. It first performs all of the gravity forcing betweeen the particles, then handles all of the collisions for the current time step, and finally does the integrator 'drift' step. Note that the for structure is parallel everywhere in fortress, so the step of the simulation happens concurrently.

# Main (a la C)

```
run() = do
println "Starting particle simulation..."
var i : ZZ32 = 0
stepsize = 0.0001
popsize = 2
steps = 1000000
pop : Particle[2] = [
Particle(Vect3D(1.0,0.0,0.0) ,Vect3D(0.0,1.0,0.0),0.1,0.000001)
Particle(Vect3D(0.0,0.0,0.0),Vect3D(0.0,0.0,0.0),0.5,1.0)]
while i < steps do
var j : ZZ32 = 0
step(pop,popsize,stepsize)
while j < popsize do
println i " " j " " pop[j].pos.x " " pop[j].pos.y " " pop[j].pos.z " " pop[j].vel.x " " pop[j].vel.y " " pop[j].vel.z " " pop[j].radius
j += 1
end
i += 1
end
end
end particlesim
```