Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
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
105
106
107
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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
#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;
}
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) {
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;
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;
}
}