// nbody.cpp - compiled by csrcpp.srp from nbody.csr #include "csrstdlib.h" // The Computer Language Benchmarks Game; // https://salsa.debian.org/benchmarksgame-team/benchmarksgame/; // originally by Kevin Carson; // modified by Tupteq, Fredrik Johansson, and Daniel Nanz; // modified by Maciej Fijalkowski; // 2to3; class Body_class : public RCBase { public: Real x, y, z; Real vx, vy, vz; Real m; void init(Real x_, Real y_, Real z_, Real vx_, Real vy_, Real vz_, Real m_) {; x = x_; y = y_; z = z_; vx = vx_; vy = vy_; vz = vz_; m = m_; } }; typedef RCPtr Body; Body Body_create(Real x_, Real y_, Real z_, Real vx_, Real vy_, Real vz_, Real m_) { Body rcptr(new Body_class); rcptr->init(x_, y_, z_, vx_, vy_, vz_, m_); return rcptr; } class Pair_class : public RCBase { public: Body b1; Body b2; void init(Body b1_, Body b2_) {; b1 = b1_; b2 = b2_; } }; typedef RCPtr Pair; Pair Pair_create(Body b1_, Body b2_) { Pair rcptr(new Pair_class); rcptr->init(b1_, b2_); return rcptr; } // typedef Bodies = array_of(Body) typedef RCPtrA, Body > Bodies; Bodies Bodies_create(Int n = 0, Body val = 0) { Bodies rcptr(new Array_class(n, val)); return rcptr; } Int len(Bodies ar) { return ar->size(); } // typedef Pairs = array_of(Pair) typedef RCPtrA, Pair > Pairs; Pairs Pairs_create(Int n = 0, Pair val = 0) { Pairs rcptr(new Array_class(n, val)); return rcptr; } Int len(Pairs ar) { return ar->size(); } // typedef String_Body_di = dict_of(String, Body) typedef RCPtrD, String, Body > String_Body_di; String_Body_di String_Body_di_create() { String_Body_di rcptr(new Dict_class); return rcptr; } Pairs combinations(Bodies l) { Pairs result = Pairs_create(0, 0); for (Int i = 0; i < len(l) - 1; i += 1) { Body x = l[i]; for (Int j = i + 1; j < len(l); j += 1) { result->append(Pair_create(l[i], l[j])); } } return result; } Real PI; Real SOLAR_MASS; Real DAYS_PER_YEAR; String_Body_di BODIES; Bodies SYSTEM; Pairs PAIRS; void advance(Real dt, Int n) { for (Int i = 0; i < n; i += 1) { for (int _i = 0; _i < PAIRS->size(); _i++) { Pair pair = (*(PAIRS))[_i]; Body b1 = pair->b1; Body b2 = pair->b2; Real dx = b1->x - b2->x; Real dy = b1->y - b2->y; Real dz = b1->z - b2->z; Real mag = dt * pow(dx * dx + dy * dy + dz * dz, -1.5); Real b1m = b1->m * mag; Real b2m = b2->m * mag; b1->vx = b1->vx - dx * b2m; b1->vy = b1->vy - dy * b2m; b1->vz = b1->vz - dz * b2m; b2->vx = b2->vx + dx * b1m; b2->vy = b2->vy + dy * b1m; b2->vz = b2->vz + dz * b1m; } for (int _i = 0; _i < SYSTEM->size(); _i++) { Body body = (*(SYSTEM))[_i]; body->x = body->x + dt * body->vx; body->y = body->y + dt * body->vy; body->z = body->z + dt * body->vz; } } } void report_energy() { Real e = 0.0; for (int _i = 0; _i < PAIRS->size(); _i++) { Pair pair = (*(PAIRS))[_i]; Body b1 = pair->b1; Body b2 = pair->b2; Real dx = b1->x - b2->x; Real dy = b1->y - b2->y; Real dz = b1->z - b2->z; e = e - (b1->m * b2->m) / pow(dx * dx + dy * dy + dz * dz, 0.5); } for (int _i = 0; _i < SYSTEM->size(); _i++) { Body body = (*(SYSTEM))[_i]; e = e + body->m * (body->vx * body->vx + body->vy * body->vy + body->vz * body->vz) / 2.0; } print(e, 0); print("\n", 0); } void offset_momentum(Body ref) { Real px = 0; Real py = 0; Real pz = 0; for (int _i = 0; _i < SYSTEM->size(); _i++) { Body body = (*(SYSTEM))[_i]; Real m = body->m; px = px - body->vx * m; py = py - body->vy * m; pz = pz - body->vz * m; } Real m = ref->m; ref->vx = px / m; ref->vy = py / m; ref->vz = pz / m; } void main_(Int n) { offset_momentum(BODIES[String("sun")]); report_energy(); advance(0.01, n); report_energy(); } int main(int argc, char **argv) { csrinit(argc, argv); //; PI = 3.14159265358979323; SOLAR_MASS = 4 * PI * PI; DAYS_PER_YEAR = 365.24; String_Body_di string_body_di_1 = String_Body_di_create(); string_body_di_1[String("sun")] = Body_create(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, SOLAR_MASS); string_body_di_1[String("jupiter")] = Body_create(4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01, 1.66007664274403694e-03 * DAYS_PER_YEAR, 7.69901118419740425e-03 * DAYS_PER_YEAR, -6.90460016972063023e-05 * DAYS_PER_YEAR, 9.54791938424326609e-04 * SOLAR_MASS); string_body_di_1[String("saturn")] = Body_create(8.34336671824457987e+00, 4.12479856412430479e+00, -4.03523417114321381e-01, -2.76742510726862411e-03 * DAYS_PER_YEAR, 4.99852801234917238e-03 * DAYS_PER_YEAR, 2.30417297573763929e-05 * DAYS_PER_YEAR, 2.85885980666130812e-04 * SOLAR_MASS); string_body_di_1[String("uranus")] = Body_create(1.28943695621391310e+01, -1.51111514016986312e+01, -2.23307578892655734e-01, 2.96460137564761618e-03 * DAYS_PER_YEAR, 2.37847173959480950e-03 * DAYS_PER_YEAR, -2.96589568540237556e-05 * DAYS_PER_YEAR, 4.36624404335156298e-05 * SOLAR_MASS); string_body_di_1[String("neptune")] = Body_create(1.53796971148509165e+01, -2.59193146099879641e+01, 1.79258772950371181e-01, 2.68067772490389322e-03 * DAYS_PER_YEAR, 1.62824170038242295e-03 * DAYS_PER_YEAR, -9.51592254519715870e-05 * DAYS_PER_YEAR, 5.15138902046611451e-05 * SOLAR_MASS) ; BODIES = string_body_di_1; SYSTEM = BODIES->values(); PAIRS = combinations(SYSTEM); print("nbody", 1); print(command_line_arguments[1], 0); print("\n", 0); main_(string_to_int(command_line_arguments[1])); exit(); return 0; }