benchmarks_game/nbody_marray.fz

# https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/nbody.html#nbody
# version using Mutable_Array instead of array
nbody_marray is

  # NYI this is broken currently

  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 Mutable_Array body, dt f64) =>
    for i in bodies.indices do
      for j in (i+1..bodies.count-1) do
        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[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[j] := body bodies[j] vx vy vz

    for i in bodies.indices do
      b := bodies[i]
      bodies[i] := 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_of_body(b body) =>
    b.mass * (b.vx * b.vx
            + b.vy * b.vy
            + b.vz * b.vz)

  energy (bodies Mutable_Array body) =>
    (bodies.indices.map f64 (i ->
        bodyA := bodies[i]

        # NYI better name for var x. unit is kg²/m
        x := (bodies.as_list.drop(i+1).map f64 (bodyB ->
            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)
            bodyA.mass * bodyB.mass / distance
          )).fold(f64.type.sum)

        - x + 0.5 * energy_of_body bodyA
      )).fold(f64.type.sum)

  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)

  mi : mutate is
  mi.go ()->
    moving_bodies := (mutate.array body).type.new mi 5 offsetted_sun
    moving_bodies[0] := offsetted_sun
    moving_bodies[1]:= jupiter
    moving_bodies[2]:= saturn
    moving_bodies[3]:= uranus
    moving_bodies[4]:= neptune

    say (energy moving_bodies)

    for _ in (0..n-1) do
      advance moving_bodies 0.01

    say (energy moving_bodies)