N-body gravity simulator. There are two identical stars, rotating around the center-of-mass, and the light probe. Requires g++ >=4.6. Run "make".

main.c

#include "physics.h"

typedef enum { SUN, SUN1, PROBE } ID;

Object objects[] =
    { 
        { SUN, 100000
        , { -2 + 0.05,     0, 0 }
        , { -13000,   -10000, 0 } 
        }

    , 
        { SUN1, 100000
        , { 2 + 0.07,      0, 0 }
        , { -11000,    10000, 0 } 
    	}

    , 
        { PROBE, 1
        , { 0, 0, 0 }
        , { -12000 + 5900, 0, 0 } 
        }
};

int main() {

    step(objects, 3, 1000000);
}

main.c++

#include "physics.h++"

enum ID { SUN, SUN1, PROBE };

vector<Object> objects =
    { 
        { SUN, 100000
        , { -2 + 0.05,     0, 0 }
        , { -13000,   -10000, 0 } 
        }

    , 
        { SUN1, 100000
        , { 2 + 0.07,      0, 0 }
        , { -11000,    10000, 0 } 
    	}

    , 
        { PROBE, 1
        , { 0, 0, 0 }
        , { -12000 + 5900, 0, 0 } 
        }
};

int main() {

    // make_distinct(objects);

    step(1000000, objects);
}

Makefile

SHELL = /bin/bash

all: test-c test-cc test-cpp

clean:
	-rm c cpp cc

test-c: c
	@echo "running c version..."
	@time ./c
	@echo ""

test-cc: cc
	@echo "running c version build with g++..."
	@time ./cc
	@echo ""

test-cpp: cpp
	@echo "running c++ version..."
	@time ./cpp
	@echo ""

c: physics.c main.c
	@gcc -std=c11 physics.c -lm main.c -O2 -o c

*.c: physics.h

cc: physics.c main.c
	@g++ -std=c++11 physics.c -lm main.c -O2 -o cc

cpp: physics.c++ main.c++
	@g++ -std=c++11 physics.c++ main.c++ -O2 -o cpp

*.cpp: physics.h

physics.c

#include "physics.h"

void make_simulation_step(Object * item, REAL dt, const Object * objects, int count) {

    REAL a[] = { 0, 0, 0 };

    advance(item, dt, acceleration(item, a, objects, count));

}

void advance(Object * item, REAL t, REAL * a) {

    for (int i = 0; i < 3; i++) {

        item->x[i] += (item->v[i] += a[i] * t) * t;

    }
}

REAL * acceleration(Object * item, REAL * accel, const Object * objects, int count) {

    for (int i = 0; i < count; i++) {

        const Object * object = objects + i;

        if (item->id == object->id) continue;

        REAL offset[3], distance_sqr = 0;

        for (int i = 0; i < 3; i++) { 

            offset[i]     = object->x[i] - item->x[i];
            distance_sqr += offset[i] * offset[i];

        }

        for (int i = 0; i < 3; i++) {

            offset[i] /= sqrt(distance_sqr);
            accel [i] += offset[i] * object->mass / distance_sqr;

        }
    }

    return accel;
}

void step(Object * objects, int count, int steps) {

    for (int j = 0; j < steps; j++) {

        for (int i = 0; i < count; i++) {

            make_simulation_step(objects + i, 1, objects, count);

        }

    }

}

void make_distinct(Object * objects, int count) {

    for (int i = 0; i < count; ++i)
    {
        objects[i].id = i;
    }

}

physics.c++

#include "physics.h++"

#define times(iter, n) for (int iter = 0; iter < n; iter++)

template <class T>
inline T sqr(T x) { return x * x; }

void Object::advance(REAL t, vector<REAL> a) {

    times(i, 3) {

        v[i] += a[i] * t;
        x[i] += v[i] * t;

    }

}

vector<REAL> Object::acceleration(const vector<Object> &objects) const {

    vector<REAL> accel = { 0, 0, 0 };

    for (auto const &object : objects) {
        
        if (id == object.id) continue;

        REAL offset[3], distance_sqr = 0;

        times(i, 3) { 

            offset[i]     = object.x[i] - x[i];
            distance_sqr += sqr(offset[i]);

        }

        times(i, 3) {

            offset[i] /= sqrt(distance_sqr);
            accel [i] += offset[i] * object.mass / distance_sqr;

        }

    }

    return accel;
}

void Object::make_simulation_step(REAL dt, const vector<Object> &objects) {

    advance(dt, acceleration(objects));

}

void step(int steps, vector<Object> objects) {

    times(_, steps) {

        for (auto &object : objects) {

            object.make_simulation_step(1, objects);

        }

    }

}

void make_distinct(vector<Object> &objects) {

    times (i, objects.size()) {

        objects[i].id = i;

    }

}

physics.h

#ifndef PHYSICS_H
#define PHYSICS_H

/*
    N-object gravity simulator (approximative).
*/

#include <unistd.h>
#include <math.h>
#include <string.h>

#define REAL double

typedef struct  Object {

    int  id;
    REAL mass;
    REAL v[3], x[3];

} Object;

void make_simulation_step(Object *, REAL dt, const Object * objects, int count);
void advance             (Object *, REAL dt, REAL * acceleration);

REAL * acceleration(Object *, REAL * result_ptr, const Object * objects, int count);

void step(Object * objects, int count, int steps);

void make_distinct(Object * objects, int count);

#endif

physics.h++

#ifndef PHYSICS_HPP
#define PHYSICS_HPP

/*
    N-object gravity simulator (approximative).
*/

#include <unistd.h>
#include <iostream>
#include <vector>
#include <cmath>
#include <string.h>

typedef double REAL;

using namespace std;

struct Object;

void step(int steps, vector<Object>);

struct Object {

    int  id;
    REAL mass;
    REAL v[3], x[3];

    void make_simulation_step(REAL dt, const vector<Object> &objects);
    void advance             (REAL dt,       vector<REAL>    a);

    vector<REAL> acceleration(const vector<Object> &objects) const;
};

void make_distinct(vector<Object> &objects);

#endif