|
1 | 1 |
|
2 | | -// Simplified N-Body |
3 | 2 | class Body { |
4 | | - init(x, y, z, mass) { |
| 3 | + init(x, y, z, vx, vy, vz, mass) { |
5 | 4 | this.x = x; |
6 | 5 | this.y = y; |
7 | 6 | this.z = z; |
| 7 | + this.vx = vx; |
| 8 | + this.vy = vy; |
| 9 | + this.vz = vz; |
8 | 10 | this.mass = mass; |
9 | | - this.vx = 0; |
10 | | - this.vy = 0; |
11 | | - this.vz = 0; |
12 | 11 | } |
13 | 12 | } |
14 | 13 |
|
15 | | -func advance(bodies, dt) { |
16 | | - let n = length(bodies); |
17 | | - for (let i = 0; i < n; i = i + 1) { |
18 | | - for (let j = i + 1; j < n; j = j + 1) { |
19 | | - let b1 = bodies[i]; |
20 | | - let b2 = bodies[j]; |
21 | | - |
22 | | - let dx = b1.x - b2.x; |
23 | | - let dy = b1.y - b2.y; |
24 | | - let dz = b1.z - b2.z; |
25 | | - |
26 | | - let dist_sq = dx*dx + dy*dy + dz*dz; |
27 | | - let mag = dt / (dist_sq * sqrt(dist_sq)); |
28 | | - |
29 | | - b1.vx = b1.vx - dx * b2.mass * mag; |
30 | | - b1.vy = b1.vy - dy * b2.mass * mag; |
31 | | - b1.vz = b1.vz - dz * b2.mass * mag; |
32 | | - |
33 | | - b2.vx = b2.vx + dx * b1.mass * mag; |
34 | | - b2.vy = b2.vy + dy * b1.mass * mag; |
35 | | - b2.vz = b2.vz + dz * b1.mass * mag; |
| 14 | +class NBodySystem { |
| 15 | + init() { |
| 16 | + this.bodies = []; |
| 17 | + // Sun |
| 18 | + list_push(this.bodies, Body(0, 0, 0, 0, 0, 0, 1.98892e30)); |
| 19 | + // Jupiter |
| 20 | + list_push(this.bodies, Body(4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01, |
| 21 | + 1.66007664274403694e-03 * 365.24, 7.69901118419740425e-03 * 365.24, -6.90460016972063023e-05 * 365.24, |
| 22 | + 9.54791938424326609e-04 * 1.98892e30)); |
| 23 | + // Additional bodies can be added... |
| 24 | + } |
| 25 | + |
| 26 | + offset_momentum() { |
| 27 | + let px = 0; |
| 28 | + let py = 0; |
| 29 | + let pz = 0; |
| 30 | + let n = len(this.bodies); |
| 31 | + for (let i = 0; i < n; i = i + 1) { |
| 32 | + let b = this.bodies[i]; |
| 33 | + px = px + b.vx * b.mass; |
| 34 | + py = py + b.vy * b.mass; |
| 35 | + pz = pz + b.vz * b.mass; |
36 | 36 | } |
| 37 | + let sol = this.bodies[0]; |
| 38 | + sol.vx = -px / sol.mass; |
| 39 | + sol.vy = -py / sol.mass; |
| 40 | + sol.vz = -pz / sol.mass; |
37 | 41 | } |
38 | | - |
39 | | - for (let i = 0; i < n; i = i + 1) { |
40 | | - let b = bodies[i]; |
41 | | - b.x = b.x + b.vx * dt; |
42 | | - b.y = b.y + b.vy * dt; |
43 | | - b.z = b.z + b.vz * dt; |
| 42 | + |
| 43 | + energy() { |
| 44 | + let e = 0.0; |
| 45 | + let n = len(this.bodies); |
| 46 | + for (let i = 0; i < n; i = i + 1) { |
| 47 | + let bi = this.bodies[i]; |
| 48 | + e = e + 0.5 * bi.mass * (bi.vx * bi.vx + bi.vy * bi.vy + bi.vz * bi.vz); |
| 49 | + for (let j = i + 1; j < n; j = j + 1) { |
| 50 | + let bj = this.bodies[j]; |
| 51 | + let dx = bi.x - bj.x; |
| 52 | + let dy = bi.y - bj.y; |
| 53 | + let dz = bi.z - bj.z; |
| 54 | + let dist = sqrt(dx*dx + dy*dy + dz*dz); |
| 55 | + e = e - (bi.mass * bj.mass) / dist; |
| 56 | + } |
| 57 | + } |
| 58 | + return e; |
| 59 | + } |
| 60 | + |
| 61 | + advance(dt) { |
| 62 | + let n = len(this.bodies); |
| 63 | + for (let i = 0; i < n; i = i + 1) { |
| 64 | + let bi = this.bodies[i]; |
| 65 | + for (let j = i + 1; j < n; j = j + 1) { |
| 66 | + let bj = this.bodies[j]; |
| 67 | + let dx = bi.x - bj.x; |
| 68 | + let dy = bi.y - bj.y; |
| 69 | + let dz = bi.z - bj.z; |
| 70 | + |
| 71 | + let dist_sq = dx*dx + dy*dy + dz*dz; |
| 72 | + let mag = dt / (dist_sq * sqrt(dist_sq)); |
| 73 | + |
| 74 | + let bim = bi.mass * mag; |
| 75 | + let bjm = bj.mass * mag; |
| 76 | + |
| 77 | + bi.vx = bi.vx - dx * bjm; |
| 78 | + bi.vy = bi.vy - dy * bjm; |
| 79 | + bi.vz = bi.vz - dz * bjm; |
| 80 | + |
| 81 | + bj.vx = bj.vx + dx * bim; |
| 82 | + bj.vy = bj.vy + dy * bim; |
| 83 | + bj.vz = bj.vz + dz * bim; |
| 84 | + } |
| 85 | + } |
| 86 | + |
| 87 | + for (let i = 0; i < n; i = i + 1) { |
| 88 | + let b = this.bodies[i]; |
| 89 | + b.x = b.x + b.vx * dt; |
| 90 | + b.y = b.y + b.vy * dt; |
| 91 | + b.z = b.z + b.vz * dt; |
| 92 | + } |
44 | 93 | } |
45 | 94 | } |
46 | 95 |
|
47 | 96 | func main() { |
48 | 97 | let start = clock(); |
49 | 98 | let system = NBodySystem(); |
50 | | - for (let i = 0; i < 10000; i = i + 1) { |
| 99 | + system.offset_momentum(); |
| 100 | + print("Initial Energy: " + to_string(system.energy())); |
| 101 | + |
| 102 | + for (let i = 0; i < 1000; i = i + 1) { |
51 | 103 | system.advance(0.01); |
52 | 104 | } |
| 105 | + |
53 | 106 | let elapsed = clock() - start; |
54 | 107 |
|
55 | | - print("Energy: " + to_string(system.energy())); |
| 108 | + print("Final Energy: " + to_string(system.energy())); |
56 | 109 | print("Time: " + to_string(elapsed)); |
57 | 110 | } |
58 | 111 |
|
|
0 commit comments