C++ implementation of NBody simulation, example code.

Simulation.cpp

float mult, dist, dX, dY, dZ;
Body *b;
Body *n;
Vector3f closing;

for (std::size_t i = 0; i < bodies->size(); i++) {
	b = bodies->at(i);
	closing.x = closing.y = closing.z = 0; // Reset force computation
	for (std::size_t j = 0; j < bodies->size(); j++) {
		n = bodies->at(j);
		if (i != j) {
			dist = n->getPosition()->dist(b->getPosition())); // Get the distance.
			if (dist > 0.0) {
				dX = (b->getPosition()->x - n->getPosition()->x) / dist; // Normalized changes
				dY = (b->getPosition()->y - n->getPosition()->y) / dist; // Exploded directional
				dZ = (b->getPosition()->z - n->getPosition()->z) / dist; // Unit vector b->n

				// Force between two masses
				mult = -GRAVITATIONAL_CONSTANT * n->getMass() * b->getMass() / (dist * dist);

				closing.add(dX * mult, dY * mult, dZ * mult);	// Add to computed force.
			}
		}
	}
	b->getVelocity()->add(closing.multiply(deltaT / b->getMass());	// Integrate accleration into velocity
	b->getPosition()->add(b->getVelocity()->clone()->multiply(deltaT));	// Integrate velocity into position
}

Vector3f.cpp

#include "Vector3f.h"
#include <math.h>

Vector3f::Vector3f(float x, float y, float z) {
	this->x = x;
	this->y = y;
	this->z = z;
}

Vector3f * Vector3f::add(float x, float y, float z) {
	this->x += x;
	this->y += y;
	this->z += z;
	return this;
}

float Vector3f::chunk(int i) {
	switch (i) {
	case 0:
		return x;
	case 1:
		return y;
	case 2:
		return z;
	default:
		return 0;
	}
}
void Vector3f::setChunk(int i, float v) {
	switch (i) {
	case 0:
		x = v;
		return;
	case 1:
		y = v;
		return;
	case 2:
		z = v;
		return;
	}
}

Vector3f * Vector3f::subtract(Vector3f * p) {
	this->x -= p->x;
	this->y -= p->y;
	this->z -= p->z;
	return this;
}

Vector3f * Vector3f::add(Vector3f * p) {
	this->x += p->x;
	this->y += p->y;
	this->z += p->z;
	return this;
}

Vector3f * Vector3f::multiply(float scalar) {
	this->x *= scalar;
	this->y *= scalar;
	this->z *= scalar;
	return this;
}

float Vector3f::dist(Vector3f * p) {
	return sqrt(distSquare(p));
}

float Vector3f::distSquare(Vector3f * p) {
	float dX = p->x - x, dY = p->y - y, dZ = p->z - z;
	return (dX * dX) + (dY * dY) + (dZ * dZ);
}

Vector3f * Vector3f::clone() {
	return new Vector3f(x, y, z);
}

Vector3f * Vector3f::normalize() {
	float mag = magnitude();
	if (mag != 0) {
		x /= mag;
		y /= mag;
		z /= mag;
	}
	return this;
}

Vector3f * Vector3f::negative() {
	x = -x;
	y = -y;
	z = -z;
	return this;
}

float Vector3f::magnitude() {
	return dist(new Vector3f(0, 0, 0));
}

float Vector3f::mag2() {
	return (x * x) + (y * y) + (z * z);
}