games/galaxy: add n-body simulator

9 files changed, 1507 insertions, 0 deletions

diff --git a/sys/man/1/galaxy b/sys/man/1/galaxy
new file mode 100644
index 000000000..37fc83377
--- /dev/null
+++ b/sys/man/1/galaxy
@@ -0,0 +1,231 @@
+.TH GALAXY 1
+.SH NAME
+galaxy, mkgalaxy \- galactic n-body simulator
+.SH SYNOPSIS
+.B games/galaxy
+[
+.I options
+] [
+.B -i
+] [
+.I file
+]
+.br
+.B games/mkgalaxy
+[
+.I options
+] [
+.B -i
+] [
+.B -f
+.I file
+]
+.I size
+.SH DESCRIPTION
+.I Galaxy
+is an n-body simulator that uses a Barnes-Hut quad-tree
+to calculate gravitational interactions.
+Typical usage is to read a galaxy file (see
+.IR galaxy (6))
+from standard input
+using the
+.B -i
+command-line option or from a
+.I file
+using the
+.B -f
+option. If no file is read then the simulator starts with an empty
+universe.
+.SS Mouse commands
+.PP
+New planetary bodies can be created with mouse button 1.
+Holding button 1 will reposition the body.
+Holding a button 1-2 chord changes the mass/size
+of the body. Holding a button 1-3 chord
+changes the initial velocity of the body. Releasing button 1
+restarts the simulator with the new body in motion. When new
+bodies are created, the simulator maintains the Galilean (inertial)
+reference frame.
+.PP
+Mouse button 3 opens a menu with the following options:
+.TP
+.B move
+Change the visible region of the simulation
+by holding button 1 and dragging. Any other mouse
+button restarts the simulation.
+.TP
+.B zoom
+Prompts for a floating point value to change the scale of the
+simulation. E.g. a value of 2 will halve the scale (zoom in)
+and a value of 0.5 will double the scale (zoom out).
+.TP
+.B speed
+Prompts for a floating point value to change the speed of
+the simulation. E.g. a value of 2 will double the speed
+of the simulation and a value of 0.5 will
+halve the speed. Accuracy is sacrificed for greater speed.
+.TP
+.B gravity
+Prompts for a floating point value to change the gravitational
+constant. E.g. a value of 2 will double the force exerted by
+gravity and a value of 0.5 will halve it.
+.TP
+.B save
+Prompts for a file name to save the current galaxy as a
+.IR galaxy (6)
+file.
+.TP
+.B load
+Prompts for a file name to load the galaxy from the
+.IR galaxy (6)
+file.
+.TP
+.B exit
+Exits the simulator.
+.SS Keyboard commands
+The following keys are recognized as commands:
+.TP
+.B a
+Show accelerations as vectors.
+.TP
+.B v
+Show velocities as vectors.
+.TP
+.B s
+Show statistics such as the number of bodies being
+simulated, the maximum depth of the quad-tree, and the
+average number of calculations made per body.
+.TP
+.B q
+Exit the simulator.
+.TP
+.B space
+Pause and unpause the simulator.
+.TP
+.B del
+Exit the simulator.
+.SS Command-line options
+Certain aspects of the galaxy simulator are controlled by
+the following options:
+.TP
+.BI -t " throttle"
+Causes the process that calculates forces to relinquish
+the processor for
+.I throttle
+milliseconds after each calculation.
+.TP
+.BI -G " gravity"
+Sets the gravitational constant to
+.I gravity.
+The default value is 1.
+.TP
+.BI -ε " softening"
+Sets the
+.I softening
+factor to prevent gravitational singularities during
+collisions or near-collisions. The default value is 500.
+.TP
+.BI -f " file"
+Reads the galaxy file
+.I file
+(see
+.IR galaxy (6)).
+.TP
+.B -i
+Reads a galaxy file from standard input.
+.SS Mkgalaxy
+.I Mkgalaxy
+is a utility to create galaxies for simulation.
+Galaxies can be assembled incrementally by reading an
+existing galaxy file from standard input with the
+.B -i
+command-line option or from a
+.I file
+with the
+.B -f
+option. Mkgalaxy then writes to standard output a
+.IR galaxy (6)
+file with a galaxy of the given
+.I size
+together with the previously read galaxy.
+Galaxies generated by mkgalaxy have characteristics
+determined by the following options:
+.TP
+.BI -d " distance"
+.I Distance
+determines the spacing between bodies.
+The default value is 100.
+.TP
+.BI -s " size"
+Bodies have the given
+.IR size .
+The default value is 25.
+.TP
+.BI -v " velocity"
+Bodies have the given
+.I velocity
+in a random direction.
+The default value is 0.
+.TP
+.BI -av " angular velocity"
+Bodies have the given
+.I "angular velocity"
+relative to the center of mass of the new galaxy being generated.
+The default value is 0.
+.TP
+.BI -gv " x,y"
+The entire galaxy being generated is given the directional velocity determined
+by the vector
+.RI ( x,y ).
+The default value is (0, 0).
+.TP
+.BI -o " x,y"
+The entire galaxy being generated is offset by the vector
+.RI ( x,y ).
+The default value is (0, 0).
+.TP
+.B -sq
+The galaxy being generated is a square. Without this option, the galaxy
+will be circular.
+.PP
+The arguments to the
+.BR -d ,
+.BR -s ,
+.BR -v ,
+and
+.B -av
+arguments have the form
+.B s
+or
+.B s±r
+where s and r are double-precision floating point numbers.
+.B S
+is the base value and
+.B r
+if given determines a range in which the value will vary randomly
+from the base.
+.SH EXAMPLES
+Two rotating circles destroy each other:
+.IP
+.EX
+games/mkgalaxy -av 100 -d 60±50 -v 10 2000 |
+games/mkgalaxy -i -av -70 -d 80±50 -v 10 -o 6000,2000 -gv -80,40 3000 |
+games/galaxy -i
+.PP
+Cool patterns made by a square galaxy:
+.IP
+.EX
+games/mkgalaxy -sq -av 20 5000 | games/galaxy -i
+.SH SOURCE
+.B /sys/src/games/galaxy
+.SH SEE ALSO
+J. Barnes & P. Hut (December 1986). "A hierarchical O(N log N) force-calculation algorithm".
+.IR Nature .
+324 (4): 446–449.
+.PP
+.IR galaxy (6)
+.SH HISTORY
+.I Galaxy
+and
+.I mkgalaxy
+first appeared in 9front (Feb, 2017).
diff --git a/sys/man/6/galaxy b/sys/man/6/galaxy
new file mode 100644
index 000000000..b73649847
--- /dev/null
+++ b/sys/man/6/galaxy
@@ -0,0 +1,41 @@
+.TH GALAXY 6
+.SH NAME
+galaxy \- representations of n-body simulations
+.SH DESCRIPTION
+Files of this format are interpreted by
+.IR galaxy (1)
+as describing the inital condition of n-body simulations
+or the saved state of simulation in progress.
+A galaxy file is a UTF stream of instruction lines. The
+instruction is given by the first space delimited word. The
+following instructions are accepted.
+.TP
+.B MKBODY
+The rest of the line must contain 5 white space delimited
+double-precision floating point numbers. They represent a
+body's x coordinate, y coordinate, x velocity component,
+y velocity component, and size respectively.
+.TP
+.B ORIG
+The rest of the line must contain 2 white space delimited
+double-precision floating point numbers. They represent
+the current location of the origin with respect to the
+view window of
+.IR galaxy (1).
+.TP
+.B DT
+The rest of the line must contain a double-precision
+floating point number which determines the time-scale
+of the simulation.
+.TP
+.B SCALE
+The rest of the line must contain a double-precision
+floating point number which determines the scale
+of the view of the simulation.
+.TP
+.B GRAV
+The rest of the line must contain a double-precision
+floating point number which determines the gravitational
+constant of the simulation.
+.SH SEE ALSO
+.IR galaxy (1)
diff --git a/sys/src/games/galaxy/body.c b/sys/src/games/galaxy/body.c
new file mode 100644
index 000000000..ea610ad56
--- /dev/null
+++ b/sys/src/games/galaxy/body.c
@@ -0,0 +1,210 @@
+#include <u.h>
+#include <libc.h>
+#include <draw.h>
+#include <bio.h>
+#include "galaxy.h"
+
+void
+glxyinit(void)
+{
+
+	glxy.a = calloc(5, sizeof(Body));
+	if(glxy.a == nil)
+		sysfatal("could not allocate glxy: %r");
+	glxy.l = 0;
+	glxy.max = 5;
+}
+
+Body*
+body(void)
+{
+	Body *b;
+
+	if(glxy.l == glxy.max) {
+		glxy.max *= 2;
+		glxy.a = realloc(glxy.a, sizeof(Body) * glxy.max);
+		if(glxy.a == nil)
+			sysfatal("could not realloc glxy: %r");
+	}
+	b = glxy.a + glxy.l++;
+	*b = ZB;
+	return b;
+}
+
+void
+drawbody(Body *b)
+{
+	Point pos, v;
+	int s;
+
+	pos.x = b->x / scale + orig.x;
+	pos.y = b->y / scale + orig.y;
+	s = b->size/scale;
+	fillellipse(screen, pos, s, s, b->col, ZP);
+	v.x = b->v.x/scale*10;
+	v.y = b->v.y/scale*10;
+	if(v.x != 0 || v.y != 0)
+		line(screen, pos, addpt(pos, v), Enddisc, Endarrow, 0, b->col, ZP);
+	flushimage(display, 1);
+}
+
+Vector
+center(void)
+{
+	Body *b;
+	Vector gc, gcv;
+	double mass;
+
+	if(glxy.l == 0)
+		return (Vector){0, 0};
+
+	gc.x = gc.y = gcv.x = gcv.y = mass = 0;
+	for(b = glxy.a; b < glxy.a+glxy.l; b++) {
+		gc.x += b->x * b->mass;
+		gc.y += b->y * b->mass;
+		gcv.x += b->v.x * b->mass;
+		gcv.y += b->v.y * b->mass;
+		mass += b->mass;
+	}
+	gc.x /= mass;
+	gc.y /= mass;
+	gcv.x /= mass;
+	gcv.y /= mass;
+	for(b = glxy.a; b < glxy.a+glxy.l; b++) {
+		b->x -= gc.x;
+		b->y -= gc.y;
+		b->v.x -= gcv.x;
+		b->v.y -= gcv.y;
+	}
+	return gc;
+}
+
+int
+Bfmt(Fmt *f)
+{
+	Body *b;
+	int r;
+
+	b = va_arg(f->args, Body*);
+
+	r = fmtprint(f, "MKBODY %g %g ", b->x, b->y);
+	if(r < 0)
+		return -1;
+
+	r = fmtprint(f, "%g %g ", b->v.x, b->v.y);
+	if(r < 0)
+		return -1;
+
+	return fmtprint(f, "%g", b->size);
+}
+
+enum {
+	MKBODY,
+	ORIG,
+	DT,
+	SCALE,
+	GRAV,
+	NOCMD,
+};
+
+int
+getcmd(char *l)
+{
+	static char *cmds[] = {
+		[MKBODY]	"MKBODY",
+		[ORIG]	"ORIG",
+		[DT]	"DT",
+		[SCALE]	"SCALE",
+		[GRAV]	"GRAV",
+	};
+	int cmd;
+
+	for(cmd = 0; cmd < nelem(cmds); cmd++) {
+		if(strcmp(l, cmds[cmd]) == 0)
+			return cmd;
+	}
+	sysfatal("getcmd: no such command %s", l);
+	return NOCMD;
+}
+
+void
+readglxy(int fd)
+{
+	static Biobuf bin;
+	char *line;
+	double f;
+	int cmd, len;
+	Body *b;
+
+	glxy.l = 0;
+	Binit(&bin, fd, OREAD);
+	for(;;) {
+		line = Brdline(&bin, ' ');
+		len = Blinelen(&bin);
+		if(line == nil) {
+			if(len == 0)
+				break;
+			sysfatal("load: malformed command");
+		}
+
+		line[len-1] = '\0';
+		cmd = getcmd(line);
+
+		line = Brdline(&bin, '\n');
+		if(line == nil) {
+			if(len == 0)
+				sysfatal("load: malformed command");
+			sysfatal("load: read error: %r");
+		}
+		len = Blinelen(&bin);
+		line[len-1] = '\0';
+
+		switch(cmd) {
+		case MKBODY:
+			b = body();
+			b->x = strtod(line, &line);
+			b->y = strtod(line, &line);
+			b->v.x = strtod(line, &line);
+			b->v.y = strtod(line, &line);
+			b->size = strtod(line, nil);
+			b->mass = b->size*b->size*b->size;
+			b->col = randcol();
+			CHECKLIM(b, f);
+			break;
+		case ORIG:
+			orig.x = strtol(line, &line, 10);
+			orig.y = strtol(line, nil, 10);
+			break;
+		case DT:
+			dt = strtod(line, nil);
+			dt² = dt*dt;
+			break;
+		case SCALE:
+			scale = strtod(line, nil);
+			break;
+		case GRAV:
+			G = strtod(line, nil);
+			break;
+		}
+	}
+	Bterm(&bin);
+}
+
+void
+writeglxy(int fd)
+{
+	static Biobuf bout;
+	Body *b;
+
+	Binit(&bout, fd, OWRITE);
+	
+	Bprint(&bout, "ORIG %d %d\n", orig.x, orig.y);
+	Bprint(&bout, "SCALE %g\n", scale);
+	Bprint(&bout, "DT %g\n", dt);
+	Bprint(&bout, "GRAV %g\n", G);
+
+	for(b = glxy.a; b < glxy.a + glxy.l; b++)
+		Bprint(&bout, "%B\n", b);
+
+	Bterm(&bout);
+}
diff --git a/sys/src/games/galaxy/galaxy.c b/sys/src/games/galaxy/galaxy.c
new file mode 100644
index 000000000..3b5d5d768
--- /dev/null
+++ b/sys/src/games/galaxy/galaxy.c
@@ -0,0 +1,616 @@
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#include <draw.h>
+#include <thread.h>
+#include <mouse.h>
+#include <cursor.h>
+#include <keyboard.h>
+#include "galaxy.h"
+
+Cursor crosscursor = {
+	{-7, -7},
+	{0x03, 0xC0, 0x03, 0xC0, 0x03, 0xC0, 0x03, 0xC0,
+	 0x03, 0xC0, 0x03, 0xC0, 0xFF, 0xFF, 0xFF, 0xFF,
+	 0xFF, 0xFF, 0xFF, 0xFF, 0x03, 0xC0, 0x03, 0xC0,
+	 0x03, 0xC0, 0x03, 0xC0, 0x03, 0xC0, 0x03, 0xC0, },
+	{0x00, 0x00, 0x01, 0x80, 0x01, 0x80, 0x01, 0x80,
+	 0x01, 0x80, 0x01, 0x80, 0x01, 0x80, 0x7F, 0xFE,
+	 0x7F, 0xFE, 0x01, 0x80, 0x01, 0x80, 0x01, 0x80,
+	 0x01, 0x80, 0x01, 0x80, 0x01, 0x80, 0x00, 0x00, }
+};
+
+Cursor pausecursor={
+	0, 0,
+	0x01, 0x80, 0x03, 0xC0, 0x07, 0xE0, 0x07, 0xe0,
+	0x07, 0xe0, 0x07, 0xe0, 0x03, 0xc0, 0x0F, 0xF0,
+	0x1F, 0xF8, 0x1F, 0xF8, 0x1F, 0xF8, 0x1F, 0xF8,
+	0x0F, 0xF0, 0x1F, 0xF8, 0x3F, 0xFC, 0x3F, 0xFC,
+
+	0x01, 0x80, 0x03, 0xC0, 0x07, 0xE0, 0x04, 0x20,
+	0x04, 0x20, 0x06, 0x60, 0x02, 0x40, 0x0C, 0x30,
+	0x10, 0x08, 0x14, 0x08, 0x14, 0x28, 0x12, 0x28,
+	0x0A, 0x50, 0x16, 0x68, 0x20, 0x04, 0x3F, 0xFC,
+};
+
+enum {
+	STK = 8192,
+	MOVE = 0,
+	ZOOM,
+	SPEED,
+	GRAV,
+	SAVE,
+	LOAD,
+	EXIT,
+	MEND,
+};
+
+Cursor *cursor;
+Mousectl *mc;
+Keyboardctl kc;
+double
+	G = 1,
+	θ = 1,
+	scale = 30,
+	ε = 500,
+	dt = .2,
+	LIM = 10,
+	dt²;
+char *file;
+int showv, showa, throttle, paused;
+
+char *menustr[] = {
+	[SAVE]	"save",
+	[LOAD]	"load",
+	[ZOOM]	"zoom",
+	[SPEED]	"speed",
+	[GRAV]	"gravity",
+	[MOVE]	"move",
+	[EXIT]	"exit",
+	[MEND]	nil
+};
+Menu menu = {
+	.item menustr
+};
+
+Image*
+randcol(void)
+{
+	static struct {
+		ulong c;
+		Image *i;
+	} cols[] = {
+		DWhite, nil,
+		DRed, nil,
+		DGreen, nil,
+		DCyan, nil,
+		DMagenta, nil,
+		DYellow, nil,
+		DPaleyellow, nil,
+		DDarkyellow, nil,
+		DDarkgreen, nil,
+		DPalegreen, nil,
+		DPalebluegreen, nil,
+		DPaleblue, nil,
+		DPalegreygreen, nil,
+		DYellowgreen, nil,
+		DGreyblue, nil,
+		DPalegreyblue, nil,
+	};
+	int r;
+
+	r = nrand(nelem(cols));
+	if(cols[r].i == nil)
+		cols[r].i = allocimage(display, Rect(0,0,1,1), screen->chan, 1, cols[r].c);
+	return cols[r].i;
+}
+
+void
+pause(int p, int pri)
+{
+	static int paused, ppri;
+
+	switch(p) {
+	default:
+		sysfatal("invalid pause value %d:", p);
+		break;
+	case 0:
+		if(pri > ppri)
+			ppri = pri;
+		if(paused)
+			break;
+		paused = 1;
+		qlock(&glxy);
+		break;
+	case 1:
+		if(!paused || pri < ppri)
+			break;
+		paused = ppri = 0;
+		qunlock(&glxy);
+		break;
+	}
+}
+
+void
+drawstats(void)
+{
+	Point p;
+	static char buf[1024];
+
+	snprint(buf, sizeof(buf), "Number of bodies: %d", glxy.l);
+	p = addpt(screen->r.min, (Point){5, 3});
+	string(screen, p, display->white, ZP, font, buf);
+
+	snprint(buf, sizeof(buf), "Avg. calculations per body: %g", avgcalcs);
+	p = addpt(p, (Point){0, font->height});
+	string(screen, p, display->white, ZP, font, buf);
+
+	snprint(buf, sizeof(buf), "Max depth of quad tree: %d", quaddepth);
+	p = addpt(p, (Point){0, font->height});
+	string(screen, p, display->white, ZP, font, buf);
+}
+
+void
+drawglxy(void)
+{
+	Point pos, va;
+	Body *b;
+	int s;
+
+	draw(screen, screen->r, display->black, 0, ZP);
+	for(b = glxy.a; b < glxy.a + glxy.l; b++) {
+		pos.x = b->x / scale + orig.x;
+		pos.y = b->y / scale + orig.y;
+		s = b->size/scale;
+		fillellipse(screen, pos, s, s, b->col, ZP);
+		if(showv) {
+			va.x = b->v.x/scale;
+			va.y = b->v.y/scale;
+			if(va.x != 0 || va.y != 0)
+				line(screen, pos, addpt(pos, va), Enddisc, Endarrow, 0, b->col, ZP);
+		}
+		if(showa) {
+			va.x = b->a.x/scale*50;
+			va.y = b->a.y/scale*50;
+			if(va.x != 0 || va.y != 0)
+				line(screen, pos, addpt(pos, va), Enddisc, Endarrow, 0, b->col, ZP);
+		}
+	}
+	STATS(drawstats();)
+	flushimage(display, 1);
+}
+
+void
+setsize(Body *b)
+{
+	Point pos, d;
+	double h;
+
+	pos.x = b->x / scale + orig.x;
+	pos.y = b->y / scale + orig.y;
+	d = subpt(mc->xy, pos);
+	h = hypot(d.x, d.y);
+	b->size = h == 0 ? scale : h*scale;
+	b->mass = b->size*b->size*b->size;
+}
+
+void
+setvel(Body *b)
+{
+	Point pos, d;
+
+	pos.x = b->x / scale + orig.x;
+	pos.y = b->y / scale + orig.y;
+	d = subpt(mc->xy, pos);
+	b->v.x = (double)d.x*scale/10;
+	b->v.y = (double)d.y*scale/10;
+}
+
+void
+setpos(Body *b)
+{
+	b->x = (mc->xy.x - orig.x) * scale;
+	b->y = (mc->xy.y - orig.y) * scale;
+}
+
+void
+dosize(Body *b)
+{
+	Point p;
+
+	p = mc->xy;
+	for(;;) {
+		setsize(b);
+		drawglxy();
+		drawbody(b);
+		readmouse(mc);
+		if(mc->buttons != 3)
+			break;
+	}
+	moveto(mc, p);
+}
+
+void
+dovel(Body *b)
+{
+	Point p;
+	p = mc->xy;
+	for(;;) {
+		setvel(b);
+		drawglxy();
+		drawbody(b);
+		readmouse(mc);
+		if(mc->buttons != 5)
+			break;
+	}
+	moveto(mc, p);
+}
+
+void
+dobody(void)
+{
+	Vector gc;
+	double f;
+	Body *b;
+
+	pause(0, 0);
+	b = body();
+	setpos(b);
+	setvel(b);
+	setsize(b);
+	b->col = randcol();
+	for(;;) {
+		drawglxy();
+		drawbody(b);
+		readmouse(mc);
+		if(!(mc->buttons & 1))
+			break;
+		if(mc->buttons == 3)
+			dosize(b);
+		else if(mc->buttons == 5)
+			dovel(b);
+		else
+			setpos(b);
+	}
+
+	CHECKLIM(b, f);
+
+	gc = center();
+	orig.x += gc.x / scale;
+	orig.y += gc.y / scale;
+
+	pause(1, 0);
+}
+
+char*
+getinput(char *info, char *sug)
+{
+	static char buf[1024];
+	static Channel *rchan;
+	char *input;
+	int r;
+
+	if(rchan == nil)
+		rchan = chancreate(sizeof(Rune), 20);
+
+	if(sug != nil)
+		strecpy(buf, buf+1024, sug);
+	else
+		buf[0] = '\0';
+
+	kc.c = rchan;
+	r = enter(info, buf, sizeof(buf), mc, &kc, nil);
+	kc.c = nil;
+	if(r < 0)
+		sysfatal("save: could not get filename: %r");
+
+	input = strdup(buf);
+	if(input == nil)
+		sysfatal("getinput: could not save input: %r");
+	return input;
+}
+
+void
+move(void)
+{
+	Point od;
+	setcursor(mc, &crosscursor);
+	for(;;) {
+		for(;;) {
+			readmouse(mc);
+			if(mc->buttons & 1)
+				break;
+			if(mc->buttons & 4) {
+				setcursor(mc, cursor);
+				return;
+			}
+		}
+		od = subpt(orig, mc->xy);
+		for(;;) {
+			readmouse(mc);
+			if(!(mc->buttons & 1))
+				break;
+			orig = addpt(od, mc->xy);
+			drawglxy();
+		}
+	}
+}
+
+void
+load(int fd)
+{
+	orig = divpt(subpt(screen->r.max, screen->r.min), 2);
+	orig = addpt(orig, screen->r.min);
+	readglxy(fd);
+	center();
+}
+
+void
+domenu(void)
+{
+	int fd;
+	char *s;
+	double z;
+
+	pause(0, 0);
+	switch(menuhit(3, mc, &menu, nil)) {
+	case SAVE:
+		s = getinput("Enter file:", file);
+		if(s == nil || *s == '\0')
+			break;
+		free(file);
+		file = s;
+		fd = create(file, OWRITE, 0666);
+		if(fd < 0)
+			sysfatal("domenu: could not create file %s: %r", file);
+		writeglxy(fd);
+		close(fd);
+		break;
+	case LOAD:
+		s = getinput("Enter file:", file);
+		if(s == nil || *s == '\0')
+			break;
+		free(file);
+		file = s;
+		fd = open(file, OREAD);
+		if(fd < 0)
+			sysfatal("domenu: could not open file %s: %r", file);
+		load(fd);
+		close(fd);
+		break;
+	case ZOOM:
+		s = getinput("Zoom multiplier:", nil);
+		if(s == nil || *s == '\0')
+			break;
+		z = strtod(s, nil);
+		free(s);
+		if(z <= 0)
+			break;
+		scale /= z;
+		break;
+	case SPEED:
+		s = getinput("Speed multiplier:", nil);
+		if(s == nil || *s == '\0')
+			break;
+		z = strtod(s, nil);
+		free(s);
+		if(z <= 0)
+			break;
+		dt *= z;
+		dt² = dt*dt;
+		break;
+	case GRAV:
+		s = getinput("Gravity multiplier:", nil);
+		if(s == nil || *s == '\0')
+			break;
+		z = strtod(s, nil);
+		free(s);
+		if(z <= 0)
+			break;
+		G *= z;
+		break;
+	case MOVE:
+		move();
+		break;
+	case EXIT:
+		threadexitsall(nil);
+		break;
+	}
+	drawglxy();
+	pause(1, 0);
+}
+
+void
+mousethread(void*)
+{
+	threadsetname("mouse");
+	for(;;) {
+		readmouse(mc);
+		switch(mc->buttons) {
+		case 1:
+			dobody();
+			break;
+		case 4:
+			domenu();
+			break;
+		}
+	}
+}
+
+void
+resizethread(void*)
+{
+	threadsetname("resize");
+	for(;;) {
+		recv(mc->resizec, nil);
+		pause(0, 0);
+		if(getwindow(display, Refnone) < 0)
+			sysfatal("resize failed: %r");
+		drawglxy();
+		pause(1, 0);
+	}
+}
+
+void
+kbdthread(void*)
+{
+	Keyboardctl *realkc;
+	Rune r;
+
+	threadsetname("keyboard");
+	if(realkc = initkeyboard(nil), realkc == nil)
+		sysfatal("kbdthread: could not initkeyboard: %r");
+
+	for(;;) {
+		recv(realkc->c, &r);
+		if(r == Kdel) {
+			closedisplay(display);
+			threadexitsall(nil);
+		}
+		if(kc.c != nil)
+			send(kc.c, &r);
+		else switch(r) {
+		case 'q':
+			closedisplay(display);
+			threadexitsall(nil);
+			break;
+		case 's':
+			stats ^= 1;
+			break;
+		case 'v':
+			showv ^= 1;
+			break;
+		case 'a':
+			showa ^= 1;
+			break;
+		case ' ':
+			paused ^= 1;
+			if(paused) {
+				cursor = &pausecursor;
+				pause(0, 1);
+			} else {
+				cursor = nil;
+				pause(1, 1);
+			}
+			setcursor(mc, cursor);
+		}
+	}
+}
+
+/* verlet barnes-hut */
+void
+simulate(void*)
+{
+	Body *b;
+	double f;
+
+	threadsetname("simulate");
+
+	for(;;) {
+		qlock(&glxy);
+
+		if(throttle)
+			sleep(throttle);
+
+		drawglxy();
+
+Again:
+		space.t = EMPTY;
+		quads.l = 0;
+		STATS(quaddepth = 0;)
+		for(b = glxy.a; b < glxy.a + glxy.l; b++) {
+			if(quadins(b, LIM) == -1) {
+				growquads();
+				goto Again;
+			}
+		}
+
+		STATS(avgcalcs = 0;)
+		for(b = glxy.a; b < glxy.a + glxy.l; b++) {
+			b->a.x = b->newa.x;
+			b->a.y = b->newa.y;
+			b->newa.x = b->newa.y = 0;
+			STATS(calcs = 0;)
+			quadcalc(space, b, LIM);
+			STATS(avgcalcs += calcs;)
+		}
+		STATS(avgcalcs /= glxy.l;)
+
+		for(b = glxy.a; b < glxy.a + glxy.l; b++) {
+			b->x += dt*b->v.x + dt²*b->a.x/2;
+			b->y += dt*b->v.y + dt²*b->a.y/2;
+			b->v.x += dt*(b->a.x + b->newa.x)/2;
+			b->v.y += dt*(b->a.y + b->newa.y)/2;
+			CHECKLIM(b, f);
+		}
+
+		qunlock(&glxy);
+	}
+}
+
+void
+usage(void)
+{
+	fprint(2, "Usage: %s [-t throttle] [-G gravity] [-ε smooth] [-i] [file]\n", argv0);
+	threadexitsall("usage");
+}
+
+void
+threadmain(int argc, char **argv)
+{
+	int doload;
+
+	doload = 0;
+	ARGBEGIN {
+	default:
+		usage();
+		break;
+	case 't':
+		throttle = strtol(EARGF(usage()), nil, 0);
+		break;
+	case 'G':
+		G = strtod(EARGF(usage()), nil);
+		break;
+	case L'ε':
+		ε = strtod(EARGF(usage()), nil);
+		break;
+	case 'i':
+		doload++;
+		break;
+	} ARGEND
+
+	if(argc > 1)
+		usage();
+
+	fmtinstall('B', Bfmt);
+
+	if(argc == 1) {
+		if(doload++)
+			usage();
+		file = strdup(argv[0]);
+		if(file == nil)
+			sysfatal("threadmain: could not save file name: %r");
+		close(0);
+		if(open(file, OREAD) != 0)
+			sysfatal("threadmain: could not open file: %r");
+	}
+
+	if(initdraw(nil, nil, "Galaxy") < 0)
+		sysfatal("initdraw failed: %r");
+	if(mc = initmouse(nil, screen), mc == nil)
+		sysfatal("initmouse failed: %r");
+
+	dt² = dt*dt;
+	orig = divpt(subpt(screen->r.max, screen->r.min), 2);
+	orig = addpt(orig, screen->r.min);
+	glxyinit();
+	quadsinit();
+	if(doload)
+		load(0);
+	close(0);
+	threadcreate(mousethread, nil, STK);
+	threadcreate(resizethread, nil, STK);
+	threadcreate(kbdthread, nil, STK);
+	proccreate(simulate, nil, STK);
+	threadexits(nil);
+}
diff --git a/sys/src/games/galaxy/galaxy.h b/sys/src/games/galaxy/galaxy.h
new file mode 100644
index 000000000..0ba6f21cb
--- /dev/null
+++ b/sys/src/games/galaxy/galaxy.h
@@ -0,0 +1,83 @@
+typedef struct QB QB;
+typedef struct Body Body;
+typedef struct Quad Quad;
+typedef struct Vector Vector;
+
+struct QB {
+	union {
+		Quad *q;
+		Body *b;
+	};
+	int t;
+};
+
+struct Vector {
+	double x, y;
+};
+
+struct Body {
+	Vector;
+	Vector v, a, newa;
+	double size, mass;
+	Image *col;
+};
+
+struct Quad {
+	Vector;
+	QB c[4];
+	double mass;
+};
+
+#pragma varargck type "B" Body*
+
+struct {
+	QLock;
+	Body *a;
+	int l, max;
+} glxy;
+
+struct {
+	Quad *a;
+	int l, max;
+} quads;
+
+#define π2 6.28318530718e0
+
+enum {
+	EMPTY,
+	QUAD,
+	BODY,
+};
+
+Point orig;
+double G, θ, scale, ε, LIM, dt, dt²;
+Body ZB;
+QB space;
+
+Image *randcol(void);
+
+Body *body(void);
+void drawbody(Body*);
+Vector center(void);
+void glxyinit(void);
+void readglxy(int);
+void writeglxy(int);
+int Bfmt(Fmt*);
+
+void quadcalc(QB, Body*, double);
+int quadins(Body*, double);
+void growquads(void);
+void quadsinit(void);
+
+int stats;
+int quaddepth;
+int calcs;
+double avgcalcs;
+
+#define STATS(e) if(stats) {e}
+
+#define CHECKLIM(b, f) \
+	if(((f) = fabs((b)->x)) > LIM/2)	\
+		LIM = (f)*2;	\
+	if(((f) = fabs((b)->y)) > LIM/2)	\
+		LIM = (f)*2
diff --git a/sys/src/games/galaxy/mkfile b/sys/src/games/galaxy/mkfile
new file mode 100644
index 000000000..91262e8bf
--- /dev/null
+++ b/sys/src/games/galaxy/mkfile
@@ -0,0 +1,14 @@
+</$objtype/mkfile
+
+TARG=galaxy mkgalaxy
+BIN=/$objtype/bin/games
+OGALAXY=galaxy.$O quad.$O body.$O
+OMKGALAXY=mkgalaxy.$O body.$O
+
+</sys/src/cmd/mkmany
+
+$O.galaxy:	$OGALAXY
+	$LD $LDFLAGS -o $target $prereq
+
+$O.mkgalaxy: $OMKGALAXY
+	$LD $LDFLAGS -o $target $prereq
diff --git a/sys/src/games/galaxy/mkgalaxy.c b/sys/src/games/galaxy/mkgalaxy.c
new file mode 100644
index 000000000..22aa32a2f
--- /dev/null
+++ b/sys/src/games/galaxy/mkgalaxy.c
@@ -0,0 +1,171 @@
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#include <draw.h>
+#include "galaxy.h"
+
+Vector o, gv;
+double
+	d = 100, drand,
+	sz = 25, szrand,
+	v, vrand,
+	av, avrand;
+int new, c = 1;
+
+void quadcalc(QB, Body*, double){}
+Image *randcol(void){ return nil; }
+
+void
+usage(void)
+{
+	fprint(2, "Usage: %s [-d dist[±r]]\n\t[-s size[±r]] [-v vel[±r]]\n\t[-av angvel[±r]] [-gv xdir,ydir]\n\t[-o xoff,yoff] [-f file]\n\t[-sq] [-i] size\n", argv0);
+	exits("usage");
+}
+
+Vector
+polar(double ang, double mag)
+{
+	Vector v;
+
+	v.x = cos(ang)*mag;
+	v.y = sin(ang)*mag;
+	return v;
+}
+
+Vector
+getvec(char *str)
+{
+	Vector v;
+
+	v.x = strtod(str, &str);
+	if(*str != ',')
+		usage();
+	v.y = strtod(str+1, nil);
+	return v;
+}
+
+double
+getvals(char *str, double *rand)
+{
+	Rune r;
+	double val;
+	int i;
+
+	val = strtod(str, &str);
+	i = chartorune(&r, str);
+	if(r == L'±')
+		*rand = strtod(str+i, nil);
+	else
+		*rand = 0;
+	return val;
+}
+
+#define RAND(r)	((r)*(frand()*2 - 1))
+
+void
+mkbodies(double lim)
+{
+	Body *b;
+	Vector p;
+	double x, y;
+
+	for(x = -lim/2; x < lim/2; x += d)
+	for(y = -lim/2; y < lim/2; y += d) {
+		p.x = x + RAND(drand);
+		p.y = y + RAND(drand);
+		if(c)
+		if(hypot(p.x, p.y) > lim/2)
+			continue;
+		b = body();
+		b->Vector = p;
+		b->v = polar(frand()*π2, v+RAND(vrand));
+		b->v.x += gv.x - p.y*(av + RAND(avrand))/1000;
+		b->v.y += gv.y + p.x*(av + RAND(avrand))/1000;
+		b->size = sz + RAND(szrand);
+	}
+}
+
+void
+main(int argc, char **argv)
+{
+	static Biobuf bout;
+	Body *b;
+	double lim;
+	int fd;
+	char *a;
+
+	srand(truerand());
+	fmtinstall('B', Bfmt);
+	glxyinit();
+
+	ARGBEGIN {
+	case 'f':
+		fd = open(EARGF(usage()), OREAD);
+		if(fd < 0)
+			sysfatal("Could not open file %s: %r", *argv);
+		readglxy(fd);
+		close(fd);
+		break;
+	case 'i':
+		readglxy(0);
+		break;
+	case 's':
+		a = EARGF(usage());
+		switch(a[0]) {
+		case 'q':
+			if(a[1] != '\0')
+				usage();
+			c = 0;
+			break;
+		default:
+			sz = getvals(a, &szrand);
+			break;
+		}
+		break;
+	case 'a':
+		a = EARGF(usage());
+		if(a[0] != 'v' || a[1] != '\0')
+			usage();
+		argc--;
+		argv++;
+		av = getvals(*argv, &avrand);
+		break;
+	case 'g':
+		a = EARGF(usage());
+		if(a[0] != 'v' || a[1] != '\0')
+			usage();
+		argc--;
+		argv++;
+		gv = getvec(*argv);
+		break;
+	case 'v':
+		v = getvals(EARGF(usage()), &vrand);
+		break;
+	case 'o':
+		o = getvec(EARGF(usage()));
+		break;
+	case 'd':
+		d = getvals(EARGF(usage()), &drand);
+		break;
+	} ARGEND
+
+	if(argc != 1)
+		usage();
+
+	new = glxy.l;
+	lim = strtod(*argv, nil);
+	mkbodies(lim);
+
+	Binit(&bout, 1, OWRITE);
+	for(b = glxy.a; b < glxy.a + new; b++)
+		Bprint(&bout, "%B\n", b);
+
+	for(b = glxy.a+new; b < glxy.a+glxy.l; b++) {
+		b->x += o.x;
+		b->y += o.y;
+		Bprint(&bout, "%B\n", b);
+	}
+	Bterm(&bout);
+
+	exits(nil);
+}
diff --git a/sys/src/games/galaxy/quad.c b/sys/src/games/galaxy/quad.c
new file mode 100644
index 000000000..e677e3628
--- /dev/null
+++ b/sys/src/games/galaxy/quad.c
@@ -0,0 +1,140 @@
+#include <u.h>
+#include <libc.h>
+#include <draw.h>
+#include "galaxy.h"
+
+void
+growquads(void)
+{
+	quads.max *= 2;
+	quads.a = realloc(quads.a, sizeof(Quad) * quads.max);
+	if(quads.a == nil)
+		sysfatal("could not realloc quads: %r");
+}
+
+Quad*
+quad(Body *b)
+{
+	Quad *q;
+
+	if(quads.l == quads.max)
+		return nil;
+
+	q = quads.a + quads.l++;
+	memset(q->c, 0, sizeof(QB)*4);
+	q->x = b->x;
+	q->y = b->y;
+	q->mass = b->mass;
+	return q;
+}
+
+int
+quadins(Body *nb, double size)
+{
+	QB *qb;
+	Quad *q;
+	Body *b;
+	double mass, qx, qy;
+	uint qxy;
+	int d;
+
+	if(space.t == EMPTY) {
+		space.t = BODY;
+		space.b = nb;
+		return 0;
+	}
+
+	d = 0;
+	qb = &space;
+	qx = 0.0;
+	qy = 0.0;
+	for(;;) {
+		if(qb->t == BODY) {
+			b = qb->b;
+			qxy = b->x < qx ? 0 : 1;
+			qxy |= b->y < qy ? 0 : 2;
+			qb->t = QUAD;
+			if((qb->q = quad(b)) == nil)
+				return -1;
+			qb->q->c[qxy].t = BODY;
+			qb->q->c[qxy].b = b;
+		}
+
+		q = qb->q;
+		mass = q->mass + nb->mass;
+		q->x = (q->x*q->mass + nb->x*nb->mass) / mass;
+		q->y = (q->y*q->mass + nb->y*nb->mass) / mass;
+		q->mass = mass;
+
+		qxy = nb->x < qx ? 0 : 1;
+		qxy |= nb->y < qy ? 0 : 2;
+		if(q->c[qxy].t == EMPTY) {
+			q->c[qxy].t = BODY;
+			q->c[qxy].b = nb;
+			STATS(if(d > quaddepth) quaddepth = d;)
+			return 0;
+		}
+
+		STATS(d++;)
+		qb = &q->c[qxy];
+		size /= 2;
+		qx += qxy&1 ? size/2 : -size/2;
+		qy += qxy&2 ? size/2 : -size/2;
+	}
+}
+
+void
+quadcalc(QB qb, Body *b, double size)
+{
+	double fx÷❨m₁m₂❩, fy÷❨m₁m₂❩, dx, dy, h, G÷h³;
+
+	for(;;) switch(qb.t) {
+	default:
+		sysfatal("quadcalc: No such type");
+	case EMPTY:
+		return;
+	case BODY:
+		if(qb.b == b)
+			return;
+		dx = qb.b->x - b->x;
+		dy = qb.b->y - b->y;
+		h = hypot(hypot(dx, dy), ε);
+		G÷h³ = G / (h*h*h);
+		fx÷❨m₁m₂❩ = dx * G÷h³;
+		fy÷❨m₁m₂❩ = dy * G÷h³;
+		b->newa.x += fx÷❨m₁m₂❩ * qb.b->mass;
+		b->newa.y += fy÷❨m₁m₂❩ * qb.b->mass;
+		STATS(calcs++;)
+		return;
+	case QUAD:
+		dx = qb.q->x - b->x;
+		dy = qb.q->y - b->y;
+		h = hypot(dx, dy);
+		if(h != 0.0 && size/h < θ) {
+			h = hypot(h, ε);
+			G÷h³ = G / (h*h*h);
+			fx÷❨m₁m₂❩ = dx * G÷h³;
+			fy÷❨m₁m₂❩ = dy * G÷h³;
+			b->newa.x += fx÷❨m₁m₂❩ * qb.q->mass;
+			b->newa.y += fy÷❨m₁m₂❩ * qb.q->mass;
+			STATS(calcs++;)
+			return;
+		}
+		size /= 2;
+		quadcalc(qb.q->c[0], b, size);
+		quadcalc(qb.q->c[1], b, size);
+		quadcalc(qb.q->c[2], b, size);
+		qb = qb.q->c[3];
+		break;	/* quadcalc(q->q[3], b, size); */
+	}
+}
+
+void
+quadsinit(void)
+{
+	quads.a = calloc(5, sizeof(Body));
+	if(quads.a == nil)
+		sysfatal("could not allocate quads: %r");
+	quads.l = 0;
+	quads.max = 5;
+}
diff --git a/sys/src/games/mkfile b/sys/src/games/mkfile
index 8ec52df05..e1c470201 100644
--- a/sys/src/games/mkfile
+++ b/sys/src/games/mkfile
@@ -25,6 +25,7 @@ DIRS=\
 	blabs\
 	c64\
 	doom\
+	galaxy\
 	gb\
 	gba\
 	mahjongg\