--- /dev/null
+.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).
--- /dev/null
+.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)
--- /dev/null
+#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);
+}
--- /dev/null
+#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);
+}
--- /dev/null
+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
--- /dev/null
+</$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
--- /dev/null
+#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);
+}
--- /dev/null
+#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;
+}
blabs\
c64\
doom\
+ galaxy\
gb\
gba\
mahjongg\
projects / plan9front.git / commitdiff