1 | // Ported based on https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/nbody-go-3.html |
2 | // Output: |
3 | // -0.169075164 |
4 | // -0.169059907 |
5 | import math |
6 | |
7 | const ( |
8 | solar_mass = 39.47841760435743197 // 4.0 * math.Pi * math.Pi |
9 | days_per_year = 365.24 |
10 | c_n = 5 |
11 | ) |
12 | |
13 | struct Position { |
14 | pub mut: |
15 | x f64 |
16 | y f64 |
17 | z f64 |
18 | } |
19 | |
20 | struct Momentum { |
21 | pub mut: |
22 | x f64 |
23 | y f64 |
24 | z f64 |
25 | m f64 |
26 | } |
27 | |
28 | struct System { |
29 | pub mut: |
30 | v []Momentum |
31 | s []Position |
32 | } |
33 | |
34 | fn advance(mut sys System, dt f64) { |
35 | for i in 0 .. c_n - 1 { |
36 | mut vx := sys.v[i].x |
37 | mut vy := sys.v[i].y |
38 | mut vz := sys.v[i].z |
39 | |
40 | for j := i + 1; j < c_n; j++ { |
41 | dx := sys.s[i].x - sys.s[j].x |
42 | dy := sys.s[i].y - sys.s[j].y |
43 | dz := sys.s[i].z - sys.s[j].z |
44 | |
45 | dsquared := dx * dx + dy * dy + dz * dz |
46 | distance := math.sqrt(dsquared) |
47 | mag := (dt / (dsquared * distance)) |
48 | mi := sys.v[i].m |
49 | |
50 | vx -= dx * sys.v[j].m * mag |
51 | vy -= dy * sys.v[j].m * mag |
52 | vz -= dz * sys.v[j].m * mag |
53 | |
54 | sys.v[j].x += dx * mi * mag |
55 | sys.v[j].y += dy * mi * mag |
56 | sys.v[j].z += dz * mi * mag |
57 | } |
58 | sys.v[i].x = vx |
59 | sys.v[i].y = vy |
60 | sys.v[i].z = vz |
61 | } |
62 | |
63 | for i in 0 .. c_n { |
64 | sys.s[i].x += dt * sys.v[i].x |
65 | sys.s[i].y += dt * sys.v[i].y |
66 | sys.s[i].z += dt * sys.v[i].z |
67 | } |
68 | } |
69 | |
70 | fn offsetmomentum(mut sys System) { |
71 | mut px := f64(0) |
72 | mut py := f64(0) |
73 | mut pz := f64(0) |
74 | |
75 | for i in 0 .. c_n { |
76 | px += sys.v[i].x * sys.v[i].m |
77 | py += sys.v[i].y * sys.v[i].m |
78 | pz += sys.v[i].z * sys.v[i].m |
79 | } |
80 | sys.v[0].x = -px / solar_mass |
81 | sys.v[0].y = -py / solar_mass |
82 | sys.v[0].z = -pz / solar_mass |
83 | } |
84 | |
85 | fn energy(sys System) f64 { |
86 | mut e := f64(0) |
87 | for i in 0 .. c_n { |
88 | e += 0.5 * sys.v[i].m * (sys.v[i].x * sys.v[i].x + sys.v[i].y * sys.v[i].y + |
89 | sys.v[i].z * sys.v[i].z) |
90 | for j := i + 1; j < c_n; j++ { |
91 | dx := sys.s[i].x - sys.s[j].x |
92 | dy := sys.s[i].y - sys.s[j].y |
93 | dz := sys.s[i].z - sys.s[j].z |
94 | distance := math.sqrt(dx * dx + dy * dy + dz * dz) |
95 | e -= (sys.v[i].m * sys.v[j].m) / distance |
96 | } |
97 | } |
98 | return e |
99 | } |
100 | |
101 | fn arr_momentum() []Momentum { |
102 | return [ |
103 | Momentum{0.0, 0.0, 0.0, solar_mass}, |
104 | Momentum{1.66007664274403694e-03 * days_per_year, 7.69901118419740425e-03 * days_per_year, -6.90460016972063023e-05 * days_per_year, 9.54791938424326609e-04 * solar_mass}, |
105 | Momentum{-2.76742510726862411e-03 * days_per_year, 4.99852801234917238e-03 * days_per_year, 2.30417297573763929e-05 * days_per_year, 2.85885980666130812e-04 * solar_mass}, |
106 | Momentum{2.96460137564761618e-03 * days_per_year, 2.37847173959480950e-03 * days_per_year, -2.96589568540237556e-05 * days_per_year, 4.36624404335156298e-05 * solar_mass}, |
107 | Momentum{2.68067772490389322e-03 * days_per_year, 1.62824170038242295e-03 * days_per_year, -9.51592254519715870e-05 * days_per_year, 5.15138902046611451e-05 * solar_mass}, |
108 | ] |
109 | } |
110 | |
111 | fn arr_position() []Position { |
112 | return [ |
113 | Position{0.0, 0.0, 0.0}, |
114 | Position{4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01}, |
115 | Position{8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01}, |
116 | Position{1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01}, |
117 | Position{1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01}, |
118 | ] |
119 | } |
120 | |
121 | fn main() { |
122 | mut sys := &System{arr_momentum(), arr_position()} |
123 | offsetmomentum(mut sys) |
124 | println('${energy(sys):.9f}') //-0.169075164 |
125 | for _ in 0 .. 50_000_000 { |
126 | advance(mut sys, 0.01) |
127 | } |
128 | println('${energy(sys):.9f}') //-0.169059907 |
129 | } |