Body.cpp

#include "Body.hpp"
#include <stdlib.h>
#include <iostream>
#include <cmath>
#include <cfloat>

namespace nbody {
	using namespace std;

	class Derivative {
		friend class Body;
	private:
		double dx[3];
		double dv[3];

		Derivative() {
			for (int i = 0; i < 3; i++) {
				this->dx[i] = 0.0;
				this->dv[i] = 0.0;
			}
		}

		~Derivative() {}
	};

	Body::Body() {
		for (int i = 0; i < 3; i++) {
			this->position[i] = 0.0;
			this->velocity[i] = 0.0;
			this->acceleration[i] = 0.0;
		}
		this->mass = 1.0;
		this->refinement = 0;
	}

	Body::Body(double positionX, double positionY, double positionZ) {
		for (int i = 0; i < 3; i++) {
			this->velocity[i] = 0.0;
			this->acceleration[i] = 0.0;
		}
		this->position[0] = positionX;
		this->position[1] = positionY;
		this->position[2] = positionZ;
		this->mass = 1.0;
		this->refinement = 0;
	}

	Body::Body(double positionX, double positionY, double positionZ, double velocityX, double velocityY, double velocityZ) {
		for (int i = 0; i < 3; i++) {
			this->acceleration[i] = 0.0;
		}
		this->position[0] = positionX;
		this->position[1] = positionY;
		this->position[2] = positionZ;
		this->velocity[0] = velocityX;
		this->velocity[1] = velocityY;
		this->velocity[2] = velocityZ;
		this->mass = 1.0;
		this->refinement = 0;
	}

	Body::Body(double positionX, double positionY, double positionZ, double weight) {
		for (int i = 0; i < 3; i++) {
			this->velocity[i] = 0.0;
			this->acceleration[i] = 0.0;
		}
		this->position[0] = positionX;
		this->position[1] = positionY;
		this->position[2] = positionZ;
		this->mass = weight;
		this->refinement = 0;
	}

	Body::Body(double positionX, double positionY, double positionZ, double velocityX, double velocityY, double velocityZ, double weight) {
		for (int i = 0; i < 3; i++) {
			this->acceleration[i] = 0.0;
		}
		this->position[0] = positionX;
		this->position[1] = positionY;
		this->position[2] = positionZ;
		this->velocity[0] = velocityX;
		this->velocity[1] = velocityY;
		this->velocity[2] = velocityZ;
		this->mass = weight;
		this->refinement = 0;
	}

	Body::~Body() {

	}

	Derivative Body::evaluate(double dt, Derivative d) {
		double velocity[3];
		Derivative output;

		for (int i = 0; i < 3; i++) {
			velocity[i] = this->velocity[i] + d.dv[i] * dt;
			output.dx[i] = velocity[i];
			output.dv[i] = this->acceleration[i];
		}
		return output;
	}

	void Body::integrate() {
		Derivative start;
		Derivative a, b, c, d;
		double dxdt[3];
		double dvdt[3];

		if (this->refinement) {
			return;
		}
		a = evaluate(0.0, start);
		b = evaluate(timestep * 0.5, a);
		c = evaluate(timestep * 0.5, b);
		d = evaluate(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]);
			this->position[i] = this->position[i] + dxdt[i] * timestep;
			this->velocity[i] = this->velocity[i] + dvdt[i] * timestep;
		}
		//this->print();
	}

	void Body::accumulateForceOnto(Body& to) {
		double distance2 = 0.0;

		distance2 += (this->position[0] - to.position[0]) * (this->position[0] - to.position[0]);
		distance2 += (this->position[1] - to.position[1]) * (this->position[1] - to.position[1]);
		distance2 += (this->position[2] - to.position[2]) * (this->position[2] - to.position[2]);
		if (fabs(distance2) < FLT_EPSILON) {
			distance2 = FLT_EPSILON;
		}
		double distance = sqrt(distance2);
		double mdist = -1.0 * ((this->mass * to.mass) / distance2);
		for (int i = 0; i < 3; i++) {
			double vec = (this->position[i] - to.position[i]) / distance;

			to.acceleration[i] += vec * mdist;
		}
	}

	void Body::resetAcceleration() {
		for (int i = 0; i < 3; i++) {
			this->acceleration[i] = 0.0;
		}
	}

	bool Body::valid() {

		for (int i = 0; i < 3; i++) {
			if (this->position[i] < -FLT_MAX || this->position[i] > FLT_MAX) {
				return false;
			}
		}
		return this->mass >= 0.0;
	}

	void Body::print() {
		cout << "Position: " << this->position[0] << " " << this->position[1] << " " << this->position[2] << endl;
		cout << "Velocity: " << this->velocity[0] << " " << this->velocity[1] << " " << this->velocity[2] << endl;
		cout << "Acceleration: " << this->acceleration[0] << " " << this->acceleration[1] << " " << this->acceleration[2] << endl;
		cout << "Mass: " << this->mass << endl;
		cout << "Refinement: " << this->refinement << endl << endl;
	}
}