|
| 1 | +// N body problem benchmark - Umka version |
| 2 | + |
| 3 | +import "../import/std.um" |
| 4 | + |
| 5 | +type Body = struct { |
| 6 | + x, y, z, vx, vy, vz, mass: real |
| 7 | +} |
| 8 | + |
| 9 | +const ( |
| 10 | + pi = 3.14159265358979323846 |
| 11 | + solarMass = 4 * pi * pi |
| 12 | + daysPerYear = 365.24 |
| 13 | +) |
| 14 | + |
| 15 | +fn offsetMomentum(b: ^Body, px, py, pz: real) { |
| 16 | + b.vx = -px / solarMass |
| 17 | + b.vy = -py / solarMass |
| 18 | + b.vz = -pz / solarMass |
| 19 | +} |
| 20 | + |
| 21 | +type System = []^Body |
| 22 | + |
| 23 | +fn NewSystem(body: ^[]Body, len: int, sys: ^System) { |
| 24 | + for i := 0; i < len; i++ { |
| 25 | + sys[i] = &body[i] |
| 26 | + } |
| 27 | + var px, py, pz: real |
| 28 | + for i := 0; i < len; i++ { |
| 29 | + body := sys[i] |
| 30 | + px += body.vx * body.mass |
| 31 | + py += body.vy * body.mass |
| 32 | + pz += body.vz * body.mass |
| 33 | + } |
| 34 | + offsetMomentum(sys[0], px, py, pz) |
| 35 | +} |
| 36 | + |
| 37 | +fn energy(sys: ^System, len: int): real { |
| 38 | + var e: real |
| 39 | + for i := 0; i < len; i++ { |
| 40 | + body := sys[i] |
| 41 | + e += 0.5 * body.mass * (body.vx*body.vx + body.vy*body.vy + body.vz*body.vz) |
| 42 | + for j := i + 1; j < len; j++ { |
| 43 | + body2 := sys[j] |
| 44 | + dx := body.x - body2.x |
| 45 | + dy := body.y - body2.y |
| 46 | + dz := body.z - body2.z |
| 47 | + distance := sqrt(dx*dx + dy*dy + dz*dz) |
| 48 | + e -= (body.mass * body2.mass) / distance |
| 49 | + } |
| 50 | + } |
| 51 | + return e |
| 52 | +} |
| 53 | + |
| 54 | +fn advance(sys: ^System, len: int, dt: real) { |
| 55 | + for i := 0; i < len; i++ { |
| 56 | + body := sys[i] |
| 57 | + for j := i + 1; j < len; j++ { |
| 58 | + body2 := sys[j] |
| 59 | + dx := body.x - body2.x |
| 60 | + dy := body.y - body2.y |
| 61 | + dz := body.z - body2.z |
| 62 | + |
| 63 | + dSquared := dx*dx + dy*dy + dz*dz |
| 64 | + distance := sqrt(dSquared) |
| 65 | + mag := dt / (dSquared * distance) |
| 66 | + |
| 67 | + body.vx -= dx * body2.mass * mag |
| 68 | + body.vy -= dy * body2.mass * mag |
| 69 | + body.vz -= dz * body2.mass * mag |
| 70 | + |
| 71 | + body2.vx += dx * body.mass * mag |
| 72 | + body2.vy += dy * body.mass * mag |
| 73 | + body2.vz += dz * body.mass * mag |
| 74 | + } |
| 75 | + } |
| 76 | + |
| 77 | + for i := 0; i < len; i++ { |
| 78 | + body := sys[i] |
| 79 | + body.x += dt * body.vx |
| 80 | + body.y += dt * body.vy |
| 81 | + body.z += dt * body.vz |
| 82 | + } |
| 83 | +} |
| 84 | + |
| 85 | + jupiter := Body{ |
| 86 | + x: 4.84143144246472090e+00, |
| 87 | + y: -1.16032004402742839e+00, |
| 88 | + z: -1.03622044471123109e-01, |
| 89 | + vx: 1.66007664274403694e-03 * daysPerYear, |
| 90 | + vy: 7.69901118419740425e-03 * daysPerYear, |
| 91 | + vz: -6.90460016972063023e-05 * daysPerYear, |
| 92 | + mass: 9.54791938424326609e-04 * solarMass} |
| 93 | + |
| 94 | + saturn := Body{ |
| 95 | + x: 8.34336671824457987e+00, |
| 96 | + y: 4.12479856412430479e+00, |
| 97 | + z: -4.03523417114321381e-01, |
| 98 | + vx: -2.76742510726862411e-03 * daysPerYear, |
| 99 | + vy: 4.99852801234917238e-03 * daysPerYear, |
| 100 | + vz: 2.30417297573763929e-05 * daysPerYear, |
| 101 | + mass: 2.85885980666130812e-04 * solarMass} |
| 102 | + |
| 103 | + uranus := Body{ |
| 104 | + x: 1.28943695621391310e+01, |
| 105 | + y: -1.51111514016986312e+01, |
| 106 | + z: -2.23307578892655734e-01, |
| 107 | + vx: 2.96460137564761618e-03 * daysPerYear, |
| 108 | + vy: 2.37847173959480950e-03 * daysPerYear, |
| 109 | + vz: -2.96589568540237556e-05 * daysPerYear, |
| 110 | + mass: 4.36624404335156298e-05 * solarMass} |
| 111 | + |
| 112 | + neptune := Body{ |
| 113 | + x: 1.53796971148509165e+01, |
| 114 | + y: -2.59193146099879641e+01, |
| 115 | + z: 1.79258772950371181e-01, |
| 116 | + vx: 2.68067772490389322e-03 * daysPerYear, |
| 117 | + vy: 1.62824170038242295e-03 * daysPerYear, |
| 118 | + vz: -9.51592254519715870e-05 * daysPerYear, |
| 119 | + mass: 5.15138902046611451e-05 * solarMass} |
| 120 | + |
| 121 | + sun := Body{ |
| 122 | + x: 0, y: 0, z: 0, vx: 0, vy: 0, vz: 0, |
| 123 | + mass: solarMass} |
| 124 | + |
| 125 | + |
| 126 | +fn main() { |
| 127 | + if std.argc() == 3 { |
| 128 | + iter := std.atoi(std.argv(2)) |
| 129 | + |
| 130 | + const num = 5 |
| 131 | + bodies := [num]Body{sun, jupiter, saturn, uranus, neptune} |
| 132 | + var system: [num]^Body |
| 133 | + NewSystem(&bodies, num, &system) |
| 134 | + printf("{.9f}\n", energy(&system, num)) |
| 135 | + for i := 0; i < iter; i++ { |
| 136 | + advance(&system, num, 0.01) |
| 137 | + } |
| 138 | + printf("{.9f}\n", energy(&system, num)) |
| 139 | + } |
| 140 | +} |
| 141 | + |
0 commit comments