|
| 1 | +package main |
| 2 | + |
| 3 | +import "core:fmt" |
| 4 | +import "core:math" |
| 5 | +import "core:runtime" |
| 6 | +import "core:strconv" |
| 7 | + |
| 8 | +PI :: 3.141592653589793 |
| 9 | +SOLAR_MASS :: 4 * PI * PI |
| 10 | +DAYS_PER_YEAR :: 365.24 |
| 11 | + |
| 12 | +Vec3 :: [3]f64 |
| 13 | +Body :: struct { pos: Vec3, velocity: Vec3, mass: f64 } |
| 14 | +Sys :: [5]Body |
| 15 | + |
| 16 | +init_body :: proc(x, y, z, vx, vy, vz, mass : f64) -> Body { |
| 17 | + return Body { [3]f64 {x, y, z}, [3]f64 {vx, vy, vz}*DAYS_PER_YEAR, mass*SOLAR_MASS } |
| 18 | +} |
| 19 | + |
| 20 | +init_sys :: proc() -> Sys { |
| 21 | + return Sys { |
| 22 | + // Sun |
| 23 | + init_body(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0), |
| 24 | + // Jupiter |
| 25 | + init_body(4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01, |
| 26 | + 1.66007664274403694e-03, 7.69901118419740425e-03, -6.90460016972063023e-05, |
| 27 | + 9.54791938424326609e-04), |
| 28 | + // Saturn |
| 29 | + init_body(8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01, |
| 30 | + -2.76742510726862411e-03, 4.99852801234917238e-03, 2.30417297573763929e-05, |
| 31 | + 2.85885980666130812e-04), |
| 32 | + // Uranus |
| 33 | + init_body(1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01, |
| 34 | + 2.96460137564761618e-03, 2.37847173959480950e-03, -2.96589568540237556e-05, |
| 35 | + 4.36624404335156298e-05), |
| 36 | + // Neptune |
| 37 | + init_body(1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01, |
| 38 | + 2.68067772490389322e-03, 1.62824170038242295e-03, -9.51592254519715870e-05, |
| 39 | + 5.15138902046611451e-05), |
| 40 | + } |
| 41 | +} |
| 42 | + |
| 43 | +advance :: proc(sys: ^Sys, dt: f64) { |
| 44 | + for b, i in sys { |
| 45 | + v := b.velocity |
| 46 | + for j in (i+1)..<5 { |
| 47 | + b2 := &sys[j] |
| 48 | + dpos := b.pos - b2.pos |
| 49 | + distance_square := sum(dpos * dpos) |
| 50 | + distance := math.sqrt_f64(distance_square) |
| 51 | + mag := dt / (distance * distance_square) |
| 52 | + v -= dpos * (mag * b2.mass) |
| 53 | + b2.velocity += dpos * (mag * b.mass) |
| 54 | + } |
| 55 | + b.velocity = v |
| 56 | + b.pos += v * dt |
| 57 | + } |
| 58 | +} |
| 59 | + |
| 60 | +offset_momentum :: proc(sys: ^Sys) { |
| 61 | + p := Vec3 {} |
| 62 | + for b, _ in sys { |
| 63 | + p -= b.velocity * b.mass |
| 64 | + } |
| 65 | + sys[0].velocity = p / SOLAR_MASS |
| 66 | +} |
| 67 | + |
| 68 | +sum :: proc (v: Vec3) -> f64 { |
| 69 | + return v[0] + v[1] + v[2] |
| 70 | +} |
| 71 | + |
| 72 | +energy :: proc(sys: ^Sys) -> f64 { |
| 73 | + e := 0.0 |
| 74 | + for b, i in sys { |
| 75 | + e += 0.5 * b.mass * sum(b.velocity * b.velocity) |
| 76 | + for j in (i+1)..<5 { |
| 77 | + b2 := sys[j] |
| 78 | + dpos := b.pos - b2.pos |
| 79 | + distance_square := sum(dpos * dpos) |
| 80 | + distance := math.sqrt_f64(distance_square) |
| 81 | + e -= b.mass * b2.mass / distance |
| 82 | + } |
| 83 | + } |
| 84 | + return e |
| 85 | +} |
| 86 | + |
| 87 | +main :: proc() { |
| 88 | + args := runtime.args__ |
| 89 | + n := strconv.parse_int(auto_cast args[1]) or_else 1000 |
| 90 | + |
| 91 | + sys := init_sys() |
| 92 | + offset_momentum(&sys) |
| 93 | + fmt.printf("%.9f\n", energy(&sys)) |
| 94 | + for _ in 0..<n { |
| 95 | + advance(&sys, 0.01) |
| 96 | + } |
| 97 | + fmt.printf("%.9f\n", energy(&sys)) |
| 98 | +} |
0 commit comments