benchmarks_game/nbody.fz
# https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/nbody.html#nbody
nbody is
PI := f64 3.141592653589793
SOLAR_MASS := f64 4 * PI * PI
DAYS_PER_YEAR := f64 365.24
body(x, y, z, vx, vy, vz, mass f64) is
# change velocity
body(b body, vx, vy, vz f64) =>
body b.x b.y b.z vx vy vz b.mass
jupiter := body 4.84143144246472090 -1.16032004402742839 -1.03622044471123109E-1 (1.66007664274403694E-3 * DAYS_PER_YEAR) (7.69901118419740425E-3 * DAYS_PER_YEAR) (-6.90460016972063023E-5 * DAYS_PER_YEAR) (9.54791938424326609E-4 * SOLAR_MASS)
saturn := body 8.34336671824457987 4.12479856412430479 -4.03523417114321381E-1 (-2.76742510726862411E-3 * DAYS_PER_YEAR) (4.99852801234917238E-3 * DAYS_PER_YEAR) (2.30417297573763929E-5 * DAYS_PER_YEAR) (2.85885980666130812E-4 * SOLAR_MASS)
uranus := body 1.28943695621391310E+1 -1.51111514016986312E+1 -2.23307578892655734E-1 (2.96460137564761618E-3 * DAYS_PER_YEAR) (2.37847173959480950E-3 * DAYS_PER_YEAR) (-2.96589568540237556E-5 * DAYS_PER_YEAR) (4.36624404335156298E-5 * SOLAR_MASS)
neptune := body 1.53796971148509165E+1 -2.59193146099879641E+1 1.79258772950371181E-1 (2.68067772490389322E-3 * DAYS_PER_YEAR) (1.62824170038242295E-3 * DAYS_PER_YEAR) (-9.51592254519715870E-5 * DAYS_PER_YEAR) (5.15138902046611451E-5 * SOLAR_MASS)
sun := body 0 0 0 0 0 0 SOLAR_MASS
advance (bodies array body, dt f64) =>
for i:=0, i + 1
while i < bodies.count
for j:=i+1, j+1
while j < bodies.count
dx := bodies[i].x - bodies[j].x
dy := bodies[i].y - bodies[j].y
dz := bodies[i].z - bodies[j].z
dSquared := dx * dx + dy * dy + dz * dz
distance := f64.type.sqrt dSquared
mag := dt / (dSquared * distance)
vx := bodies[i].vx - dx * bodies[j].mass * mag
vy := bodies[i].vy - dy * bodies[j].mass * mag
vz := bodies[i].vz - dz * bodies[j].mass * mag
bodies.put(i, body bodies[i] vx vy vz)
vx := bodies[j].vx + dx * bodies[i].mass * mag
vy := bodies[j].vy + dy * bodies[i].mass * mag
vz := bodies[j].vz + dz * bodies[i].mass * mag
bodies.put(j, body bodies[i] vx vy vz)
bodies.map body (b -> body (b.x + dt * b.vx) (b.y + dt * b.vy) (b.z + dt * b.vz) b.vx b.vy b.vz b.mass)
energy (bodies array body) =>
e := mut 0.0
for i:=0, i+1
while i < bodies.count
bodyA := bodies[i]
e <- e.get + 0.5 * bodyA.mass * (bodyA.vx * bodyA.vx
+ bodyA.vy * bodyA.vy
+ bodyA.vz * bodyA.vz)
for j:=i+1, j+1
while j < bodies.count
bodyB := bodies[j]
dx := bodyA.x - bodyB.x
dy := bodyA.y - bodyB.y
dz := bodyA.z - bodyB.z
distance := f64.type.sqrt (dx*dx + dy*dy + dz*dz)
e <- e.get - (bodyA.mass * bodyB.mass) / distance
e.get
n := 1
start_constellation := [ sun, jupiter, saturn, uranus, neptune ]
px := start_constellation
.as_stream
.map f64 (body -> body.vx * body.mass)
.fold f64.type.sum
py := start_constellation
.as_stream
.map f64 (body -> body.vy * body.mass)
.fold f64.type.sum
pz := start_constellation
.as_stream
.map f64 (body -> body.vz * body.mass)
.fold f64.type.sum
offsetted_sun := body start_constellation[0] (-px / SOLAR_MASS) (-py / SOLAR_MASS) (-pz / SOLAR_MASS)
moving_bodies := mut [offsetted_sun, jupiter, saturn, uranus, neptune]
say (energy moving_bodies.get)
for i:=0, i+1
while i<n
moving_bodies <- advance moving_bodies.get 0.01
say (energy moving_bodies.get)
last changed: 2023-11-14