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) {
for (int i = 0; i < 3; i++) {
body.acceleration[i] = 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;
}
}
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
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) {
char buf[128];
sprintf(buf, "%f", val);
for (int i = 0; i < strlen(buf); i++) {
if (!((buf[i] >= '0' && buf[i] <= '9') || buf[i] == '.' || buf[i] == '-')) {
return false;
}
}
return true;
}
bool valid(Body body) {
for (int i = 0; i < 3; i++) {
if (!validDouble(body.position[i])) return false;
if (!validDouble(body.velocity[i])) return false;
if (!validDouble(body.acceleration[i])) 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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
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);
}
infile.close();
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) {
for (vector<Body>::iterator it = bodies.begin(); it != bodies.end(); it++) {
if (!valid(*it)) {
return false;
}
}
return true;
}
}