An error occurred while loading the file. Please try again.
-
Paul Heinzlreiter authored28dfbdb7
Body.cpp 4.05 KiB
#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;
}
}