Skip to content

Instantly share code, notes, and snippets.

@skeeto skeeto/.gitignore
Created Jul 30, 2012

Embed
What would you like to do?
Throwaway N-body Sim
galaxy
output.*
*.png

Simple N-body Solver

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

Videos

Todo

  • Collisions/merging
  • Workers (multi-threading)
  • SDL display?
  • RK4 integration?
#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;
}
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
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
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.*
#!/usr/bin/octave -qf
ganim(stdin);
#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;
}
#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
#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);
}
#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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.