games/galaxy: add n-body simulator
spew
parent
412b7501e4
commit
5aabf85d7c
|
|
@ -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).
|
||||
|
|
@ -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)
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -25,6 +25,7 @@ DIRS=\
|
|||
blabs\
|
||||
c64\
|
||||
doom\
|
||||
galaxy\
|
||||
gb\
|
||||
gba\
|
||||
mahjongg\
|
||||
|
|
|
|||
Loading…
Reference in New Issue