nbody.ino with correct time measurement

nbody2.ino

/*
 An implementation of the n-body benchmark for Arduino clones.

 Expected output with TIME_STEPS = 50000000:

 -0.169075164
 -0.169059907
 # time: ...

 N-body benchmark (Java7 #2 1 Core):
 http://benchmarksgame.alioth.debian.org/u32/benchmark.php?test=nbody&lang=java
*/

const int TIME_STEPS = 50000000; // command line argument in benchmarksgame
const double DELTA_T = 0.01; // as in benchmarksgame

const double MY_PI = 3.141592653589793;
const double SOLAR_MASS = 4 * MY_PI * MY_PI;
const double DAYS_PER_YEAR = 365.24;

// There are no structs in the Arduino programming language
// http://arduino.cc/en/Reference/HomePage
// Using arrays to keep track of planet variables.
// All arrays have N_PLANETS elements.
const int N_PLANETS = 5;

// Array indices are:
const int JUPITER = 0;
const int SATURN = 1;
const int URANUS = 2;
const int NEPTUNE = 3;
const int SUN = 4;

// Global state of the n-body system
double x[N_PLANETS];
double y[N_PLANETS];
double z[N_PLANETS];
double vx[N_PLANETS];
double vy[N_PLANETS];
double vz[N_PLANETS];
double mass[N_PLANETS];

void init_jupiter() {
    x[JUPITER] = 4.84143144246472090e+00;
    y[JUPITER] = -1.16032004402742839e+00;
    z[JUPITER] = -1.03622044471123109e-01;
    vx[JUPITER] = 1.66007664274403694e-03 * DAYS_PER_YEAR;
    vy[JUPITER] = 7.69901118419740425e-03 * DAYS_PER_YEAR;
    vz[JUPITER] = -6.90460016972063023e-05 * DAYS_PER_YEAR;
    mass[JUPITER] = 9.54791938424326609e-04 * SOLAR_MASS;
}

void init_saturn() {
    x[SATURN] = 8.34336671824457987e+00;
    y[SATURN] = 4.12479856412430479e+00;
    z[SATURN] = -4.03523417114321381e-01;
    vx[SATURN] = -2.76742510726862411e-03 * DAYS_PER_YEAR;
    vy[SATURN] = 4.99852801234917238e-03 * DAYS_PER_YEAR;
    vz[SATURN] = 2.30417297573763929e-05 * DAYS_PER_YEAR;
    mass[SATURN] = 2.85885980666130812e-04 * SOLAR_MASS;
}

void init_uranus() {
    x[URANUS] = 1.28943695621391310e+01;
    y[URANUS] = -1.51111514016986312e+01;
    z[URANUS] = -2.23307578892655734e-01;
    vx[URANUS] = 2.96460137564761618e-03 * DAYS_PER_YEAR;
    vy[URANUS] = 2.37847173959480950e-03 * DAYS_PER_YEAR;
    vz[URANUS] = -2.96589568540237556e-05 * DAYS_PER_YEAR;
    mass[URANUS] = 4.36624404335156298e-05 * SOLAR_MASS;
}

void init_neptune() {
    x[NEPTUNE] = 1.53796971148509165e+01;
    y[NEPTUNE] = -2.59193146099879641e+01;
    z[NEPTUNE] = 1.79258772950371181e-01;
    vx[NEPTUNE] = 2.68067772490389322e-03 * DAYS_PER_YEAR;
    vy[NEPTUNE] = 1.62824170038242295e-03 * DAYS_PER_YEAR;
    vz[NEPTUNE] = -9.51592254519715870e-05 * DAYS_PER_YEAR;
    mass[NEPTUNE] = 5.15138902046611451e-05 * SOLAR_MASS;
}

void init_sun() {
    x[SUN] = 0.0;
    y[SUN] = 0.0;
    z[SUN] = 0.0;
    vx[SUN] = 0.0;
    vy[SUN] = 0.0;
    vz[SUN] = 0.0;
    mass[SUN] = SOLAR_MASS;
}

void init_momentum() {
    double px = 0.0;
    double py = 0.0;
    double pz = 0.0;
    // total momentum
    for (int i=0; i < N_PLANETS; i++) {
        px += vx[i]*mass[i];
        py += vy[i]*mass[i];
        pz += vz[i]*mass[i];
    }
    // offset sun momentum
    vx[SUN] = -px/SOLAR_MASS;
    vy[SUN] = -py/SOLAR_MASS;
    vz[SUN] = -pz/SOLAR_MASS;
}

void init_system() {
    init_jupiter();
    init_saturn();
    init_uranus();
    init_neptune();
    init_sun();
    init_momentum();
}

double energy() {
    double dx, dy, dz, distance;
    double e = 0.0;
    for (int i=0; i < N_PLANETS; i++) {
        e += 0.5 * mass[i] * (vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
        for (int j=i+1; j < N_PLANETS; j++) {
            dx = x[i] - x[j];
            dy = y[i] - y[j];
            dz = z[i] - z[j];
            distance = sqrt(dx*dx + dy*dy + dz*dz);
            e -= (mass[i]*mass[j]) / distance;
        }
    }
    return e;
}

double advance(double dt) {
    // update velocities
    for (int i=0; i < N_PLANETS; i++) {
        for (int j=i+1; j < N_PLANETS; j++) {
            double dx = x[i] - x[j];
            double dy = y[i] - y[j];
            double dz = z[i] - z[j];
            double dist2 = dx*dx + dy*dy + dz*dz;
            double dist = sqrt(dist2);
            double mag = dt / (dist2 * dist);
            double ki = mass[i] * mag;
            double kj = mass[j] * mag;
            vx[i] -= dx * kj;
            vy[i] -= dy * kj;
            vz[i] -= dz * kj;
            vx[j] += dx * ki;
            vy[j] += dy * ki;
            vz[j] += dz * ki;
        }
    }
    // update positions
    for (int i=0; i < N_PLANETS; i++) {
        x[i] += dt * vx[i];
        y[i] += dt * vy[i];
        z[i] += dt * vz[i];
    }
}

void setup() {
    // put your setup code here, to run once:
    Serial.begin(9600);
}

void loop() {
    // put your main code here, to run repeatedly:
    unsigned long start = millis();
    init_system();
    double e_start = energy();
    for (int i=0; i < TIME_STEPS; i++) {
        advance(DELTA_T);
    }
    double e_stop = energy();
    unsigned long stop = millis();
    Serial.println(e_start, 9);
    Serial.println(e_stop, 9);
    Serial.print("# time: ");
    Serial.print(stop-start);
    Serial.print(" ms / ");
    Serial.print(TIME_STEPS);
    Serial.println(" iterations\n");
    delay(1000);
}