Throwaway N-body Sim

.gitignore

galaxy
output.*
*.png

README.md

Simple N-body Solver

This is a dead simple implementation using Euler method integration, just for fun and to make some cool videos.

Videos

• http://nullprogram.s3.amazonaws.com/galaxy/attempt-1.webm
• http://nullprogram.s3.amazonaws.com/galaxy/attempt-2.webm
• http://nullprogram.s3.amazonaws.com/galaxy/attempt-3.webm
• http://nullprogram.s3.amazonaws.com/galaxy/v10-z6.webm
• http://nullprogram.s3.amazonaws.com/galaxy/v10-z9.webm

Todo

• Collisions/merging
• Workers (multi-threading)
• SDL display?
• RK4 integration?

galaxy.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <getopt.h>
#include "vec.h"
#include "random.h"

/* Parameters. */
char *progname = "";
unsigned num_stars = 2;
double G = 10.0;
unsigned num_frames = 0;
int quiet = 0;
double galaxy_radius = 100.0;
double initial_velocity = 0.0;

typedef struct {
    vec3 pos, vel, acc;
    double mass;
} star;

star *create_galaxy(unsigned count, double radius)
{
    star *stars = calloc(count, sizeof(star));
    unsigned i;
    vec3 up = {1.0, 0, 0};
    for (i = 0; i < count; i++) {
        stars[i].pos.x =  r_normal() * radius;
        stars[i].pos.y =  r_normal() * radius;
        stars[i].pos.z =  r_normal() * radius;
        stars[i].mass = r_normal() / 10.0 + 1.0;
        if (stars[i].mass <= 0.1) stars[i].mass = 0.1;
        stars[i].vel = norm(cross(stars[i].pos, up));
        double v = r_normal() * initial_velocity / 10.0 + initial_velocity;
        stars[i].vel = vmul(stars[i].vel, v);
    }
    return stars;
}

int main(int argc, char **argv)
{
    int c;
    while ((c = getopt (argc, argv, "n:r:v:G:s:f:q")) != -1) {
        switch (c) {
        case 'n':
            num_stars = atoi(optarg);
            break;
        case 'r':
            galaxy_radius = atof(optarg);
            break;
        case 'v':
            initial_velocity = atof(optarg);
            break;
        case 'G':
            G = atof(optarg);
            break;
        case 's':
            srand(atoi(optarg));
            break;
        case 'f':
            num_frames = atoi(optarg);
            break;
        case 'q':
            quiet = 1;
            break;
        default:
            abort();
        }
    }

    star *stars = create_galaxy(num_stars, galaxy_radius);
    unsigned frame = 0;
    while (num_frames == 0 || frame < num_frames) {
        unsigned i, j;
        /* Update velocity. */
        for (i = 0; i < num_stars; i++) {
            for (j = i + 1; j < num_stars; j++) {
                double r = dist(&stars[i].pos, &stars[j].pos);
                vec3 dir = norm(vdiff(stars[j].pos, stars[i].pos));
                vec3 g = vmul(dir, G / (r * r));
                stars[i].acc = vsum(stars[i].acc, vmul(g,  stars[j].mass));
                stars[j].acc = vsum(stars[j].acc, vmul(g, -stars[i].mass));
            }
        }
        /* Update position. */
        for (i = 0; i < num_stars; i++) {
            stars[i].vel.x += stars[i].acc.x;
            stars[i].vel.y += stars[i].acc.y;
            stars[i].vel.z += stars[i].acc.z;
            stars[i].pos.x += stars[i].vel.x;
            stars[i].pos.y += stars[i].vel.y;
            stars[i].pos.z += stars[i].vel.z;
            stars[i].acc = ZERO;
            printf("%f, %f, %f%s",
                   stars[i].pos.x, stars[i].pos.y, stars[i].pos.z,
                   i == num_stars - 1 ? "\n" : ", ");
        }
        if (!quiet)
            fprintf(stderr, "frame: %u\n", frame);
        frame++;
    }

    return 0;
}

ganim.m

function ganim(fid)
  warning("off", "Octave:broadcast");
  i = 0;
  while ~feof(fid)
    gshow(fscanf(fid, "%f,")');
    title(sprintf("frame %d", i++));
    #drawnow();
    print(sprintf("frame-%08d.png", i), "-dpng", "-r100");
  end
end

gshow.m

function gshow(frame)
  frame = reshape(frame, 3, size(frame, 2) / 3)';
  plot3(frame(:, 1), frame(:, 2), frame(:, 3), ".");
  axis("equal");
  m = mean(frame);
  d = sqrt(sum((frame - m) .^ 2, 2));
  select = median(d);
  s = std(d(d < select)) * 9;
  xlim([-s s] + m(1));
  ylim([-s s] + m(2));
  zlim([-s s] + m(3));
end

Makefile

CFLAGS = -W -Wall -Wextra -O3 -g -ansi
LDLIBS = -lm

galaxy : galaxy.c vec.c vec.h random.c random.h

run : galaxy
	./$< -n 2000 -f 200 | gzip -c > output.txt.gz

clean :
	$(RM) galaxy *.o output.*

process.m

#!/usr/bin/octave -qf

ganim(stdin);

random.c

#include <stdlib.h>
#include <math.h>
#include <inttypes.h>
#include "random.h"

double r_uniform()
{
    union {
        uint16_t i[4];
        uint64_t l;
    } r;
    r.i[0] = rand() & UINT16_MAX;
    r.i[1] = rand() & UINT16_MAX;
    r.i[2] = rand() & UINT16_MAX;
    r.i[3] = rand() & UINT16_MAX;
    return r.l * 1.0 / UINT64_MAX;
}

/* Leftover values from r_normal(). */
int pool_set = 0;
double pool;

/* Box-Muller transform. */
double r_normal()
{
    /* Use leftover value. */
    if (pool_set) {
        pool_set = 0;
        return pool;
    }

    /* Generate a pair of values. */
    double x, y, r;
    do {
        x = 2 * r_uniform() - 1;
        y = 2 * r_uniform() - 1;
        r = x * x + y * y;
    } while (r > 1.0);
    double d = sqrt(-2.0 * log(r) / r);
    pool = x * d;
    pool_set = 1;
    return y * d;
}

random.h

#ifndef RANDOM_H
#define RANDOM_H

/* Generate a uniform random number in [0, 1). */
double r_uniform();

/* Generate a normally-distributed random number with mean 0 and
 * standard deviation 1. */
double r_normal();

#endif

vec.c

#include <stdio.h>
#include <math.h>
#include "vec.h"

vec3 ZERO = {0, 0, 0};

double dist(vec3 *a, vec3 *b)
{
    return sqrt(pow(a->x - b->x, 2) +
                pow(a->y - b->y, 2) +
                pow(a->z - b->z, 2));
}

vec3 vdiv(vec3 v, double d)
{
    v.x /= d;
    v.y /= d;
    v.z /= d;
    return v;
}

vec3 vmul(vec3 v, double d)
{
    v.x *= d;
    v.y *= d;
    v.z *= d;
    return v;
}

vec3 vdiff(vec3 a, vec3 b)
{
    a.x -= b.x;
    a.y -= b.y;
    a.z -= b.z;
    return a;
}

vec3 vsum(vec3 a, vec3 b)
{
    a.x += b.x;
    a.y += b.y;
    a.z += b.z;
    return a;
}

vec3 cross(vec3 a, vec3 b)
{
    vec3 r;
    r.x = a.y * b.z + a.z * b.y;
    r.y = a.z * b.x + a.x * b.z;
    r.z = a.x * b.y + a.y * b.x;
    return r;
}

vec3 norm(vec3 v)
{
    return vdiv(v, dist(&ZERO, &v));
}

void print_vec3(vec3 v)
{
    printf("(%f, %f, %f)\n", v.x, v.y, v.z);
}

vec.h

#ifndef VEC_H
#define VEC_H

typedef struct {
    double x, y, z;
} vec3;

extern vec3 ZERO;

double dist(vec3 *a, vec3 *b);

vec3 vdiv(vec3 v, double d);
vec3 vmul(vec3 v, double d);

vec3 vdiff(vec3 a, vec3 b);
vec3 vsum(vec3 a, vec3 b);

vec3 cross(vec3 a, vec3 b);
vec3 norm(vec3 v);

void print_vec3(vec3 v);

#endif