1/* 2 * Copyright (c) 2022-2024 Huawei Device Co., Ltd. 3 * Licensed under the Apache License, Version 2.0 (the "License"); 4 * you may not use this file except in compliance with the License. 5 * You may obtain a copy of the License at 6 * 7 * http://www.apache.org/licenses/LICENSE-2.0 8 * 9 * Unless required by applicable law or agreed to in writing, software 10 * distributed under the License is distributed on an "AS IS" BASIS, 11 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 * See the License for the specific language governing permissions and 13 * limitations under the License. 14 */ 15 16import { PI } from "std/math"; 17 18// readonly expected : double = -1.3524862408537381; 19class Body { 20 readonly SOLAR_MASS : double = 4 * PI * PI; 21 x : double ; 22 y : double ; 23 z : double ; 24 vx : double ; 25 vy : double ; 26 vz : double ; 27 mass : double ; 28 public constructor(x : double, y : double, z : double, vx : double, vy : double, vz : double, mass : double) { 29 this.x = x; 30 this.y = y; 31 this.z = z; 32 this.vx = vx; 33 this.vy = vy; 34 this.vz = vz; 35 this.mass = mass; 36 } 37 38 public offsetMomentum(px : double, py : double, pz : double): Body { 39 this.vx = -px / this.SOLAR_MASS; 40 this.vy = -py / this.SOLAR_MASS; 41 this.vz = -pz / this.SOLAR_MASS; 42 return this; 43 } 44} 45 46class NBodySystem { 47 bodies : Body[] ; 48 public constructor(bodies : Body[]) { 49 this.bodies = bodies; 50 let px : double = 0.0; 51 let py : double = 0.0; 52 let pz : double = 0.0; 53 let size : int = this.bodies.length; 54 for (let i : int = 0; i < size; i++) { 55 let b : Body = this.bodies[i]; 56 let m : double = b.mass; 57 px += b.vx * m; 58 py += b.vy * m; 59 pz += b.vz * m; 60 } 61 this.bodies[0].offsetMomentum(px, py, pz); 62 } 63 64 public advance(dt : double): void { 65 let dx : double ; 66 let dy : double ; 67 let dz : double ; 68 let distance : double ; 69 let mag : double ; 70 let size : int = this.bodies.length; 71 for (let i : int = 0; i < size; i++) { 72 let bodyi : Body = this.bodies[i]; 73 for (let j : int = i + 1; j < size; j++) { 74 let bodyj : Body = this.bodies[j]; 75 dx = bodyi.x - bodyj.x; 76 dy = bodyi.y - bodyj.y; 77 dz = bodyi.z - bodyj.z; 78 distance = sqrt(dx * dx + dy * dy + dz * dz); 79 mag = dt / (distance * distance * distance); 80 bodyi.vx -= dx * bodyj.mass * mag; 81 bodyi.vy -= dy * bodyj.mass * mag; 82 bodyi.vz -= dz * bodyj.mass * mag; 83 bodyj.vx += dx * bodyi.mass * mag; 84 bodyj.vy += dy * bodyi.mass * mag; 85 bodyj.vz += dz * bodyi.mass * mag; 86 } 87 } 88 for (let i : int = 0; i < size; i++) { 89 let body : Body = this.bodies[i]; 90 body.x += dt * body.vx; 91 body.y += dt * body.vy; 92 body.z += dt * body.vz; 93 } 94 } 95 public energy(): double { 96 let dx : double ; 97 let dy : double ; 98 let dz : double ; 99 let distance : double ; 100 let e : double = 0.0; 101 let size : int = this.bodies.length; 102 for (let i : int = 0; i < size; i++) { 103 let bodyi : Body = this.bodies[i]; 104 e += 0.5 * bodyi.mass * (bodyi.vx * bodyi.vx + bodyi.vy * bodyi.vy + bodyi.vz * bodyi.vz); 105 for (let j : int = i + 1; j < size; j++) { 106 let bodyj : Body = this.bodies[j]; 107 dx = bodyi.x - bodyj.x; 108 dy = bodyi.y - bodyj.y; 109 dz = bodyi.z - bodyj.z; 110 distance = sqrt(dx * dx + dy * dy + dz * dz); 111 e -= (bodyi.mass * bodyj.mass) / distance; 112 } 113 } 114 return e; 115 } 116 } 117 118export class AccessNBody { 119 static readonly SOLAR_MASS : double = 4 * PI * PI; 120 static readonly DAYS_PER_YEAR : double = 365.24; 121 n1 : int = 3; 122 n2 : int = 24; 123 124 static jupiter(): Body { 125 return new Body(4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01, 1.66007664274403694e-03 * AccessNBody.DAYS_PER_YEAR, 7.69901118419740425e-03 * AccessNBody.DAYS_PER_YEAR, -6.90460016972063023e-05 * AccessNBody.DAYS_PER_YEAR, 9.54791938424326609e-04 * AccessNBody.SOLAR_MASS); 126 } 127 static saturn(): Body { 128 return new Body(8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01, -2.76742510726862411e-03 * AccessNBody.DAYS_PER_YEAR, 4.99852801234917238e-03 * AccessNBody.DAYS_PER_YEAR, 2.30417297573763929e-05 * AccessNBody.DAYS_PER_YEAR, 2.85885980666130812e-04 * AccessNBody.SOLAR_MASS); 129 } 130 static uranus(): Body { 131 return new Body(1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01, 2.96460137564761618e-03 * AccessNBody.DAYS_PER_YEAR, 2.37847173959480950e-03 * AccessNBody.DAYS_PER_YEAR, -2.96589568540237556e-05 * AccessNBody.DAYS_PER_YEAR, 4.36624404335156298e-05 * AccessNBody.SOLAR_MASS); 132 } 133 static neptune(): Body { 134 return new Body(1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01, 2.68067772490389322e-03 * AccessNBody.DAYS_PER_YEAR, 1.62824170038242295e-03 * AccessNBody.DAYS_PER_YEAR, -9.51592254519715870e-05 * AccessNBody.DAYS_PER_YEAR, 5.15138902046611451e-05 * AccessNBody.SOLAR_MASS); 135 } 136 static sun(): Body { 137 return new Body(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, AccessNBody.SOLAR_MASS); 138 } 139 readonly expected : double = -1.3524862408537381; 140 public run(): void { 141 let ret : double = 0; 142 for (let n : int = this.n1; n <= this.n2; n *= 2) { 143 let bodies : NBodySystem = new NBodySystem([AccessNBody.sun(), AccessNBody.jupiter(), AccessNBody.saturn(), AccessNBody.uranus(), AccessNBody.neptune()]); 144 let max : int = n * 100; 145 ret += bodies.energy(); 146 for (let i : int = 0; i < max; i++) { 147 bodies.advance(0.01); 148 } 149 ret += bodies.energy(); 150 } 151 152 assert ret == this.expected: "Incorrect result"; 153 } 154} 155 156function main(): void { 157 let a = new AccessNBody; 158 a.run(); 159} 160