Newer
Older
#include "Body.hpp"
#include <stdlib.h>
#include <iostream>
#include <sstream>
#include <fstream>
#include <cmath>
#include <cfloat>
#include <climits>
#include <cstring>
#include <cstdio>
namespace nbody {
using namespace std;
void resetAcceleration(Body& body) {
fill(body.acceleration.begin(), body.acceleration.end(), 0.0);
Derivative evaluate(Body body, double dt, Derivative d) {
double velocity[3];
Derivative output;
for (int i = 0; i < 3; i++) {
velocity[i] = body.velocity[i] + d.dv[i] * dt;
output.dx[i] = velocity[i];
output.dv[i] = body.acceleration[i];
}
return output;
}
//integrating acceleration -> velocity -> position
void integrate(Body& body) {
Derivative start;
Derivative a, b, c, d;
double dxdt[3];
double dvdt[3];
if (body.refinement) {
return;
}
a = evaluate(body, 0.0, start);
b = evaluate(body, timestep * 0.5, a);
c = evaluate(body, timestep * 0.5, b);
d = evaluate(body, timestep, c);
for (int i = 0; i < 3; i++) {
dxdt[i] = 1.0 / 6.0 * (a.dx[i] + 2.0f * (b.dx[i] + c.dx[i]) + d.dx[i]);
dvdt[i] = 1.0 / 6.0 * (a.dv[i] + 2.0f * (b.dv[i] + c.dv[i]) + d.dv[i]);
body.position[i] = body.position[i] + dxdt[i] * timestep;
body.velocity[i] = body.velocity[i] + dvdt[i] * timestep;
}
}
void accumulateForceOntoBody(Body& to, Body from) {
double distance2 = 0.0;
distance2 += (from.position[0] - to.position[0]) * (from.position[0] - to.position[0]);
distance2 += (from.position[1] - to.position[1]) * (from.position[1] - to.position[1]);
distance2 += (from.position[2] - to.position[2]) * (from.position[2] - to.position[2]);
if (fabs(distance2) < FLT_EPSILON) {
distance2 = FLT_EPSILON;
}
double distance = sqrt(distance2);
double mdist = -1.0 * ((from.mass * to.mass) / distance2);
for (int i = 0; i < 3; i++) {
double vec = (from.position[i] - to.position[i]) / distance;
to.acceleration[i] += vec * mdist;
}
}
bool validDouble(double val) {
return !isnan(val) && isfinite(val);
bool validBody(Body body) {
if (!all_of(body.position.begin(), body.position.end(), validDouble)) return false;
if (!all_of(body.velocity.begin(), body.velocity.end(), validDouble)) return false;
if (!all_of(body.acceleration.begin(), body.acceleration.end(), validDouble)) return false;
if (!validDouble(body.mass)) return false;
return body.mass >= 0.0;
}
void printBody(int parallelId, Body body) {
cout << parallelId << " " << body.id << " Position: " << body.position[0] << " " << body.position[1] << " " << body.position[2] << endl;
cout << parallelId << " " << body.id << " Velocity: " << body.velocity[0] << " " << body.velocity[1] << " " << body.velocity[2] << endl;
cout << parallelId << " " << body.id << " Acceleration: " << body.acceleration[0] << " " << body.acceleration[1] << " " << body.acceleration[2] << endl;
cout << parallelId << " " << body.id << " Mass: " << body.mass << endl;
cout << parallelId << " " << body.id << " Refinement: " << body.refinement << endl << endl;
}
vector<Body> dubinskiParse(string filename) {
vector<Body> result;
string line;
ifstream infile(filename.c_str(), ifstream::in);
double mass, px, py, pz, vx, vy, vz;
unsigned long id = 0;
while (std::getline(infile, line)) {
istringstream iss(line);
Body b;
//not all input files have velocity, so initialize properly
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
vx = vy = vz = 0.0;
iss >> mass >> px >> py >> pz >> vx >> vy >> vz;
b.position[0] = px;
b.position[1] = py;
b.position[2] = pz;
b.velocity[0] = vx;
b.velocity[1] = vy;
b.velocity[2] = vz;
b.acceleration[0] = 0.0;
b.acceleration[1] = 0.0;
b.acceleration[2] = 0.0;
b.refinement = 0;
b.mass = mass;
b.id = id++;
result.push_back(b);
}
return result;
}
void printBodies(int parallelId, vector<Body> bodies) {
for (vector<Body>::iterator it = bodies.begin(); it != bodies.end(); it++) {
printBody(parallelId, *it);
}
}
bool valid(vector<Body> bodies) {
return all_of(bodies.begin(), bodies.end(), validBody);