Add scat. Temporary fix to rc r.e. note groups.

18 files changed, 5659 insertions, 0 deletions

diff --git a/src/cmd/scat/bitinput.c b/src/cmd/scat/bitinput.c
new file mode 100644
index 00000000..b4bd286b
--- /dev/null
+++ b/src/cmd/scat/bitinput.c
@@ -0,0 +1,76 @@
+#include	<u.h>
+#include	<libc.h>
+#include	<bio.h>
+#include	"sky.h"
+
+static int hufvals[] = {
+	 1,  1,  1,  1,  1,  1,  1,  1,
+	 2,  2,  2,  2,  2,  2,  2,  2,
+	 4,  4,  4,  4,  4,  4,  4,  4,
+	 8,  8,  8,  8,  8,  8,  8,  8,
+	 3,  3,  3,  3,  5,  5,  5,  5,
+	10, 10, 10, 10, 12, 12, 12, 12,
+	15, 15, 15, 15,  6,  6,  7,  7,
+	 9,  9, 11, 11, 13, 13,  0, 14,
+};
+
+static int huflens[] = {
+	3, 3, 3, 3, 3, 3, 3, 3,
+	3, 3, 3, 3, 3, 3, 3, 3,
+	3, 3, 3, 3, 3, 3, 3, 3,
+	3, 3, 3, 3, 3, 3, 3, 3,
+	4, 4, 4, 4, 4, 4, 4, 4,
+	4, 4, 4, 4, 4, 4, 4, 4,
+	4, 4, 4, 4, 5, 5, 5, 5,
+	5, 5, 5, 5, 5, 5, 6, 6,
+};
+
+static	int	buffer;
+static	int	bits_to_go;		/* Number of bits still in buffer */
+
+void
+start_inputing_bits(void)
+{
+	bits_to_go = 0;
+}
+
+int
+input_huffman(Biobuf *infile)
+{
+	int c;
+
+	if(bits_to_go < 6) {
+		c = Bgetc(infile);
+		if(c < 0) {
+			fprint(2, "input_huffman: unexpected EOF\n");
+			exits("format");
+		}
+		buffer = (buffer<<8) | c;
+		bits_to_go += 8;
+	}
+	c = (buffer >> (bits_to_go-6)) & 0x3f;
+	bits_to_go -= huflens[c];
+	return hufvals[c];
+}
+
+int
+input_nybble(Biobuf *infile)
+{
+	int c;
+
+	if(bits_to_go < 4) {
+		c = Bgetc(infile);
+		if(c < 0){
+			fprint(2, "input_nybble: unexpected EOF\n");
+			exits("format");
+		}
+		buffer = (buffer<<8) | c;
+		bits_to_go += 8;
+	}
+
+	/*
+	 * pick off the first 4 bits
+	 */
+	bits_to_go -= 4;
+	return (buffer>>bits_to_go) & 0x0f;
+}
diff --git a/src/cmd/scat/desc.c b/src/cmd/scat/desc.c
new file mode 100644
index 00000000..968229cc
--- /dev/null
+++ b/src/cmd/scat/desc.c
@@ -0,0 +1,327 @@
+char *desctab[][2]={
+	"!!!",	"(magnificent or otherwise interesting object)",
+	"!!",	"(superb)",
+	"!",	"(remarkable)",
+	"1st",	"first",
+	"2nd",	"second",
+	"3rd",	"third",
+	"4th",	"fourth",
+	"5th",	"fifth",
+	"6th",	"sixth",
+	"7th",	"seventh",
+	"8th",	"eighth",
+	"9th",	"ninth",
+	"B",	"bright",
+	"Borealis",	"Borealis",
+	"C",	"compressed",
+	"Car",	"Car",
+	"Cas",	"Cas",
+	"Cen",	"Cen",
+	"Cl",	"cluster",
+	"Cyg",	"Cyg",
+	"D",	"double",
+	"Dra",	"Dra",
+	"Dumbbell",	"Dumbbell",
+	"E",	"extended",
+	"F",	"faint",
+	"L",	"large",
+	"Lyra",	"Lyra",
+	"M",	"(in the) middle",
+	"Merope",	"Merope",
+	"Milky",	"Milky",
+	"Mon",	"Mon",
+	"N",	"(to a) nucleus",
+	"Neb",	"Nebula",
+	"Nor",	"Nor",
+	"Nucl",	"nucleus",
+	"Nuclei",	"nuclei",
+	"P",	"poor in stars",
+	"PN",	"planetary nebula",
+	"Polarissima",	"Polarissima",
+	"Praesepe",	"Praesepe",
+	"Psc",	"Psc",
+	"R",	"round",
+	"RA",	"RA",
+	"RR",	"exactly round",
+	"Ri",	"rich in stars",
+	"S",	"small",
+	"Sco",	"Sco",
+	"S*",	"small star",
+	"Way",	"Way",
+	"ab",	"about",
+	"about",	"about",
+	"alm",	"almost",
+	"alpha",	"α",
+	"am",	"among",
+	"annul",	"annular",
+	"att",	"attached",
+	"b",	"brighter",
+	"beautiful",	"beautiful",
+	"bet",	"between",
+	"beta",	"β",
+	"bf",	"brightest to f side",
+	"biN",	"binuclear",
+	"bifid",	"bifid",
+	"bifurcated",	"bifurcated",
+	"black",	"black",
+	"blue",	"blue",
+	"bn",	"brightest to n side",
+	"border",	"border",
+	"bp",	"brightest to p side",
+	"branch",	"branch",
+	"branched",	"branched",
+	"branches",	"branches",
+	"bright",	"bright",
+	"brighter",	"brighter",
+	"brightest",	"brightest",
+	"brightness",	"brightness",
+	"brush",	"brush",
+	"bs",	"brightest to s side",
+	"but",	"but",
+	"by",	"by",
+	"c",	"considerably",
+	"centre",	"centre",
+	"certain",	"certain",
+	"chev",	"chevelure",
+	"chief",	"chief",
+	"chiefly",	"chiefly",
+	"circle",	"circle",
+	"close",	"close",
+	"cloud",	"cloud",
+	"cluster",	"cluster",
+	"clusters",	"clusters",
+	"co",	"coarse, coarsely",
+	"com",	"cometic",
+	"comet",	"comet",
+	"cometary",	"cometary",
+	"comp",	"companion",
+	"condens",	"condensations",
+	"condensed",	"condensed",
+	"conn",	"connected",
+	"connected",	"connected",
+	"connecting",	"connecting",
+	"cont",	"in contact",
+	"corner",	"corner",
+	"curved",	"curved",
+	"d",	"diameter",
+	"decl",	"declination",
+	"def",	"defined",
+	"defect",	"defect",
+	"deg",	"°",
+	"delta",	"δ",
+	"dense",	"dense",
+	"densest",	"densest",
+	"descr",	"description",
+	"description",	"description",
+	"dif",	"diffuse",
+	"different",	"different",
+	"diffic",	"difficult",
+	"difficult",	"difficult",
+	"diffuse",	"diffuse",
+	"diffused",	"diffused",
+	"disc",	"disc",
+	"dist",	"distant",
+	"distant",	"distant",
+	"distinct",	"distinct",
+	"double",	"double",
+	"doubtful",	"doubtful",
+	"dozen",	"dozen",
+	"e",	"extremely",
+	"each",	"each",
+	"edge",	"edge",
+	"ee",	"most extremely",
+	"ellipt",	"elliptical",
+	"elliptic",	"elliptical",
+	"end",	"end",
+	"ends",	"ends",
+	"epsilon",	"ε",
+	"equilateral",	"equilateral",
+	"er",	"easily resolvable",
+	"eta",	"η",
+	"evidently",	"evidently",
+	"exc",	"excentric",
+	"excen",	"excentric",
+	"excent",	"excentric",
+	"excentric",	"excentric",
+	"extends",	"extends",
+	"f",	"following",
+	"faint",	"faint",
+	"fainter",	"fainter",
+	"faintest",	"faintest",
+	"falcate",	"falcate",
+	"fan",	"fan",
+	"farther",	"farther",
+	"field",	"field",
+	"fine",	"fine",
+	"forming",	"forming",
+	"forms",	"forms",
+	"found",	"found",
+	"from",	"from",
+	"g",	"gradually",
+	"gamma",	"γ",
+	"gaseous",	"gaseous",
+	"gl",	"gradually a little",
+	"glob. cl.",	"globular cluster",
+	"gr",	"group",
+	"great",	"great",
+	"greater",	"greater",
+	"group",	"group",
+	"groups",	"groups",
+	"i",	"irregular",
+	"iF",	"irregular figure",
+	"if",	"if",
+	"in",	"in",
+	"indistinct",	"indistinct",
+	"incl",	"including",
+	"inv",	"involved",
+	"iota",	"ι",
+	"irr",	"irregular",
+	"is",	"is",
+	"it",	"it",
+	"kappa",	"κ",
+	"l",	"little, long",
+	"lC",	"little compressed",
+	"lE",	"little extended",
+	"lambda",	"λ",
+	"larger",	"larger",
+	"last",	"last",
+	"lb",	"little brighter",
+	"least",	"least",
+	"like",	"like",
+	"line",	"in a line",
+	"little",	"little",
+	"long",	"long",
+	"looks",	"looks",
+	"looped",	"looped",
+	"m",	"magnitude",
+	"mE",	"much extended",
+	"mag",	"mag",
+	"makes",	"makes",
+	"many",	"many",
+	"mb",	"much brighter",
+	"more",	"more",
+	"mottled",	"mottled",
+	"mu",	"μ",
+	"mult",	"multiple",
+	"n",	"north",
+	"narrow",	"narrow",
+	"near",	"near",
+	"nearly",	"nearly",
+	"neb",	"nebula",
+	"nebs",	"nebulous",
+	"nebula",	"nebula",
+	"nebulosity",	"nebulosity",
+	"nebulous",	"nebulous",
+	"neby",	"nebulosity",
+	"nf",	"north following",
+	"no",	"no",
+	"nonexistent",	"nonexistent",
+	"not",	"not",
+	"np",	"north preceding",
+	"nr",	"near",
+	"ns",	"north-south",
+	"nu",	"ν",
+	"omega",	"ω",
+	"p",	"preceding",
+	"pB",	"pretty bright",
+	"pC",	"pretty compressed",
+	"pF",	"pretty faint",
+	"pL",	"pretty large",
+	"pR",	"pretty round",
+	"pS",	"pretty small",
+	"parallel",	"parallel",
+	"part",	"part",
+	"partly",	"partly",
+	"patch",	"patch",
+	"patches",	"patches",
+	"perhaps",	"perhaps",
+	"perpendicular",	"perpendicular",
+	"pf",	"preceding-following",
+	"pg",	"pretty gradually",
+	"photosphere",	"photosphere",
+	"pi",	"π",
+	"place",	"place",
+	"plate",	"plate",
+	"plan",	"planetary nebula",
+	"pointed",	"pointed",
+	"portion",	"portion",
+	"pos",	"position angle",
+	"possibly",	"possibly",
+	"prob",	"probably",
+	"probably",	"probably",
+	"ps",	"pretty suddenly",
+	"r",	"mottled",
+	"requires",	"requires",
+	"resolved",	"resolved",
+	"rho",	"ρ",
+	"ring",	"ring",
+	"rough",	"rough",
+	"rr",	"some stars seen",
+	"rrr",	"clearly consisting of stars",
+	"ruby",	"ruby",
+	"s",	"south",
+	"same",	"same",
+	"sb",	"suddenly brighter",
+	"sc",	"scattered",
+	"second",	"second",
+	"seems",	"seems",
+	"seen",	"seen",
+	"segment",	"segment",
+	"semi",	"semi",
+	"sev",	"several",
+	"several",	"several",
+	"sf",	"south following",
+	"shape",	"shape",
+	"shaped",	"shaped",
+	"sharp",	"sharp",
+	"sigma",	"σ",
+	"sl",	"suddenly a little",
+	"slightly",	"slightly",
+	"small",	"small",
+	"south",	"south",
+	"sp",	"south preceding",
+	"spectrum",	"spectrum",
+	"spindle",	"spindle",
+	"spir",	"spiral",
+	"spiral",	"spiral",
+	"st 9...",	"stars of mag. 9 and fainter",
+	"st 9...13",	"stars of mag. 9 to 13",
+	"st",	"stars",
+	"stell",	"stellar, pointlike",
+	"stellar",	"stellar",
+	"straight",	"straight",
+	"streak",	"streak",
+	"strongly",	"strongly",
+	"surrounded",	"surrounded",
+	"surrounds",	"surrounds",
+	"susp",	"suspected",
+	"suspected",	"suspected",
+	"tau",	"τ",
+	"theta",	"θ",
+	"trap",	"trapezium",
+	"trapez",	"trapezium",
+	"trapezium",	"trapezium",
+	"triN",	"trinuclear",
+	"v",	"very",
+	"var",	"variable",
+	"variable",	"variable",
+	"verification",	"verification",
+	"verified",	"verified",
+	"very",	"very",
+	"vl",	"very little",
+	"vm",	"very much",
+	"vs",	"very suddenly",
+	"vv",	"very very",
+	"zeta",	"ζ",
+	0,	0,
+};
+
+/*&
+	"ST9",	"stars from the 9th magnitude downward",
+	"ST9...13",	"stars from 9th to 13th magnitude",
+	"()",	"items questioned by Dreyer enclosed in parentheses",
+	*10	star of 10th mag
+	"*",	"a star (or stars)",
+	"**",	"double star",
+	"***",	"triple star",
+*/
diff --git a/src/cmd/scat/display.c b/src/cmd/scat/display.c
new file mode 100644
index 00000000..11642e59
--- /dev/null
+++ b/src/cmd/scat/display.c
@@ -0,0 +1,88 @@
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#include <draw.h>
+#include "sky.h"
+
+void
+displaypic(Picture *pic)
+{
+	int p[2];
+	int i, n;
+	uchar *a;
+	
+
+	if(pipe(p) < 0){
+		fprint(2, "pipe failed: %r\n");
+		return;
+	}
+	switch(rfork(RFPROC|RFFDG)){
+	case -1:
+		fprint(2, "fork failed: %r\n");
+		return;
+
+	case 0:
+		close(p[1]);
+		dup(p[0], 0);
+		close(p[0]);
+	//	execl("/bin/page", "page", "-w", 0);
+		execlp("img", "img", 0);
+		fprint(2, "exec failed: %r\n");
+		exits("exec");
+
+	default:
+		close(p[0]);
+		fprint(p[1], "%11s %11d %11d %11d %11d ",
+			"k8", pic->minx, pic->miny, pic->maxx, pic->maxy);
+		n = (pic->maxx-pic->minx)*(pic->maxy-pic->miny);
+		/* release the memory as we hand it off; this could be a big piece of data */
+		a = pic->data;
+		while(n > 0){
+			i = 8192 - (((int)a)&8191);
+			if(i > n)
+				i = n;
+			if(write(p[1], a, i)!=i)
+				fprint(2, "write error: %r\n");
+		//	if(i == 8192)	/* page aligned */
+		//		segfree(a, i);
+			n -= i;
+			a += i;
+		}
+		free(pic->data);
+		free(pic);
+		close(p[1]);
+		break;
+	}
+}
+
+void
+displayimage(Image *im)
+{
+	int p[2];	
+
+	if(pipe(p) < 0){
+		fprint(2, "pipe failed: %r\n");
+		return;
+	}
+	switch(rfork(RFPROC|RFFDG)){
+	case -1:
+		fprint(2, "fork failed: %r\n");
+		return;
+
+	case 0:
+		close(p[1]);
+		dup(p[0], 0);
+		close(p[0]);
+		execlp("img", "img", 0);
+	//	execl("/bin/page", "page", "-w", 0);
+		fprint(2, "exec failed: %r\n");
+		exits("exec");
+
+	default:
+		close(p[0]);
+		writeimage(p[1], im, 0);
+		freeimage(im);
+		close(p[1]);
+		break;
+	}
+}
diff --git a/src/cmd/scat/dssread.c b/src/cmd/scat/dssread.c
new file mode 100644
index 00000000..06f0135e
--- /dev/null
+++ b/src/cmd/scat/dssread.c
@@ -0,0 +1,113 @@
+#include	<u.h>
+#include	<libc.h>
+#include	<bio.h>
+#include	"sky.h"
+
+static	void	dodecode(Biobuf*, Pix*, int, int, uchar*);
+static	long	getlong(uchar*);
+int	debug;
+
+Img*
+dssread(char *file)
+{
+	int nx, ny, scale, sumall;
+	Pix  *p, *pend;
+	uchar buf[21];
+	Biobuf *bp;
+	Img *ip;
+
+	if(debug)
+		Bprint(&bout, "reading %s\n", file);
+	bp = Bopen(file, OREAD);
+	if(bp == 0)
+		return 0;
+	if(Bread(bp, buf, sizeof(buf)) != sizeof(buf) ||
+	   buf[0] != 0xdd || buf[1] != 0x99){
+		werrstr("bad format");
+		return 0;
+	}
+	nx = getlong(buf+2);
+	ny = getlong(buf+6);
+	scale = getlong(buf+10);
+	sumall = getlong(buf+14);
+	if(debug)
+		fprint(2, "%s: nx=%d, ny=%d, scale=%d, sumall=%d, nbitplanes=%d,%d,%d\n",
+			file, nx, ny, scale, sumall, buf[18], buf[19], buf[20]);
+	ip = malloc(sizeof(Img) + (nx*ny-1)*sizeof(int));
+	if(ip == 0){
+		Bterm(bp);
+		werrstr("no memory");
+		return 0;
+	}
+	ip->nx = nx;
+	ip->ny = ny;
+	dodecode(bp, ip->a, nx, ny, buf+18);
+	ip->a[0] = sumall;	/* sum of all pixels */
+	Bterm(bp);
+	if(scale > 1){
+		p = ip->a;
+		pend = &ip->a[nx*ny];
+		while(p < pend)
+			*p++ *= scale;
+	}
+	hinv(ip->a, nx, ny);
+	return ip;
+}
+
+static
+void
+dodecode(Biobuf *infile, Pix  *a, int nx, int ny, uchar *nbitplanes)
+{
+	int nel, nx2, ny2, bits, mask;
+	Pix *aend, px;
+
+	nel = nx*ny;
+	nx2 = (nx+1)/2;
+	ny2 = (ny+1)/2;
+	memset(a, 0, nel*sizeof(*a));
+
+	/*
+	 * Initialize bit input
+	 */
+	start_inputing_bits();
+
+	/*
+	 * read bit planes for each quadrant
+	 */
+	qtree_decode(infile, &a[0],          ny, nx2,  ny2,  nbitplanes[0]);
+	qtree_decode(infile, &a[ny2],        ny, nx2,  ny/2, nbitplanes[1]);
+	qtree_decode(infile, &a[ny*nx2],     ny, nx/2, ny2,  nbitplanes[1]);
+	qtree_decode(infile, &a[ny*nx2+ny2], ny, nx/2, ny/2, nbitplanes[2]);
+
+	/*
+	 * make sure there is an EOF symbol (nybble=0) at end
+	 */
+	if(input_nybble(infile) != 0){
+		fprint(2, "dodecode: bad bit plane values\n");
+		exits("format");
+	}
+
+	/*
+	 * get the sign bits
+	 */
+	aend = &a[nel];
+	mask = 0;
+	bits = 0;;
+	for(; a<aend; a++) {
+		if(px = *a) {
+			if(mask == 0) {
+				mask = 0x80;
+				bits = Bgetc(infile);
+			}
+			if(mask & bits)
+				*a = -px;
+			mask >>= 1;
+		}
+	}
+}
+
+static
+long	getlong(uchar *p)
+{
+	return (p[0]<<24) | (p[1]<<16) | (p[2]<<8) | p[3];
+}
diff --git a/src/cmd/scat/header.c b/src/cmd/scat/header.c
new file mode 100644
index 00000000..519d98f9
--- /dev/null
+++ b/src/cmd/scat/header.c
@@ -0,0 +1,247 @@
+#include	<u.h>
+#include	<libc.h>
+#include	<bio.h>
+#include	"sky.h"
+
+struct
+{
+	char	name[9];
+	char	offset;
+} Hproto[] =
+{
+	"ppo1",		Pppo1,
+	"ppo2",		Pppo2,
+	"ppo3",		Pppo3,
+	"ppo4",		Pppo4,
+	"ppo5",		Pppo5,
+	"ppo6",		Pppo6,
+
+	"amdx1",	Pamdx1,
+	"amdx2",	Pamdx2,
+	"amdx3",	Pamdx3,
+	"amdx4",	Pamdx4,
+	"amdx5",	Pamdx5,
+	"amdx6",	Pamdx6,
+	"amdx7",	Pamdx7,
+	"amdx8",	Pamdx8,
+	"amdx9",	Pamdx9,
+	"amdx10",	Pamdx10,
+	"amdx11",	Pamdx11,
+	"amdx12",	Pamdx12,
+	"amdx13",	Pamdx13,
+	"amdx14",	Pamdx14,
+	"amdx15",	Pamdx15,
+	"amdx16",	Pamdx16,
+	"amdx17",	Pamdx17,
+	"amdx18",	Pamdx18,
+	"amdx19",	Pamdx19,
+	"amdx20",	Pamdx20,
+
+	"amdy1",	Pamdy1,
+	"amdy2",	Pamdy2,
+	"amdy3",	Pamdy3,
+	"amdy4",	Pamdy4,
+	"amdy5",	Pamdy5,
+	"amdy6",	Pamdy6,
+	"amdy7",	Pamdy7,
+	"amdy8",	Pamdy8,
+	"amdy9",	Pamdy9,
+	"amdy10",	Pamdy10,
+	"amdy11",	Pamdy11,
+	"amdy12",	Pamdy12,
+	"amdy13",	Pamdy13,
+	"amdy14",	Pamdy14,
+	"amdy15",	Pamdy15,
+	"amdy16",	Pamdy16,
+	"amdy17",	Pamdy17,
+	"amdy18",	Pamdy18,
+	"amdy19",	Pamdy19,
+	"amdy20",	Pamdy20,
+
+	"pltscale",	Ppltscale,
+	"xpixelsz",	Pxpixelsz,
+	"ypixelsz",	Pypixelsz,
+
+	"pltrah",	Ppltrah,
+	"pltram",	Ppltram,
+	"pltras",	Ppltras,
+	"pltdecd",	Ppltdecd,
+	"pltdecm",	Ppltdecm,
+	"pltdecs",	Ppltdecs,
+
+};
+
+Header*
+getheader(char *rgn)
+{
+	char rec[81], name[81], value[81];
+	char *p;
+	Biobuf *bin;
+	Header hd, *h;
+	int i, j, decsn, dss;
+
+	dss = 0;
+	sprint(rec, "/lib/sky/dssheaders/%s.hhh", rgn);
+	bin = Bopen(unsharp(rec), OREAD);
+/*
+	if(bin == 0) {
+		dss = 102;
+		sprint(rec, "/n/juke/dss/dss.102/headers/%s.hhh", rgn);
+		bin = Bopen(rec, OREAD);
+	}
+	if(bin == 0) {
+		dss = 61;
+		sprint(rec, "/n/juke/dss/dss.061/headers/%s.hhh", rgn);
+		bin = Bopen(rec, OREAD);
+	}
+*/
+	if(bin == 0) {
+		fprint(2, "cannot open %s\n", rgn);
+		exits("file");
+	}
+	if(debug)
+		Bprint(&bout, "reading %s\n", rec);
+	if(dss)
+		Bprint(&bout, "warning: reading %s from jukebox\n", rec);
+
+	memset(&hd, 0, sizeof(hd));
+	j = 0;
+	decsn = 0;
+	for(;;) {
+		if(dss) {
+			if(Bread(bin, rec, 80) != 80)
+				break;
+			rec[80] = 0;
+		} else {
+			p = Brdline(bin, '\n');
+			if(p == 0)
+				break;
+			p[Blinelen(bin)-1] = 0;
+			strcpy(rec, p);
+		}
+
+		p = strchr(rec, '/');
+		if(p)
+			*p = 0;
+		p = strchr(rec, '=');
+		if(p == 0)
+			continue;
+		*p++ = 0;
+		if(getword(name, rec) == 0)
+			continue;
+		if(getword(value, p) == 0)
+			continue;
+		if(strcmp(name, "pltdecsn") == 0) {
+			if(strchr(value, '-'))
+				decsn = 1;
+			continue;
+		}
+		for(i=0; i<nelem(Hproto); i++) {
+			j++;
+			if(j >= nelem(Hproto))
+				j = 0;
+			if(strcmp(name, Hproto[j].name) == 0) {
+				hd.param[(uchar)Hproto[j].offset] = atof(value);
+				break;
+			}
+		}
+	}
+	Bterm(bin);
+
+	hd.param[Ppltra] = RAD(hd.param[Ppltrah]*15 +
+		hd.param[Ppltram]/4 + hd.param[Ppltras]/240);
+	hd.param[Ppltdec] = RAD(hd.param[Ppltdecd] +
+		hd.param[Ppltdecm]/60 + hd.param[Ppltdecs]/3600);
+	if(decsn)
+		hd.param[Ppltdec] = -hd.param[Ppltdec];
+	hd.amdflag = 0;
+	for(i=Pamdx1; i<=Pamdx20; i++)
+		if(hd.param[i] != 0) {
+			hd.amdflag = 1;
+			break;
+		}
+	h = malloc(sizeof(*h));
+	*h = hd;
+	return h;
+}
+
+void
+getplates(void)
+{
+	char rec[81], *q;
+	Plate *p;
+	Biobuf *bin;
+	int c, i, dss;
+
+	dss = 0;
+	sprint(rec, "/lib/sky/dssheaders/lo_comp.lis");
+	bin = Bopen(rec, OREAD);
+	if(bin == 0) {
+		dss = 102;
+		sprint(rec, "%s/headers/lo_comp.lis", dssmount(dss));
+		bin = Bopen(rec, OREAD);
+	}
+	if(bin == 0) {
+		dss = 61;
+		sprint(rec, "%s/headers/lo_comp.lis", dssmount(dss));
+		bin = Bopen(rec, OREAD);
+	}
+	if(bin == 0) {
+		fprint(2, "can't open lo_comp.lis; try 9fs juke\n");
+		exits("open");
+	}
+	if(debug)
+		Bprint(&bout, "reading %s\n", rec);
+	if(dss)
+		Bprint(&bout, "warning: reading %s from jukebox\n", rec);
+	for(nplate=0;;) {
+		if(dss) {
+			if(Bread(bin, rec, 80) != 80)
+				break;
+			rec[80] = 0;
+		} else {
+			q = Brdline(bin, '\n');
+			if(q == 0)
+				break;
+			q[Blinelen(bin)-1] = 0;
+			strcpy(rec, q);
+		}
+		if(rec[0] == '#')
+			continue;
+		if(nplate < nelem(plate)) {
+			p = &plate[nplate];
+			memmove(p->rgn, rec+0, 5);
+			if(p->rgn[4] == ' ')
+				p->rgn[4] = 0;
+			for(i=0; c=p->rgn[i]; i++)
+				if(c >= 'A' && c <= 'Z')
+					p->rgn[i] += 'a'-'A';
+			p->ra = RAD(atof(rec+12)*15 +
+				atof(rec+15)/4 +
+				atof(rec+18)/240);
+			p->dec = RAD(atof(rec+26) +
+				atof(rec+29)/60 +
+				atof(rec+32)/3600);
+			if(rec[25] == '-')
+				p->dec = -p->dec;
+			p->disk = atoi(rec+53);
+			if(p->disk == 0)
+				continue;
+		}
+		nplate++;
+	}
+	Bterm(bin);
+
+	if(nplate >= nelem(plate))
+		fprint(2, "nplate too small %d %d\n", nelem(plate), nplate);
+	if(debug)
+		Bprint(&bout, "%d plates\n", nplate);
+}
+
+char*
+dssmount(int dskno)
+{
+	Bprint(&bout, "not mounting dss\n");
+	return "/n/dss";
+}
+
diff --git a/src/cmd/scat/hinv.c b/src/cmd/scat/hinv.c
new file mode 100644
index 00000000..c76b00ea
--- /dev/null
+++ b/src/cmd/scat/hinv.c
@@ -0,0 +1,231 @@
+#include	<u.h>
+#include	<libc.h>
+#include	<bio.h>
+#include	"sky.h"
+
+static	void	unshuffle(Pix*, int, int, Pix*);
+static	void	unshuffle1(Pix*, int, Pix*);
+
+void
+hinv(Pix *a, int nx, int ny)
+{
+	int nmax, log2n, h0, hx, hy, hc, i, j, k;
+	int nxtop, nytop, nxf, nyf, c;
+	int oddx, oddy;
+	int shift;
+	int s10, s00;
+	Pix *tmp;
+
+	/*
+	 * log2n is log2 of max(nx, ny) rounded up to next power of 2
+	 */
+	nmax = ny;
+	if(nx > nmax)
+		nmax = nx;
+	log2n = log(nmax)/LN2 + 0.5;
+	if(nmax > (1<<log2n))
+		log2n++;
+
+	/*
+	 * get temporary storage for shuffling elements
+	 */
+	tmp = (Pix*)malloc(((nmax+1)/2) * sizeof(*tmp));
+	if(tmp == nil) {
+		fprint(2, "hinv: insufficient memory\n");
+		exits("memory");
+	}
+
+	/*
+	 * do log2n expansions
+	 *
+	 * We're indexing a as a 2-D array with dimensions (nx,ny).
+	 */
+	shift = 1;
+	nxtop = 1;
+	nytop = 1;
+	nxf = nx;
+	nyf = ny;
+	c = 1<<log2n;
+	for(k = log2n-1; k>=0; k--) {
+		/*
+		 * this somewhat cryptic code generates the sequence
+		 * ntop[k-1] = (ntop[k]+1)/2, where ntop[log2n] = n
+		 */
+		c = c>>1;
+		nxtop = nxtop<<1;
+		nytop = nytop<<1;
+		if(nxf <= c)
+			nxtop--;
+		else
+			nxf -= c;
+		if(nyf <= c)
+			nytop--;
+		else
+			nyf -= c;
+
+		/*
+		 * halve divisors on last pass
+		 */
+		if(k == 0)
+			shift = 0;
+
+		/*
+		 * unshuffle in each dimension to interleave coefficients
+		 */
+		for(i = 0; i<nxtop; i++)
+			unshuffle1(&a[ny*i], nytop, tmp);
+		for(j = 0; j<nytop; j++)
+			unshuffle(&a[j], nxtop, ny, tmp);
+		oddx = nxtop & 1;
+		oddy = nytop & 1;
+		for(i = 0; i<nxtop-oddx; i += 2) {
+			s00 = ny*i;			/* s00 is index of a[i,j]	*/
+			s10 = s00+ny;			/* s10 is index of a[i+1,j]	*/
+			for(j = 0; j<nytop-oddy; j += 2) {
+				/*
+				 * Multiply h0,hx,hy,hc by 2 (1 the last time through).
+				 */
+				h0 = a[s00  ] << shift;
+				hx = a[s10  ] << shift;
+				hy = a[s00+1] << shift;
+				hc = a[s10+1] << shift;
+
+				/*
+				 * Divide sums by 4 (shift right 2 bits).
+				 * Add 1 to round -- note that these values are always
+				 * positive so we don't need to do anything special
+				 * for rounding negative numbers.
+				 */
+				a[s10+1] = (h0 + hx + hy + hc + 2) >> 2;
+				a[s10  ] = (h0 + hx - hy - hc + 2) >> 2;
+				a[s00+1] = (h0 - hx + hy - hc + 2) >> 2;
+				a[s00  ] = (h0 - hx - hy + hc + 2) >> 2;
+				s00 += 2;
+				s10 += 2;
+			}
+			if(oddy) {
+				/*
+				 * do last element in row if row length is odd
+				 * s00+1, s10+1 are off edge
+				 */
+				h0 = a[s00  ] << shift;
+				hx = a[s10  ] << shift;
+				a[s10  ] = (h0 + hx + 2) >> 2;
+				a[s00  ] = (h0 - hx + 2) >> 2;
+			}
+		}
+		if(oddx) {
+			/*
+			 * do last row if column length is odd
+			 * s10, s10+1 are off edge
+			 */
+			s00 = ny*i;
+			for(j = 0; j<nytop-oddy; j += 2) {
+				h0 = a[s00  ] << shift;
+				hy = a[s00+1] << shift;
+				a[s00+1] = (h0 + hy + 2) >> 2;
+				a[s00  ] = (h0 - hy + 2) >> 2;
+				s00 += 2;
+			}
+			if(oddy) {
+				/*
+				 * do corner element if both row and column lengths are odd
+				 * s00+1, s10, s10+1 are off edge
+				 */
+				h0 = a[s00  ] << shift;
+				a[s00  ] = (h0 + 2) >> 2;
+			}
+		}
+	}
+	free(tmp);
+}
+
+static
+void
+unshuffle(Pix *a, int n, int n2, Pix *tmp)
+{
+	int i;
+	int nhalf, twon2, n2xnhalf;
+	Pix *p1, *p2, *pt;
+
+	twon2 = n2<<1;
+	nhalf = (n+1)>>1;
+	n2xnhalf = n2*nhalf;		/* pointer to a[i] */
+
+	/*
+	 * copy 2nd half of array to tmp
+	 */
+	pt = tmp;
+	p1 = &a[n2xnhalf];		/* pointer to a[i] */
+	for(i=nhalf; i<n; i++) {
+		*pt = *p1;
+		pt++;
+		p1 += n2;
+	}
+
+	/*
+	 * distribute 1st half of array to even elements
+	 */
+	p2 = &a[n2xnhalf];		/* pointer to a[i] */
+	p1 = &a[n2xnhalf<<1];		/* pointer to a[2*i] */
+	for(i=nhalf-1; i>=0; i--) {
+		p1 -= twon2;
+		p2 -= n2;
+		*p1 = *p2;
+	}
+
+	/*
+	 * now distribute 2nd half of array (in tmp) to odd elements
+	 */
+	pt = tmp;
+	p1 = &a[n2];			/* pointer to a[i] */
+	for(i=1; i<n; i+=2) {
+		*p1 = *pt;
+		p1 += twon2;
+		pt++;
+	}
+}
+
+static
+void
+unshuffle1(Pix *a, int n, Pix *tmp)
+{
+	int i;
+	int nhalf;
+	Pix *p1, *p2, *pt;
+
+	nhalf = (n+1) >> 1;
+
+	/*
+	 * copy 2nd half of array to tmp
+	 */
+	pt = tmp;
+	p1 = &a[nhalf];				/* pointer to a[i]			*/
+	for(i=nhalf; i<n; i++) {
+		*pt = *p1;
+		pt++;
+		p1++;
+	}
+
+	/*
+	 * distribute 1st half of array to even elements
+	 */
+	p2 = &a[nhalf];		/* pointer to a[i]			*/
+	p1 = &a[nhalf<<1];		/* pointer to a[2*i]		*/
+	for(i=nhalf-1; i>=0; i--) {
+		p1 -= 2;
+		p2--;
+		*p1 = *p2;
+	}
+
+	/*
+	 * now distribute 2nd half of array (in tmp) to odd elements
+	 */
+	pt = tmp;
+	p1 = &a[1];					/* pointer to a[i]			*/
+	for(i=1; i<n; i+=2) {
+		*p1 = *pt;
+		p1 += 2;
+		pt++;
+	}
+}
diff --git a/src/cmd/scat/image.c b/src/cmd/scat/image.c
new file mode 100644
index 00000000..e61f0106
--- /dev/null
+++ b/src/cmd/scat/image.c
@@ -0,0 +1,158 @@
+#include	<u.h>
+#include	<libc.h>
+#include	<bio.h>
+#include	"sky.h"
+
+char	rad28[] = "0123456789abcdefghijklmnopqr";
+
+Picture*
+image(Angle ra, Angle dec, Angle wid, Angle hig)
+{
+	Pix *p;
+	uchar *b, *up;
+	int i, j, sx, sy, x, y;
+	char file[50];
+	Picture *pic;
+	Img* ip;
+	int lowx, lowy, higx, higy;
+	int slowx, slowy, shigx, shigy;
+	Header *h;
+	Angle d, bd;
+	Plate *pp, *bp;
+
+	if(gam.gamma == 0)
+		gam.gamma = -1;
+	if(gam.max == gam.min) {
+		gam.max = 17600;
+		gam.min = 2500;
+	}
+	gam.absgamma = gam.gamma;
+	gam.neg = 0;
+	if(gam.absgamma < 0) {
+		gam.absgamma = -gam.absgamma;
+		gam.neg = 1;
+	}
+	gam.mult1 = 1. / (gam.max - gam.min);
+	gam.mult2 = 255. * gam.mult1;
+
+	if(nplate == 0)
+		getplates();
+
+	bp = 0;
+	bd = 0;
+	for(i=0; i<nplate; i++) {
+		pp = &plate[i];
+		d = dist(ra, dec, pp->ra, pp->dec);
+		if(bp == 0 || d < bd) {
+			bp = pp;
+			bd = d;
+		}
+	}
+
+	if(debug)
+		Bprint(&bout, "best plate: %s %s disk %d %s\n",
+			hms(bp->ra), dms(bp->dec),
+			bp->disk, bp->rgn);
+
+	h = getheader(bp->rgn);
+	xypos(h, ra, dec, 0, 0);
+	if(wid <= 0 || hig <= 0) {
+		lowx = h->x;
+		lowy = h->y;
+		lowx = (lowx/500) * 500;
+		lowy = (lowy/500) * 500;
+		higx = lowx + 500;
+		higy = lowy + 500;
+	} else {
+		lowx = h->x - wid*ARCSECONDS_PER_RADIAN*1000 /
+			(h->param[Pxpixelsz]*h->param[Ppltscale]*2);
+		lowy = h->y - hig*ARCSECONDS_PER_RADIAN*1000 /
+			(h->param[Pypixelsz]*h->param[Ppltscale]*2);
+		higx = h->x + wid*ARCSECONDS_PER_RADIAN*1000 /
+			(h->param[Pxpixelsz]*h->param[Ppltscale]*2);
+		higy = h->y + hig*ARCSECONDS_PER_RADIAN*1000 /
+			(h->param[Pypixelsz]*h->param[Ppltscale]*2);
+	}
+	free(h);
+
+	if(lowx < 0) lowx = 0;
+	if(higx < 0) higx = 0;
+	if(lowy < 0) lowy = 0;
+	if(higy < 0) higy = 0;
+	if(lowx > 14000) lowx = 14000;
+	if(higx > 14000) higx = 14000;
+	if(lowy > 14000) lowy = 14000;
+	if(higy > 14000) higy = 14000;
+
+	if(debug)
+		Bprint(&bout, "xy on plate: %d,%d %d,%d\n",
+			lowx,lowy, higx, higy);
+
+	if(lowx >= higx || lowy >=higy) {
+		Bprint(&bout, "no image found\n");
+		return 0;
+	}
+
+	b = malloc((higx-lowx)*(higy-lowy)*sizeof(*b));
+	if(b == 0) {
+ emalloc:
+		fprint(2, "malloc error\n");
+		return 0;
+	}
+	memset(b, 0, (higx-lowx)*(higy-lowy)*sizeof(*b));
+
+	slowx = lowx/500;
+	shigx = (higx-1)/500;
+	slowy = lowy/500;
+	shigy = (higy-1)/500;
+
+	for(sx=slowx; sx<=shigx; sx++)
+	for(sy=slowy; sy<=shigy; sy++) {
+		if(sx < 0 || sx >= nelem(rad28) || sy < 0 || sy >= nelem(rad28)) {
+			fprint(2, "bad subplate %d %d\n", sy, sx);
+			free(b);
+			return 0;
+		}
+		sprint(file, "%s/%s/%s.%c%c",
+			dssmount(bp->disk),
+			bp->rgn, bp->rgn,
+			rad28[sy],
+			rad28[sx]);
+
+		ip = dssread(file);
+		if(ip == 0) {
+			fprint(2, "can't read %s: %r\n", file);
+			free(b);
+			return 0;
+		}
+
+		x = sx*500;
+		y = sy*500;
+		for(j=0; j<ip->ny; j++) {
+			if(y+j < lowy || y+j >= higy)
+				continue;
+			p = &ip->a[j*ip->ny];
+			up = b + (higy - (y+j+1))*(higx-lowx) + (x - lowx);
+			for(i=0; i<ip->nx; i++) {
+				if(x+i >= lowx && x+i < higx)
+					*up = dogamma(*p);
+				up++;
+				p += 1;
+			}
+		}
+		free(ip);
+	}
+
+	pic = malloc(sizeof(Picture));
+	if(pic == 0){
+		free(b);
+		goto emalloc;
+	}
+	pic->minx = lowx;
+	pic->miny = lowy;
+	pic->maxx = higx;
+	pic->maxy = higy;
+	pic->data = b;
+	strcpy(pic->name, bp->rgn);
+	return pic;
+}
diff --git a/src/cmd/scat/mkfile b/src/cmd/scat/mkfile
new file mode 100644
index 00000000..c0258a83
--- /dev/null
+++ b/src/cmd/scat/mkfile
@@ -0,0 +1,31 @@
+<$PLAN9/src/mkhdr
+
+TARG=scat
+OFILES=scat.$O\
+	bitinput.$O\
+	desc.$O\
+	display.$O\
+	dssread.$O\
+	header.$O\
+	hinv.$O\
+	image.$O\
+	patch.$O\
+	plot.$O\
+	posn.$O\
+	prose.$O\
+	qtree.$O\
+	util.$O\
+
+HFILES=sky.h
+CFLAGS=$CFLAGS -I../map
+LDFLAGS=$LDFLAGS -L$X11/lib -lX11
+
+SHORTLIB=draw bio 9
+LIB=../map/libmap/libmap.a
+<$PLAN9/src/mkone
+
+scat.$O:	strings.c
+
+$LIB:V:
+	cd ../map/libmap; mk
+
diff --git a/src/cmd/scat/patch.c b/src/cmd/scat/patch.c
new file mode 100644
index 00000000..27d04edb
--- /dev/null
+++ b/src/cmd/scat/patch.c
@@ -0,0 +1,101 @@
+#include <u.h>
+#include <libc.h>
+#include	<bio.h>
+#include "sky.h"
+
+/*
+ * dec varies from -89 to 89, inclusive.
+ * ra varies depending on dec; each patch is about 1 square degree.
+ *
+ * Northern hemisphere (0<=dec<=89):
+ *	from  0<=dec<=59, ra is every 4m,  360 values
+ *	from 60<=dec<=69, ra is every 8m,  180 values
+ *	from 70<=dec<=79, ra is every 12m, 120 values
+ *	from 80<=dec<=84, ra is every 24m,  60 values
+ *	at dec=85 and 86, ra is every 48m,  30 values
+ *	at dec=87,        ra is every 60m,  24 values
+ *	at dec=88,        ra is every 120m, 12 values
+ *	at dec=89,        ra is 12h,	     1 value
+ *
+ * Total number of patches in northern hemisphere is therefore:
+ * 	360*60+180*10+120*10+60*5+30*2+24*1+12*1+1 = 24997
+ * Total number of patches is therefore
+ *	2*24997-360 = 49634	(dec=0 has been counted twice)
+ * (cf. 41253 square degrees in the sky)
+ */
+
+void
+radec(int p, int *rah, int *ram, int *deg)
+{
+	*deg = (p&255)-90;
+	p >>= 8;
+	*rah = p/15;
+	*ram = (p%15)*4;
+	if(*deg<0)
+		(*deg)++;
+}
+
+long
+patcha(Angle ra, Angle dec)
+{
+	ra = DEG(ra);
+	dec = DEG(dec);
+	if(dec >= 0)
+		return patch(floor(ra/15), ((int)floor(ra*4))%60, floor(dec));
+	dec = -dec;
+	return patch(floor(ra/15), ((int)floor(ra*4))%60, -floor(dec));
+}
+
+#define round scatround
+
+char round[91]={	/* extra 0 is to offset the array */
+	/*  0 */    0,	 1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
+	/* 10 */	 1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
+	/* 20 */	 1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
+	/* 30 */	 1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
+	/* 40 */	 1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
+	/* 50 */	 1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
+	/* 60 */	 2,  2,  2,  2,  2,  2,  2,  2,  2,  2,
+	/* 70 */	 3,  3,  3,  3,  3,  3,  3,  3,  3,  3,
+	/* 80 */	 6,  6,  6,  6,  6, 12, 12, 15, 30, -1,
+	/* 90 */
+};
+
+long
+patch(int rah, int ram, int deg)
+{
+	int ra, dec;
+
+	/*
+	 * patches go from lower limit <= position < upper limit.
+	 * therefore dec ' " can be ignored; always inc. dec degrees.
+	 * the computed angle is then the upper limit (ignoring sign).
+	 * when done, +ve values are shifted down so 90 (0 degrees) is a value;
+	 */
+	if(rah<0 || rah>=24 || ram<0 || abs(deg)>=90){
+		fprint(2, "scat: patch: bad ra or dec %dh%dm %d\n", rah, ram, deg);
+		abort();
+	}
+	if(deg < 0)
+		deg--;
+	else if(deg < 90)
+		deg++;
+	dec = deg+90;
+	deg = abs(deg);
+	if(deg<1 || deg>90){
+		fprint(2, "scat: patch: panic %dh%dm %d\n", rah, ram, deg);
+		abort();
+	}
+	if(deg == 90)
+		ra = 180;
+	else{
+		ra = 15*rah+ram/4;
+		ra -= ra%round[deg];
+	}
+	/* close the hole at 0 */
+	if(dec > 90)
+		--dec;
+	if(ra >= 360)
+		ra -= 360;
+	return (ra<<8)|dec;
+}
diff --git a/src/cmd/scat/plate.h b/src/cmd/scat/plate.h
new file mode 100644
index 00000000..4d10bb5b
--- /dev/null
+++ b/src/cmd/scat/plate.h
@@ -0,0 +1,138 @@
+#define	RAD(x)	((x)*PI_180)
+#define	DEG(x)	((x)/PI_180)
+#define ARCSECONDS_PER_RADIAN	(DEG(1)*3600)
+#define input_nybble(infile)    input_nbits(infile,4)
+
+typedef float	Angle;	/* in radians */
+
+enum
+{
+	/*
+	 * parameters for plate
+	 */
+	Pppo1	= 0,
+	Pppo2,
+	Pppo3,
+	Pppo4,
+	Pppo5,
+	Pppo6,
+	Pamdx1,
+	Pamdx2,
+	Pamdx3,
+	Pamdx4,
+	Pamdx5,
+	Pamdx6,
+	Pamdx7,
+	Pamdx8,
+	Pamdx9,
+	Pamdx10,
+	Pamdx11,
+	Pamdx12,
+	Pamdx13,
+	Pamdx14,
+	Pamdx15,
+	Pamdx16,
+	Pamdx17,
+	Pamdx18,
+	Pamdx19,
+	Pamdx20,
+	Pamdy1,
+	Pamdy2,
+	Pamdy3,
+	Pamdy4,
+	Pamdy5,
+	Pamdy6,
+	Pamdy7,
+	Pamdy8,
+	Pamdy9,
+	Pamdy10,
+	Pamdy11,
+	Pamdy12,
+	Pamdy13,
+	Pamdy14,
+	Pamdy15,
+	Pamdy16,
+	Pamdy17,
+	Pamdy18,
+	Pamdy19,
+	Pamdy20,
+	Ppltscale,
+	Pxpixelsz,
+	Pypixelsz,
+	Ppltra,
+	Ppltrah,
+	Ppltram,
+	Ppltras,
+	Ppltdec,
+	Ppltdecd,
+	Ppltdecm,
+	Ppltdecs,
+	Pnparam,
+};
+
+typedef	struct	Plate	Plate;
+struct	Plate
+{
+	char	rgn[7];
+	char	disk;
+	Angle	ra;
+	Angle	dec;
+};
+
+typedef	struct	Header	Header;
+struct	Header
+{
+	float	param[Pnparam];
+	int	amdflag;
+
+	float	x;
+	float	y;
+	float	xi;
+	float	eta;
+};
+typedef	long	Type;
+
+typedef struct	Image	Image;
+struct	Image
+{
+	int	nx;
+	int	ny;	/* ny is the fast-varying dimension */
+	Type	a[1];
+};
+
+int	nplate;
+Plate	plate[2000];		/* needs to go to 2000 when the north comes */
+double	PI_180;
+double	TWOPI;
+int	debug;
+struct
+{
+	float	min;
+	float	max;
+	float	del;
+	double	gamma;
+	int	neg;
+} gam;
+
+char*	hms(Angle);
+char*	dms(Angle);
+double	xsqrt(double);
+Angle	dist(Angle, Angle, Angle, Angle);
+Header*	getheader(char*);
+char*	getword(char*, char*);
+void	amdinv(Header*, Angle, Angle, float, float);
+void	ppoinv(Header*, Angle, Angle);
+void	xypos(Header*, Angle, Angle, float, float);
+void	traneqstd(Header*, Angle, Angle);
+Angle	getra(char*);
+Angle	getdec(char*);
+void	getplates(void);
+
+Image*	dssread(char*);
+void	hinv(Type*, int, int);
+int	input_bit(Biobuf*);
+int	input_nbits(Biobuf*, int);
+void	qtree_decode(Biobuf*, Type*, int, int, int, int);
+void	start_inputing_bits(void);
+Bitmap*	image(Angle, Angle, Angle, Angle);
+int	dogamma(int);
diff --git a/src/cmd/scat/plot.c b/src/cmd/scat/plot.c
new file mode 100644
index 00000000..e87985ff
--- /dev/null
+++ b/src/cmd/scat/plot.c
@@ -0,0 +1,952 @@
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#include <draw.h>
+#include <event.h>
+#include <ctype.h>
+#include "map.h"
+#undef	RAD
+#undef	TWOPI
+#include "sky.h"
+
+	/* convert to milliarcsec */
+static long	c = MILLIARCSEC;	/* 1 degree */
+static long	m5 = 1250*60*60;	/* 5 minutes of ra */
+
+DAngle	ramin;
+DAngle	ramax;
+DAngle	decmin;
+DAngle	decmax;
+int		folded;
+
+Image	*grey;
+Image	*alphagrey;
+Image	*green;
+Image	*lightblue;
+Image	*lightgrey;
+Image	*ocstipple;
+Image	*suncolor;
+Image	*mooncolor;
+Image	*shadowcolor;
+Image	*mercurycolor;
+Image	*venuscolor;
+Image	*marscolor;
+Image	*jupitercolor;
+Image	*saturncolor;
+Image	*uranuscolor;
+Image	*neptunecolor;
+Image	*plutocolor;
+Image	*cometcolor;
+
+Planetrec	*planet;
+
+long	mapx0, mapy0;
+long	mapra, mapdec;
+double	mylat, mylon, mysid;
+double	mapscale;
+double	maps;
+int (*projection)(struct place*, double*, double*);
+
+char *fontname = "/lib/font/bit/lucida/unicode.6.font";
+
+/* types Coord and Loc correspond to types in map(3) thus:
+   Coord == struct coord;
+   Loc == struct place, except longitudes are measured
+     positive east for Loc and west for struct place
+*/
+
+typedef struct Xyz Xyz;
+typedef struct coord Coord;
+typedef struct Loc Loc;
+
+struct Xyz{
+	double x,y,z;
+};
+
+struct Loc{
+	Coord nlat, elon;
+};
+
+Xyz north = { 0, 0, 1 };
+
+int
+plotopen(void)
+{
+	if(display != nil)
+		return 1;
+	if(initdraw(drawerror, nil, nil) < 0){
+		fprint(2, "initdisplay failed: %r\n");
+		return -1;
+	}
+	grey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x777777FF);
+	lightgrey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0xAAAAAAFF);
+	alphagrey = allocimage(display, Rect(0, 0, 2, 2), RGBA32, 1, 0x77777777);
+	green = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x00AA00FF);
+	lightblue = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x009EEEFF);
+	ocstipple = allocimage(display, Rect(0, 0, 2, 2), CMAP8, 1, 0xAAAAAAFF);
+	draw(ocstipple, Rect(0, 0, 1, 1), display->black, nil, ZP);
+	draw(ocstipple, Rect(1, 1, 2, 2), display->black, nil, ZP);
+
+	suncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFF77FF);
+	mooncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAAAFF);
+	shadowcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x00000055);
+	mercurycolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFAAAAFF);
+	venuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);
+	marscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFF5555FF);
+	jupitercolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFFAAFF);
+	saturncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAFFAAFF);
+	uranuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77DDDDFF);
+	neptunecolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77FF77FF);
+	plutocolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x7777FFFF);
+	cometcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);
+	font = openfont(display, fontname);
+	if(font == nil)
+		fprint(2, "warning: no font %s: %r\n", fontname);
+	return 1;
+}
+
+/*
+ * Function heavens() for setting up observer-eye-view
+ * sky charts, plus subsidiary functions.
+ * Written by Doug McIlroy.
+ */
+
+/* heavens(nlato, wlono, nlatc, wlonc) sets up a map(3)-style
+   coordinate change (by calling orient()) and returns a
+   projection function heavens for calculating a star map
+   centered on nlatc,wlonc and rotated so it will appear
+   upright as seen by a ground observer at nlato,wlono.
+   all coordinates (north latitude and west longitude)
+   are in degrees.
+*/
+
+Coord
+mkcoord(double degrees)
+{
+	Coord c;
+
+	c.l = degrees*PI/180;
+	c.s = sin(c.l);
+	c.c = cos(c.l);
+	return c;
+}
+
+Xyz
+ptov(struct place p)
+{
+	Xyz v;
+
+	v.z = p.nlat.s;
+	v.x = p.nlat.c * p.wlon.c;
+	v.y = -p.nlat.c * p.wlon.s;
+	return v;
+}
+
+double
+dot(Xyz a, Xyz b)
+{
+	return a.x*b.x + a.y*b.y + a.z*b.z;
+}
+
+Xyz
+cross(Xyz a, Xyz b)
+{
+	Xyz v;
+
+	v.x = a.y*b.z - a.z*b.y;
+	v.y = a.z*b.x - a.x*b.z;
+	v.z = a.x*b.y - a.y*b.x;
+	return v;
+}
+
+double
+len(Xyz a)
+{
+	return sqrt(dot(a, a));
+}
+
+/* An azimuth vector from a to b is a tangent to
+   the sphere at a pointing along the (shorter)
+   great-circle direction to b.  It lies in the
+   plane of the vectors a and b, is perpendicular
+   to a, and has a positive compoent along b,
+   provided a and b span a 2-space.  Otherwise
+   it is null.  (aXb)Xa, where X denotes cross
+   product, is such a vector. */
+
+Xyz
+azvec(Xyz a, Xyz b)
+{
+	return cross(cross(a,b), a);
+}
+
+/* Find the azimuth of point b as seen
+   from point a, given that a is distinct
+   from either pole and from b */
+
+double
+azimuth(Xyz a, Xyz b)
+{
+	Xyz aton = azvec(a,north);
+	Xyz atob = azvec(a,b);
+	double lenaton = len(aton);
+	double lenatob = len(atob);
+	double az = acos(dot(aton,atob)/(lenaton*lenatob));
+
+	if(dot(b,cross(north,a)) < 0)
+		az = -az;
+	return az;
+}
+
+/* Find the rotation (3rd argument of orient() in the
+   map projection package) for the map described
+   below.  The return is radians; it must be converted
+   to degrees for orient().
+*/
+
+double
+rot(struct place center, struct place zenith)
+{
+	Xyz cen = ptov(center);
+	Xyz zen = ptov(zenith);
+
+	if(cen.z > 1-FUZZ) 	/* center at N pole */
+		return PI + zenith.wlon.l;
+	if(cen.z < FUZZ-1)		/* at S pole */
+		return -zenith.wlon.l;
+	if(fabs(dot(cen,zen)) > 1-FUZZ)	/* at zenith */
+		return 0;
+	return azimuth(cen, zen);
+}
+
+int (*
+heavens(double zlatdeg, double zlondeg, double clatdeg, double clondeg))(struct place*, double*, double*)
+{
+	struct place center;
+	struct place zenith;
+
+	center.nlat = mkcoord(clatdeg);
+	center.wlon = mkcoord(clondeg);
+	zenith.nlat = mkcoord(zlatdeg);
+	zenith.wlon = mkcoord(zlondeg);
+	projection = stereographic();
+	orient(clatdeg, clondeg, rot(center, zenith)*180/PI);
+	return projection;
+}
+
+int
+maptoxy(long ra, long dec, double *x, double *y)
+{
+	double lat, lon;
+	struct place pl;
+
+	lat = angle(dec);
+	lon = angle(ra);
+	pl.nlat.l = lat;
+	pl.nlat.s = sin(lat);
+	pl.nlat.c = cos(lat);
+	pl.wlon.l = lon;
+	pl.wlon.s = sin(lon);
+	pl.wlon.c = cos(lon);
+	normalize(&pl);
+	return projection(&pl, x, y);
+}
+
+/* end of 'heavens' section */
+
+int
+setmap(long ramin, long ramax, long decmin, long decmax, Rectangle r, int zenithup)
+{
+	int c;
+	double minx, maxx, miny, maxy;
+
+	c = 1000*60*60;
+	mapra = ramax/2+ramin/2;
+	mapdec = decmax/2+decmin/2;
+
+	if(zenithup)
+		projection = heavens(mylat, mysid, (double)mapdec/c, (double)mapra/c);
+	else{
+		orient((double)mapdec/c, (double)mapra/c, 0.);
+		projection = stereographic();
+	}
+	mapx0 = (r.max.x+r.min.x)/2;
+	mapy0 = (r.max.y+r.min.y)/2;
+	maps = ((double)Dy(r))/(double)(decmax-decmin);
+	if(maptoxy(ramin, decmin, &minx, &miny) < 0)
+		return -1;
+	if(maptoxy(ramax, decmax, &maxx, &maxy) < 0)
+		return -1;
+	/*
+	 * It's tricky to get the scale right.  This noble attempt scales
+	 * to fit the lengths of the sides of the rectangle.
+	 */
+	mapscale = Dx(r)/(minx-maxx);
+	if(mapscale > Dy(r)/(maxy-miny))
+		mapscale = Dy(r)/(maxy-miny);
+	/*
+	 * near the pole It's not a rectangle, though, so this scales
+	 * to fit the corners of the trapezoid (a triangle at the pole).
+	 */
+	mapscale *= (cos(angle(mapdec))+1.)/2;
+	if(maxy  < miny){
+		/* upside down, between zenith and pole */
+		fprint(2, "reverse plot\n");
+		mapscale = -mapscale;
+	}
+	return 1;
+}
+
+Point
+map(long ra, long dec)
+{
+	double x, y;
+	Point p;
+
+	if(maptoxy(ra, dec, &x, &y) > 0){
+		p.x = mapscale*x + mapx0;
+		p.y = mapscale*-y + mapy0;
+	}else{
+		p.x = -100;
+		p.y = -100;
+	}
+	return p;
+}
+
+int
+dsize(int mag)	/* mag is 10*magnitude; return disc size */
+{
+	double d;
+
+	mag += 25;	/* make mags all positive; sirius is -1.6m */
+	d = (130-mag)/10;
+	/* if plate scale is huge, adjust */
+	if(maps < 100.0/MILLIARCSEC)
+		d *= .71;
+	if(maps < 50.0/MILLIARCSEC)
+		d *= .71;
+	return d;
+}
+
+void
+drawname(Image *scr, Image *col, char *s, int ra, int dec)
+{
+	Point p;
+
+	if(font == nil)
+		return;
+	p = addpt(map(ra, dec), Pt(4, -1));	/* font has huge ascent */
+	string(scr, p, col, ZP, font, s);
+}
+
+int
+npixels(DAngle diam)
+{
+	Point p0, p1;
+
+	diam = DEG(angle(diam)*MILLIARCSEC);	/* diam in milliarcsec */
+	diam /= 2;	/* radius */
+	/* convert to plate scale */
+	/* BUG: need +100 because we sometimes crash in library if we plot the exact center */
+	p0 = map(mapra+100, mapdec);
+	p1 = map(mapra+100+diam, mapdec);
+	return hypot(p0.x-p1.x, p0.y-p1.y);
+}
+
+void
+drawdisc(Image *scr, DAngle semidiam, int ring, Image *color, Point pt, char *s)
+{
+	int d;
+
+	d = npixels(semidiam*2);
+	if(d < 5)
+		d = 5;
+	fillellipse(scr, pt, d, d, color, ZP);
+	if(ring == 1)
+		ellipse(scr, pt, 5*d/2, d/2, 0, color, ZP);
+	if(ring == 2)
+		ellipse(scr, pt, d, d, 0, grey, ZP);
+	if(s){
+		d = stringwidth(font, s);
+		pt.x -= d/2;
+		pt.y -= font->height/2 + 1;
+		string(scr, pt, display->black, ZP, font, s);
+	}
+}
+
+void
+drawplanet(Image *scr, Planetrec *p, Point pt)
+{
+	if(strcmp(p->name, "sun") == 0){
+		drawdisc(scr, p->semidiam, 0, suncolor, pt, nil);
+		return;
+	}
+	if(strcmp(p->name, "moon") == 0){
+		drawdisc(scr, p->semidiam, 0, mooncolor, pt, nil);
+		return;
+	}
+	if(strcmp(p->name, "shadow") == 0){
+		drawdisc(scr, p->semidiam, 2, shadowcolor, pt, nil);
+		return;
+	}
+	if(strcmp(p->name, "mercury") == 0){
+		drawdisc(scr, p->semidiam, 0, mercurycolor, pt, "m");
+		return;
+	}
+	if(strcmp(p->name, "venus") == 0){
+		drawdisc(scr, p->semidiam, 0, venuscolor, pt, "v");
+		return;
+	}
+	if(strcmp(p->name, "mars") == 0){
+		drawdisc(scr, p->semidiam, 0, marscolor, pt, "M");
+		return;
+	}
+	if(strcmp(p->name, "jupiter") == 0){
+		drawdisc(scr, p->semidiam, 0, jupitercolor, pt, "J");
+		return;
+	}
+	if(strcmp(p->name, "saturn") == 0){
+		drawdisc(scr, p->semidiam, 1, saturncolor, pt, "S");
+		
+		return;
+	}
+	if(strcmp(p->name, "uranus") == 0){
+		drawdisc(scr, p->semidiam, 0, uranuscolor, pt, "U");
+		
+		return;
+	}
+	if(strcmp(p->name, "neptune") == 0){
+		drawdisc(scr, p->semidiam, 0, neptunecolor, pt, "N");
+		
+		return;
+	}
+	if(strcmp(p->name, "pluto") == 0){
+		drawdisc(scr, p->semidiam, 0, plutocolor, pt, "P");
+		
+		return;
+	}
+	if(strcmp(p->name, "comet") == 0){
+		drawdisc(scr, p->semidiam, 0, cometcolor, pt, "C");
+		return;
+	}
+
+	pt.x -= stringwidth(font, p->name)/2;
+	pt.y -= font->height/2;
+	string(scr, pt, grey, ZP, font, p->name);
+}
+
+void
+tolast(char *name)
+{
+	int i, nlast;
+	Record *r, rr;
+
+	/* stop early to simplify inner loop adjustment */
+	nlast = 0;
+	for(i=0,r=rec; i<nrec-nlast; i++,r++)
+		if(r->type == Planet)
+			if(name==nil || strcmp(r->planet.name, name)==0){
+				rr = *r;
+				memmove(rec+i, rec+i+1, (nrec-i-1)*sizeof(Record));
+				rec[nrec-1] = rr;
+				--i;
+				--r;
+				nlast++;
+			}
+}
+
+int
+bbox(long extrara, long extradec, int quantize)
+{
+	int i;
+	Record *r;
+	int ra, dec;
+	int rah, ram, d1, d2;
+	double r0;
+
+	ramin = 0x7FFFFFFF;
+	ramax = -0x7FFFFFFF;
+	decmin = 0x7FFFFFFF;
+	decmax = -0x7FFFFFFF;
+
+	for(i=0,r=rec; i<nrec; i++,r++){
+		if(r->type == Patch){
+			radec(r->index, &rah, &ram, &dec);
+			ra = 15*rah+ram/4;
+			r0 = c/cos(dec*PI/180);
+			ra *= c;
+			dec *= c;
+			if(dec == 0)
+				d1 = c, d2 = c;
+			else if(dec < 0)
+				d1 = c, d2 = 0;
+			else
+				d1 = 0, d2 = c;
+		}else if(r->type==SAO || r->type==NGC || r->type==Planet || r->type==Abell){
+			ra = r->ngc.ra;
+			dec = r->ngc.dec;
+			d1 = 0, d2 = 0, r0 = 0;
+		}else
+			continue;
+		if(dec+d2+extradec > decmax)
+			decmax = dec+d2+extradec;
+		if(dec-d1-extradec < decmin)
+			decmin = dec-d1-extradec;
+		if(folded){
+			ra -= 180*c;
+			if(ra < 0)
+				ra += 360*c;
+		}
+		if(ra+r0+extrara > ramax)
+			ramax = ra+r0+extrara;
+		if(ra-extrara < ramin)
+			ramin = ra-extrara;
+	}
+
+	if(decmax > 90*c)
+		decmax = 90*c;
+	if(decmin < -90*c)
+		decmin = -90*c;
+	if(ramin < 0)
+		ramin += 360*c;
+	if(ramax >= 360*c)
+		ramax -= 360*c;
+
+	if(quantize){
+		/* quantize to degree boundaries */
+		ramin -= ramin%m5;
+		if(ramax%m5 != 0)
+			ramax += m5-(ramax%m5);
+		if(decmin > 0)
+			decmin -= decmin%c;
+		else
+			decmin -= c - (-decmin)%c;
+		if(decmax > 0){
+			if(decmax%c != 0)
+				decmax += c-(decmax%c);
+		}else if(decmax < 0){
+			if(decmax%c != 0)
+				decmax += ((-decmax)%c);
+		}
+	}
+
+	if(folded){
+		if(ramax-ramin > 270*c){
+			fprint(2, "ra range too wide %ld°\n", (ramax-ramin)/c);
+			return -1;
+		}
+	}else if(ramax-ramin > 270*c){
+		folded = 1;
+		return bbox(0, 0, quantize);
+	}
+
+	return 1;
+}
+
+int
+inbbox(DAngle ra, DAngle dec)
+{
+	int min;
+
+	if(dec < decmin)
+		return 0;
+	if(dec > decmax)
+		return 0;
+	min = ramin;
+	if(ramin > ramax){	/* straddling 0h0m */
+		min -= 360*c;
+		if(ra > 180*c)
+			ra -= 360*c;
+	}
+	if(ra < min)
+		return 0;
+	if(ra > ramax)
+		return 0;
+	return 1;
+}
+
+int
+gridra(long mapdec)
+{
+	mapdec = abs(mapdec)+c/2;
+	mapdec /= c;
+	if(mapdec < 60)
+		return m5;
+	if(mapdec < 80)
+		return 2*m5;
+	if(mapdec < 85)
+		return 4*m5;
+	return 8*m5;
+}
+
+#define	GREY	(nogrey? display->white : grey)
+
+void
+plot(char *flags)
+{
+	int i, j, k;
+	char *t;
+	long x, y;
+	int ra, dec;
+	int m;
+	Point p, pts[10];
+	Record *r;
+	Rectangle rect, r1;
+	int dx, dy, nogrid, textlevel, nogrey, zenithup;
+	Image *scr;
+	char *name, buf[32];
+	double v;
+
+	if(plotopen() < 0)
+		return;
+	nogrid = 0;
+	nogrey = 0;
+	textlevel = 1;
+	dx = 512;
+	dy = 512;
+	zenithup = 0;
+	for(;;){
+		if(t = alpha(flags, "nogrid")){
+			nogrid = 1;
+			flags = t;
+			continue;
+		}
+		if((t = alpha(flags, "zenith")) || (t = alpha(flags, "zenithup")) ){
+			zenithup = 1;
+			flags = t;
+			continue;
+		}
+		if((t = alpha(flags, "notext")) || (t = alpha(flags, "nolabel")) ){
+			textlevel = 0;
+			flags = t;
+			continue;
+		}
+		if((t = alpha(flags, "alltext")) || (t = alpha(flags, "alllabel")) ){
+			textlevel = 2;
+			flags = t;
+			continue;
+		}
+		if(t = alpha(flags, "dx")){
+			dx = strtol(t, &t, 0);
+			if(dx < 100){
+				fprint(2, "dx %d too small (min 100) in plot\n", dx);
+				return;
+			}
+			flags = skipbl(t);
+			continue;
+		}
+		if(t = alpha(flags, "dy")){
+			dy = strtol(t, &t, 0);
+			if(dy < 100){
+				fprint(2, "dy %d too small (min 100) in plot\n", dy);
+				return;
+			}
+			flags = skipbl(t);
+			continue;
+		}
+		if((t = alpha(flags, "nogrey")) || (t = alpha(flags, "nogray"))){
+			nogrey = 1;
+			flags = skipbl(t);
+			continue;
+		}
+		if(*flags){
+			fprint(2, "syntax error in plot\n");
+			return;
+		}
+		break;
+	}
+	flatten();
+	folded = 0;
+
+	if(bbox(0, 0, 1) < 0)
+		return;
+	if(ramax-ramin<100 || decmax-decmin<100){
+		fprint(2, "plot too small\n");
+		return;
+	}
+
+	scr = allocimage(display, Rect(0, 0, dx, dy), RGB24, 0, DBlack);
+	if(scr == nil){
+		fprint(2, "can't allocate image: %r\n");
+		return;
+	}
+	rect = scr->r;
+	rect.min.x += 16;
+	rect = insetrect(rect, 40);
+	if(setmap(ramin, ramax, decmin, decmax, rect, zenithup) < 0){
+		fprint(2, "can't set up map coordinates\n");
+		return;
+	}
+	if(!nogrid){
+		for(x=ramin; ; ){
+			for(j=0; j<nelem(pts); j++){
+				/* use double to avoid overflow */
+				v = (double)j / (double)(nelem(pts)-1);
+				v = decmin + v*(decmax-decmin);
+				pts[j] = map(x, v);
+			}
+			bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
+			ra = x;
+			if(folded){
+				ra -= 180*c;
+				if(ra < 0)
+					ra += 360*c;
+			}
+			p = pts[0];
+			p.x -= stringwidth(font, hm5(angle(ra)))/2;
+			string(scr, p, GREY, ZP, font, hm5(angle(ra)));
+			p = pts[nelem(pts)-1];
+			p.x -= stringwidth(font, hm5(angle(ra)))/2;
+			p.y -= font->height;
+			string(scr, p, GREY, ZP, font, hm5(angle(ra)));
+			if(x == ramax)
+				break;
+			x += gridra(mapdec);
+			if(x > ramax)
+				x = ramax;
+		}
+		for(y=decmin; y<=decmax; y+=c){
+			for(j=0; j<nelem(pts); j++){
+				/* use double to avoid overflow */
+				v = (double)j / (double)(nelem(pts)-1);
+				v = ramin + v*(ramax-ramin);
+				pts[j] = map(v, y);
+			}
+			bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
+			p = pts[0];
+			p.x += 3;
+			p.y -= font->height/2;
+			string(scr, p, GREY, ZP, font, deg(angle(y)));
+			p = pts[nelem(pts)-1];
+			p.x -= 3+stringwidth(font, deg(angle(y)));
+			p.y -= font->height/2;
+			string(scr, p, GREY, ZP, font, deg(angle(y)));
+		}
+	}
+	/* reorder to get planets in front of stars */
+	tolast(nil);
+	tolast("moon");		/* moon is in front of everything... */
+	tolast("shadow");	/* ... except the shadow */
+
+	for(i=0,r=rec; i<nrec; i++,r++){
+		dec = r->ngc.dec;
+		ra = r->ngc.ra;
+		if(folded){
+			ra -= 180*c;
+			if(ra < 0)
+				ra += 360*c;
+		}
+		if(textlevel){
+			name = nameof(r);
+			if(name==nil && textlevel>1 && r->type==SAO){
+				snprint(buf, sizeof buf, "SAO%ld", r->index);
+				name = buf;
+			}
+			if(name)
+				drawname(scr, nogrey? display->white : alphagrey, name, ra, dec);
+		}
+		if(r->type == Planet){
+			drawplanet(scr, &r->planet, map(ra, dec));
+			continue;
+		}
+		if(r->type == SAO){
+			m = r->sao.mag;
+			if(m == UNKNOWNMAG)
+				m = r->sao.mpg;
+			if(m == UNKNOWNMAG)
+				continue;
+			m = dsize(m);
+			if(m < 3)
+				fillellipse(scr, map(ra, dec), m, m, nogrey? display->white : lightgrey, ZP);
+			else{
+				ellipse(scr, map(ra, dec), m+1, m+1, 0, display->black, ZP);
+				fillellipse(scr, map(ra, dec), m, m, display->white, ZP);
+			}
+			continue;
+		}
+		if(r->type == Abell){
+			ellipse(scr, addpt(map(ra, dec), Pt(-3, 2)), 2, 1, 0, lightblue, ZP);
+			ellipse(scr, addpt(map(ra, dec), Pt(3, 2)), 2, 1, 0, lightblue, ZP);
+			ellipse(scr, addpt(map(ra, dec), Pt(0, -2)), 1, 2, 0, lightblue, ZP);
+			continue;
+		}
+		switch(r->ngc.type){
+		case Galaxy:
+			j = npixels(r->ngc.diam);
+			if(j < 4)
+				j = 4;
+			if(j > 10)
+				k = j/3;
+			else
+				k = j/2;
+			ellipse(scr, map(ra, dec), j, k, 0, lightblue, ZP);
+			break;
+
+		case PlanetaryN:
+			p = map(ra, dec);
+			j = npixels(r->ngc.diam);
+			if(j < 3)
+				j = 3;
+			ellipse(scr, p, j, j, 0, green, ZP);
+			line(scr, Pt(p.x, p.y+j+1), Pt(p.x, p.y+j+4),
+				Endsquare, Endsquare, 0, green, ZP);
+			line(scr, Pt(p.x, p.y-(j+1)), Pt(p.x, p.y-(j+4)),
+				Endsquare, Endsquare, 0, green, ZP);
+			line(scr, Pt(p.x+j+1, p.y), Pt(p.x+j+4, p.y),
+				Endsquare, Endsquare, 0, green, ZP);
+			line(scr, Pt(p.x-(j+1), p.y), Pt(p.x-(j+4), p.y),
+				Endsquare, Endsquare, 0, green, ZP);
+			break;
+
+		case DiffuseN:
+		case NebularCl:
+			p = map(ra, dec);
+			j = npixels(r->ngc.diam);
+			if(j < 4)
+				j = 4;
+			r1.min = Pt(p.x-j, p.y-j);
+			r1.max = Pt(p.x+j, p.y+j);
+			if(r->ngc.type != DiffuseN)
+				draw(scr, r1, ocstipple, ocstipple, ZP);
+			line(scr, Pt(p.x-j, p.y-j), Pt(p.x+j, p.y-j),
+				Endsquare, Endsquare, 0, green, ZP);
+			line(scr, Pt(p.x-j, p.y+j), Pt(p.x+j, p.y+j),
+				Endsquare, Endsquare, 0, green, ZP);
+			line(scr, Pt(p.x-j, p.y-j), Pt(p.x-j, p.y+j),
+				Endsquare, Endsquare, 0, green, ZP);
+			line(scr, Pt(p.x+j, p.y-j), Pt(p.x+j, p.y+j),
+				Endsquare, Endsquare, 0, green, ZP);
+			break;
+
+		case OpenCl:
+			p = map(ra, dec);
+			j = npixels(r->ngc.diam);
+			if(j < 4)
+				j = 4;
+			fillellipse(scr, p, j, j, ocstipple, ZP);
+			break;
+
+		case GlobularCl:
+			j = npixels(r->ngc.diam);
+			if(j < 4)
+				j = 4;
+			p = map(ra, dec);
+			ellipse(scr, p, j, j, 0, lightgrey, ZP);
+			line(scr, Pt(p.x-(j-1), p.y), Pt(p.x+j, p.y),
+				Endsquare, Endsquare, 0, lightgrey, ZP);
+			line(scr, Pt(p.x, p.y-(j-1)), Pt(p.x, p.y+j),
+				Endsquare, Endsquare, 0, lightgrey, ZP);
+			break;
+
+		}
+	}
+	flushimage(display, 1);
+	displayimage(scr);
+}
+
+int
+runcommand(char *command, int p[2])
+{
+	switch(rfork(RFPROC|RFFDG|RFNOWAIT)){
+	case -1:
+		return -1;
+	default:
+		break;
+	case 0:
+		dup(p[1], 1);
+		close(p[0]);
+		close(p[1]);
+		execlp("rc", "rc", "-c", command, nil);
+		fprint(2, "can't exec %s: %r", command);
+		exits(nil);
+	}
+	return 1;
+}
+
+void
+parseplanet(char *line, Planetrec *p)
+{
+	char *fld[6];
+	int i, nfld;
+	char *s;
+
+	if(line[0] == '\0')
+		return;
+	line[10] = '\0';	/* terminate name */
+	s = strrchr(line, ' ');
+	if(s == nil)
+		s = line;
+	else
+		s++;
+	strcpy(p->name, s);
+	for(i=0; s[i]!='\0'; i++)
+		p->name[i] = tolower(s[i]);
+	p->name[i] = '\0';
+	nfld = getfields(line+11, fld, nelem(fld), 1, " ");
+	p->ra = dangle(getra(fld[0]));
+	p->dec = dangle(getra(fld[1]));
+	p->az = atof(fld[2])*MILLIARCSEC;
+	p->alt = atof(fld[3])*MILLIARCSEC;
+	p->semidiam = atof(fld[4])*1000;
+	if(nfld > 5)
+		p->phase = atof(fld[5]);
+	else
+		p->phase = 0;
+}
+
+void
+astro(char *flags, int initial)
+{
+	int p[2];
+	int i, n, np;
+	char cmd[256], buf[4096], *lines[20], *fld[10];
+
+	snprint(cmd, sizeof cmd, "astro -p %s", flags);
+	if(pipe(p) < 0){
+		fprint(2, "can't pipe: %r\n");
+		return;
+	}
+	if(runcommand(cmd, p) < 0){
+		close(p[0]);
+		close(p[1]);
+		fprint(2, "can't run astro: %r");
+		return;
+	}
+	close(p[1]);
+	n = readn(p[0], buf, sizeof buf-1);
+	if(n <= 0){
+		fprint(2, "no data from astro\n");
+		return;
+	}
+	if(!initial)
+		Bwrite(&bout, buf, n);
+	buf[n] = '\0';
+	np = getfields(buf, lines, nelem(lines), 0, "\n");
+	if(np <= 1){
+		fprint(2, "astro: not enough output\n");
+		return;
+	}
+	Bprint(&bout, "%s\n", lines[0]);
+	Bflush(&bout);
+	/* get latitude and longitude */
+	if(getfields(lines[0], fld, nelem(fld), 1, " ") < 8)
+		fprint(2, "astro: can't read longitude: too few fields\n");
+	else{
+		mysid = getra(fld[5])*180./PI;
+		mylat = getra(fld[6])*180./PI;
+		mylon = getra(fld[7])*180./PI;
+	}
+	/*
+	 * Each time we run astro, we generate a new planet list
+	 * so multiple appearances of a planet may exist as we plot
+	 * its motion over time.
+	 */
+	planet = malloc(NPlanet*sizeof planet[0]);
+	if(planet == nil){
+		fprint(2, "astro: malloc failed: %r\n");
+		exits("malloc");
+	}
+	memset(planet, 0, NPlanet*sizeof planet[0]);
+	for(i=1; i<np; i++)
+		parseplanet(lines[i], &planet[i-1]);
+}
diff --git a/src/cmd/scat/posn.c b/src/cmd/scat/posn.c
new file mode 100644
index 00000000..5cce14d7
--- /dev/null
+++ b/src/cmd/scat/posn.c
@@ -0,0 +1,247 @@
+#include	<u.h>
+#include	<libc.h>
+#include	<bio.h>
+#include	"sky.h"
+
+void
+amdinv(Header *h, Angle ra, Angle dec, float mag, float col)
+{
+	int i, max_iterations;
+	float tolerance;
+	float object_x, object_y, delta_x, delta_y, f, fx, fy, g, gx, gy;
+	float x1, x2, x3, x4;
+	float y1, y2, y3, y4;
+
+	/*
+	 *  Initialize
+	 */
+	max_iterations = 50;
+	tolerance = 0.0000005;
+
+	/*
+	 *  Convert RA and Dec to St.coords
+	 */
+	traneqstd(h, ra, dec);
+
+	/*
+	 *  Set initial value for x,y
+	 */
+	object_x = h->xi/h->param[Ppltscale];
+	object_y = h->eta/h->param[Ppltscale];
+
+	/*
+	 *  Iterate by Newtons method
+	 */
+	for(i = 0; i < max_iterations; i++) {
+		/*
+		 *  X plate model
+		 */
+		x1 = object_x;
+		x2 = x1 * object_x;
+		x3 = x1 * object_x;
+		x4 = x1 * object_x;
+
+		y1 = object_y;
+		y2 = y1 * object_y;
+		y3 = y1 * object_y;
+		y4 = y1 * object_y;
+
+		f =
+			h->param[Pamdx1] * x1 +
+			h->param[Pamdx2] * y1 +
+		 	h->param[Pamdx3] +
+			h->param[Pamdx4] * x2 +
+			h->param[Pamdx5] * x1*y1 +
+			h->param[Pamdx6] * y2 +
+		   	h->param[Pamdx7] * (x2+y2) +
+			h->param[Pamdx8] * x2*x1 +
+			h->param[Pamdx9] * x2*y1 +
+			h->param[Pamdx10] * x1*y2 +
+			h->param[Pamdx11] * y3 +
+			h->param[Pamdx12] * x1* (x2+y2) +
+			h->param[Pamdx13] * x1 * (x2+y2) * (x2+y2) +
+			h->param[Pamdx14] * mag +
+			h->param[Pamdx15] * mag*mag +
+			h->param[Pamdx16] * mag*mag*mag +
+			h->param[Pamdx17] * mag*x1 +
+			h->param[Pamdx18] * mag * (x2+y2) +
+			h->param[Pamdx19] * mag*x1 * (x2+y2) +
+			h->param[Pamdx20] * col;
+
+		/*
+		 *  Derivative of X model wrt x
+		 */
+		fx =
+			h->param[Pamdx1] +
+			h->param[Pamdx4] * 2*x1 +
+			h->param[Pamdx5] * y1 +
+			h->param[Pamdx7] * 2*x1 +
+			h->param[Pamdx8] * 3*x2 +
+			h->param[Pamdx9] * 2*x1*y1 +
+			h->param[Pamdx10] * y2 +
+			h->param[Pamdx12] * (3*x2+y2) +
+			h->param[Pamdx13] * (5*x4 + 6*x2*y2 + y4) +
+			h->param[Pamdx17] * mag +
+			h->param[Pamdx18] * mag*2*x1 +
+			h->param[Pamdx19] * mag*(3*x2+y2);
+
+		/*
+		 *  Derivative of X model wrt y
+		 */
+		fy =
+			h->param[Pamdx2] +
+			h->param[Pamdx5] * x1 +
+			h->param[Pamdx6] * 2*y1 +
+			h->param[Pamdx7] * 2*y1 +
+			h->param[Pamdx9] * x2 +
+			h->param[Pamdx10] * x1*2*y1 +
+			h->param[Pamdx11] * 3*y2 +
+			h->param[Pamdx12] * 2*x1*y1 +
+			h->param[Pamdx13] * 4*x1*y1* (x2+y2) +
+			h->param[Pamdx18] * mag*2*y1 +
+			h->param[Pamdx19] * mag*2*x1*y1;
+		/*
+		 *  Y plate model
+		 */
+		g =
+			h->param[Pamdy1] * y1 +
+			h->param[Pamdy2] * x1 +
+			h->param[Pamdy3] +
+			h->param[Pamdy4] * y2 +
+			h->param[Pamdy5] * y1*x1 +
+			h->param[Pamdy6] * x2 +
+			h->param[Pamdy7] * (x2+y2) +
+			h->param[Pamdy8] * y3 +
+			h->param[Pamdy9] * y2*x1 +
+			h->param[Pamdy10] * y1*x3 +
+			h->param[Pamdy11] * x3 +
+			h->param[Pamdy12] * y1 * (x2+y2) +
+			h->param[Pamdy13] * y1 * (x2+y2) * (x2+y2) +
+			h->param[Pamdy14] * mag +
+			h->param[Pamdy15] * mag*mag +
+			h->param[Pamdy16] * mag*mag*mag +
+			h->param[Pamdy17] * mag*y1 +
+			h->param[Pamdy18] * mag * (x2+y2) +
+			h->param[Pamdy19] * mag*y1 * (x2+y2) +
+			h->param[Pamdy20] * col;
+
+		/*
+		 *  Derivative of Y model wrt x
+		 */
+		gx =
+			h->param[Pamdy2] +
+			h->param[Pamdy5] * y1 +
+			h->param[Pamdy6] * 2*x1 +
+			h->param[Pamdy7] * 2*x1 +
+			h->param[Pamdy9] * y2 +
+			h->param[Pamdy10] * y1*2*x1 +
+			h->param[Pamdy11] * 3*x2 +
+			h->param[Pamdy12] * 2*x1*y1 +
+			h->param[Pamdy13] * 4*x1*y1 * (x2+y2) +
+			h->param[Pamdy18] * mag*2*x1 +
+			h->param[Pamdy19] * mag*y1*2*x1;
+
+		/*
+		 *  Derivative of Y model wrt y
+		 */
+		gy =
+			h->param[Pamdy1] +
+			h->param[Pamdy4] * 2*y1 +
+			h->param[Pamdy5] * x1 +
+			h->param[Pamdy7] * 2*y1 +
+			h->param[Pamdy8] * 3*y2 +
+			h->param[Pamdy9] * 2*y1*x1 +
+			h->param[Pamdy10] * x2 +
+			h->param[Pamdy12] * 3*y2 +
+			h->param[Pamdy13] * (5*y4 + 6*x2*y2 + x4) +
+			h->param[Pamdy17] * mag +
+			h->param[Pamdy18] * mag*2*y1 +
+			h->param[Pamdy19] * mag*(x2 + 3*y2);
+
+		f = f - h->xi;
+		g = g - h->eta;
+		delta_x = (-f*gy+g*fy) / (fx*gy-fy*gx);
+		delta_y = (-g*fx+f*gx) / (fx*gy-fy*gx);
+		object_x = object_x + delta_x;
+		object_y = object_y + delta_y;
+		if((fabs(delta_x) < tolerance) && (fabs(delta_y) < tolerance))
+			break;
+	}
+
+	/*
+	 *  Convert mm from plate center to pixels
+	 */
+	h->x = (h->param[Pppo3]-object_x*1000.0)/h->param[Pxpixelsz];
+	h->y = (h->param[Pppo6]+object_y*1000.0)/h->param[Pypixelsz];
+}
+
+void
+ppoinv(Header *h, Angle ra, Angle dec)
+{
+
+	/*
+	 *  Convert RA and Dec to standard coords.
+	 */
+	traneqstd(h, ra, dec);
+
+	/*
+	 *  Convert st.coords from arcsec to radians
+	 */
+	h->xi  /= ARCSECONDS_PER_RADIAN;
+	h->eta /= ARCSECONDS_PER_RADIAN;
+
+	/*
+	 *  Compute PDS coordinates from solution
+	 */
+	h->x =
+		h->param[Pppo1]*h->xi +
+		h->param[Pppo2]*h->eta +
+		h->param[Pppo3];
+	h->y =
+		h->param[Pppo4]*h->xi +
+		h->param[Pppo5]*h->eta +
+		h->param[Pppo6];
+	/*
+	 * Convert x,y from microns to pixels
+	 */
+	h->x /= h->param[Pxpixelsz];
+	h->y /= h->param[Pypixelsz];
+
+}
+
+void
+traneqstd(Header *h, Angle object_ra, Angle object_dec)
+{
+	float div;
+
+	/*
+	 *  Find divisor
+	 */
+	div =
+		(sin(object_dec)*sin(h->param[Ppltdec]) +
+		cos(object_dec)*cos(h->param[Ppltdec]) *
+		cos(object_ra - h->param[Ppltra]));
+
+	/*
+	 *  Compute standard coords and convert to arcsec
+	 */
+	h->xi = cos(object_dec) *
+		sin(object_ra - h->param[Ppltra]) *
+		ARCSECONDS_PER_RADIAN/div;
+
+	h->eta = (sin(object_dec)*cos(h->param[Ppltdec])-
+		cos(object_dec)*sin(h->param[Ppltdec])*
+		cos(object_ra - h->param[Ppltra]))*
+		ARCSECONDS_PER_RADIAN/div;
+
+}
+
+void
+xypos(Header *h, Angle ra, Angle dec, float mag, float col)
+{
+	if (h->amdflag) {
+		amdinv(h, ra, dec, mag, col);
+	} else {
+		ppoinv(h, ra, dec);
+	}
+}
diff --git a/src/cmd/scat/prose.c b/src/cmd/scat/prose.c
new file mode 100644
index 00000000..17a09750
--- /dev/null
+++ b/src/cmd/scat/prose.c
@@ -0,0 +1,168 @@
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#include "sky.h"
+
+extern	Biobuf	bout;
+
+char*
+append(char *p, char *s)
+{
+	while(*s)
+		*p++ = *s++;
+	return p;
+}
+
+int
+matchlen(char *a, char *b)
+{
+	int n;
+
+	for(n=0; *a==*b; a++, b++, n++)
+		if(*a == 0)
+			return n;
+	if(*a == 0)
+		return n;
+	return 0;
+}
+
+char*
+prose(char *s, char *desc[][2], short index[])
+{
+	static char buf[512];
+	char *p=buf;
+	int i, j, k, max;
+
+	j = 0;
+	while(*s){
+		if(p >= buf+sizeof buf)
+			abort();
+		if(*s == ' '){
+			if(p>buf && p[-1]!=' ')
+				*p++ = ' ';
+			s++;
+			continue;
+		}
+		if(*s == ','){
+			*p++ = ';', s++;
+			continue;
+		}
+		if(s[0]=='M' && '0'<=s[1] && s[1]<='9'){	/* Messier tag */
+			*p++ = *s++;
+			continue;	/* below will copy the number */
+		}
+		if((i=index[(uchar)*s]) == -1){
+	Dup:
+			switch(*s){
+			default:
+				while(*s && *s!=',' && *s!=' ')
+					*p++=*s++;
+				break;
+
+			case '0': case '1': case '2': case '3': case '4':
+			case '5': case '6': case '7': case '8': case '9':
+				while('0'<=*s && *s<='9')
+					*p++ = *s++;
+				if(*s=='\'' || *s=='s')
+					*p++ = *s++;
+				break;
+
+			case '(': case ')':
+			case '\'': case '"':
+			case '&': case '-': case '+':
+				*p++ = *s++;
+				break;
+
+			case '*':
+				if('0'<=s[1] && s[1]<='9'){
+					int flag=0;
+					s++;
+				Pnumber:
+					while('0'<=*s && *s<='9')
+						*p++=*s++;
+					if(s[0] == '-'){
+						*p++ = *s++;
+						flag++;
+						goto Pnumber;
+					}
+					if(s[0]==',' && s[1]==' ' && '0'<=s[2] && s[2]<='9'){
+						*p++ = *s++;
+						s++;	/* skip blank */
+						flag++;
+						goto Pnumber;
+					}
+					if(s[0] == '.'){
+						if(s[1]=='.' && s[2]=='.'){
+							*p++ = '-';
+							s += 3;
+							flag++;
+							goto Pnumber;
+						}
+						*p++ = *s++;
+						goto Pnumber;
+					}
+					p = append(p, "m star");
+					if(flag)
+						*p++ = 's';
+					*p++ = ' ';
+					break;
+				}
+				if(s[1] == '*'){
+					if(s[2] == '*'){
+						p = append(p, "triple star ");
+						s += 3;
+					}else{
+						p = append(p, "double star ");
+						s += 2;
+					}
+					break;
+				}
+				p = append(p, "star ");
+				s++;
+				break;
+			}
+			continue;
+		}
+		for(max=-1; desc[i][0] && desc[i][0][0]==*s; i++){
+			k = matchlen(desc[i][0], s);
+			if(k > max)
+				max = k, j = i;
+		}
+		if(max == 0)
+			goto Dup;
+		s += max;
+		for(k=0; desc[j][1][k]; k++)
+			*p++=desc[j][1][k];
+		if(*s == ' ')
+			*p++ = *s++;
+		else if(*s == ',')
+			*p++ = ';', s++;
+		else
+			*p++ = ' ';
+	}
+	*p = 0;
+	return buf;
+}
+
+void
+prdesc(char *s, char *desc[][2], short index[])
+{
+	int c, j;
+
+	if(index[0] == 0){
+		index[0] = 1;
+		for(c=1, j=0; c<128; c++)
+			if(desc[j][0]==0 || desc[j][0][0]>c)
+				index[c] = -1;
+			else if(desc[j][0][0] == c){
+				index[c] = j;
+				while(desc[j][0] && desc[j][0][0] == c)
+					j++;
+				if(j >= NINDEX){
+					fprint(2, "scat: internal error: too many prose entries\n");
+					exits("NINDEX");
+				}
+			}
+	}
+	Bprint(&bout, "\t%s [%s]\n", prose(s, desc, index), s);
+}
diff --git a/src/cmd/scat/qtree.c b/src/cmd/scat/qtree.c
new file mode 100644
index 00000000..a3549c65
--- /dev/null
+++ b/src/cmd/scat/qtree.c
@@ -0,0 +1,278 @@
+#include	<u.h>
+#include	<libc.h>
+#include	<bio.h>
+#include	"sky.h"
+
+static void	qtree_expand(Biobuf*, uchar*, int, int, uchar*);
+static void	qtree_copy(uchar*, int, int, uchar*, int);
+static void	qtree_bitins(uchar*, int, int, Pix*, int, int);
+static void	read_bdirect(Biobuf*, Pix*, int, int, int, uchar*, int);
+
+void
+qtree_decode(Biobuf *infile, Pix *a, int n, int nqx, int nqy, int nbitplanes)
+{
+	int log2n, k, bit, b, nqmax;
+	int nx,ny,nfx,nfy,c;
+	int nqx2, nqy2;
+	unsigned char *scratch;
+
+	/*
+	 * log2n is log2 of max(nqx,nqy) rounded up to next power of 2
+	 */
+	nqmax = nqy;
+	if(nqx > nqmax)
+		nqmax = nqx;
+	log2n = log(nqmax)/LN2+0.5;
+	if (nqmax > (1<<log2n))
+		log2n++;
+
+	/*
+	 * allocate scratch array for working space
+	 */
+	nqx2 = (nqx+1)/2;
+	nqy2 = (nqy+1)/2;
+	scratch = (uchar*)malloc(nqx2*nqy2);
+	if(scratch == nil) {
+		fprint(2, "qtree_decode: insufficient memory\n");
+		exits("memory");
+	}
+
+	/*
+	 * now decode each bit plane, starting at the top
+	 * A is assumed to be initialized to zero
+	 */
+	for(bit = nbitplanes-1; bit >= 0; bit--) {
+		/*
+		 * Was bitplane was quadtree-coded or written directly?
+		 */
+		b = input_nybble(infile);
+		if(b == 0) {
+			/*
+			 * bit map was written directly
+			 */
+			read_bdirect(infile, a, n, nqx, nqy, scratch, bit);
+		} else
+		if(b != 0xf) {
+			fprint(2, "qtree_decode: bad format code %x\n",b);
+			exits("format");
+		} else {
+			/*
+			 * bitmap was quadtree-coded, do log2n expansions
+			 * read first code
+			 */
+
+			scratch[0] = input_huffman(infile);
+
+			/*
+			 * do log2n expansions, reading codes from file as necessary
+			 */
+			nx = 1;
+			ny = 1;
+			nfx = nqx;
+			nfy = nqy;
+			c = 1<<log2n;
+			for(k = 1; k<log2n; k++) {
+				/*
+				 * this somewhat cryptic code generates the sequence
+				 * n[k-1] = (n[k]+1)/2 where n[log2n]=nqx or nqy
+				 */
+				c = c>>1;
+				nx = nx<<1;
+				ny = ny<<1;
+				if(nfx <= c)
+					nx--;
+				else
+					nfx -= c;
+				if(nfy <= c)
+					ny--;
+				else
+					nfy -= c;
+				qtree_expand(infile, scratch, nx, ny, scratch);
+			}
+
+			/*
+			 * copy last set of 4-bit codes to bitplane bit of array a
+			 */
+			qtree_bitins(scratch, nqx, nqy, a, n, bit);
+		}
+	}
+	free(scratch);
+}
+
+/*
+ * do one quadtree expansion step on array a[(nqx+1)/2,(nqy+1)/2]
+ * results put into b[nqx,nqy] (which may be the same as a)
+ */
+static
+void
+qtree_expand(Biobuf *infile, uchar *a, int nx, int ny, uchar *b)
+{
+	uchar *b1;
+
+	/*
+	 * first copy a to b, expanding each 4-bit value
+	 */
+	qtree_copy(a, nx, ny, b, ny);
+
+	/*
+	 * now read new 4-bit values into b for each non-zero element
+	 */
+	b1 = &b[nx*ny];
+	while(b1 > b) {
+		b1--;
+		if(*b1 != 0)
+			*b1 = input_huffman(infile);
+	}
+}
+
+/*
+ * copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding
+ * each value to 2x2 pixels
+ * a,b may be same array
+ */
+static
+void
+qtree_copy(uchar *a, int nx, int ny, uchar *b, int n)
+{
+	int i, j, k, nx2, ny2;
+	int s00, s10;
+
+	/*
+	 * first copy 4-bit values to b
+	 * start at end in case a,b are same array
+	 */
+	nx2 = (nx+1)/2;
+	ny2 = (ny+1)/2;
+	k = ny2*(nx2-1) + ny2-1;	/* k   is index of a[i,j] */
+	for (i = nx2-1; i >= 0; i--) {
+		s00 = 2*(n*i+ny2-1);	/* s00 is index of b[2*i,2*j] */
+		for (j = ny2-1; j >= 0; j--) {
+			b[s00] = a[k];
+			k -= 1;
+			s00 -= 2;
+		}
+	}
+
+	/*
+	 * now expand each 2x2 block
+	 */
+	for(i = 0; i<nx-1; i += 2) {
+		s00 = n*i;				/* s00 is index of b[i,j] */
+		s10 = s00+n;				/* s10 is index of b[i+1,j] */
+		for(j = 0; j<ny-1; j += 2) {
+			b[s10+1] =  b[s00]     & 1;
+			b[s10  ] = (b[s00]>>1) & 1;
+			b[s00+1] = (b[s00]>>2) & 1;
+			b[s00  ] = (b[s00]>>3) & 1;
+			s00 += 2;
+			s10 += 2;
+		}
+		if(j < ny) {
+			/*
+			 * row size is odd, do last element in row
+			 * s00+1, s10+1 are off edge
+			 */
+			b[s10  ] = (b[s00]>>1) & 1;
+			b[s00  ] = (b[s00]>>3) & 1;
+		}
+	}
+	if(i < nx) {
+		/*
+		 * column size is odd, do last row
+		 * s10, s10+1 are off edge
+		 */
+		s00 = n*i;
+		for (j = 0; j<ny-1; j += 2) {
+			b[s00+1] = (b[s00]>>2) & 1;
+			b[s00  ] = (b[s00]>>3) & 1;
+			s00 += 2;
+		}
+		if(j < ny) {
+			/*
+			 * both row and column size are odd, do corner element
+			 * s00+1, s10, s10+1 are off edge
+			 */
+			b[s00  ] = (b[s00]>>3) & 1;
+		}
+	}
+}
+
+/*
+ * Copy 4-bit values from a[(nx+1)/2,(ny+1)/2] to b[nx,ny], expanding
+ * each value to 2x2 pixels and inserting into bitplane BIT of B.
+ * A,B may NOT be same array (it wouldn't make sense to be inserting
+ * bits into the same array anyway.)
+ */
+static
+void
+qtree_bitins(uchar *a, int nx, int ny, Pix *b, int n, int bit)
+{
+	int i, j;
+	Pix *s00, *s10;
+	Pix px;
+
+	/*
+	 * expand each 2x2 block
+	 */
+	for(i=0; i<nx-1; i+=2) {
+		s00 = &b[n*i];			/* s00 is index of b[i,j] */
+		s10 = s00+n;			/* s10 is index of b[i+1,j] */
+		for(j=0; j<ny-1; j+=2) {
+			px = *a++;
+			s10[1] |= ( px     & 1) << bit;
+			s10[0] |= ((px>>1) & 1) << bit;
+			s00[1] |= ((px>>2) & 1) << bit;
+			s00[0] |= ((px>>3) & 1) << bit;
+			s00 += 2;
+			s10 += 2;
+		}
+		if(j < ny) {
+			/*
+			 * row size is odd, do last element in row
+			 * s00+1, s10+1 are off edge
+			 */
+			px = *a++;
+			s10[0] |= ((px>>1) & 1) << bit;
+			s00[0] |= ((px>>3) & 1) << bit;
+		}
+	}
+	if(i < nx) {
+		/*
+		 * column size is odd, do last row
+		 * s10, s10+1 are off edge
+		 */
+		s00 = &b[n*i];
+		for(j=0; j<ny-1; j+=2) {
+			px = *a++;
+			s00[1] |= ((px>>2) & 1) << bit;
+			s00[0] |= ((px>>3) & 1) << bit;
+			s00 += 2;
+		}
+		if(j < ny) {
+			/*
+			 * both row and column size are odd, do corner element
+			 * s00+1, s10, s10+1 are off edge
+			 */
+			s00[0] |= ((*a>>3) & 1) << bit;
+		}
+	}
+}
+
+static
+void
+read_bdirect(Biobuf *infile, Pix *a, int n, int nqx, int nqy, uchar *scratch, int bit)
+{
+	int i;
+
+	/*
+	 * read bit image packed 4 pixels/nybble
+	 */
+	for(i = 0; i < ((nqx+1)/2) * ((nqy+1)/2); i++) {
+		scratch[i] = input_nybble(infile);
+	}
+
+	/*
+	 * insert in bitplane BIT of image A
+	 */
+	qtree_bitins(scratch, nqx, nqy, a, n, bit);
+}
diff --git a/src/cmd/scat/scat.c b/src/cmd/scat/scat.c
new file mode 100644
index 00000000..9579706a
--- /dev/null
+++ b/src/cmd/scat/scat.c
@@ -0,0 +1,1671 @@
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#include <draw.h>
+#include <event.h>
+#include "sky.h"
+#include "strings.c"
+
+enum
+{
+	NNGC=7840,	/* number of NGC numbers [1..NNGC] */
+	NIC = 5386,	/* number of IC numbers */
+	NNGCrec=NNGC+NIC,	/* number of records in the NGC catalog (including IC's, starting at NNGC */
+	NMrec=122,	/* number of M records */
+	NM=110,		/* number of M numbers */
+	NAbell=2712,	/* number of records in the Abell catalog */
+	NName=1000,	/* number of prose names; estimated maximum (read from editable text file) */
+	NBayer=1517,	/* number of bayer entries */
+	NSAO=258998,	/* number of SAO stars */
+	MAXcon=1932,	/* maximum number of patches in a constellation */
+	Ncon=88,	/* number of constellations */
+	Npatch=92053,	/* highest patch number */
+};
+
+char		ngctype[NNGCrec];
+Mindexrec	mindex[NMrec];
+Namerec		name[NName];
+Bayerec		bayer[NBayer];
+long		con[MAXcon];
+ushort		conindex[Ncon+1];
+long		patchaddr[Npatch+1];
+
+Record	*rec;
+Record	*orec;
+Record	*cur;
+
+char	*dir;
+int	saodb;
+int	ngcdb;
+int	abelldb;
+int	ngctypedb;
+int	mindexdb;
+int	namedb;
+int	bayerdb;
+int	condb;
+int	conindexdb;
+int	patchdb;
+char	parsed[3];
+long	nrec;
+long	nreca;
+long	norec;
+long	noreca;
+
+Biobuf	bin;
+Biobuf	bout;
+
+void
+main(int argc, char *argv[])
+{
+	char *line;
+
+	dir = unsharp(DIR);
+	Binit(&bin, 0, OREAD);
+	Binit(&bout, 1, OWRITE);
+	if(argc != 1)
+		dir = argv[1];
+	astro("", 1);
+	while(line = Brdline(&bin, '\n')){
+		line[Blinelen(&bin)-1] = 0;
+		lookup(line, 1);
+		Bflush(&bout);
+	}
+	if(display != nil){
+		closedisplay(display);
+		/* automatic refresh of rio window is triggered by mouse */
+	//	close(open("/dev/mouse", OREAD));
+	}
+	return;
+}
+
+void
+reset(void)
+{
+	nrec = 0;
+	cur = rec;
+}
+
+void
+grow(void)
+{
+	nrec++;
+	if(nreca < nrec){
+		nreca = nrec+50;
+		rec = realloc(rec, nreca*sizeof(Record));
+		if(rec == 0){
+			fprint(2, "scat: realloc fails\n");
+			exits("realloc");
+		}
+	}
+	cur = rec+nrec-1;
+}
+
+void
+copy(void)
+{
+	if(noreca < nreca){
+		noreca = nreca;
+		orec = realloc(orec, nreca*sizeof(Record));
+		if(orec == 0){
+			fprint(2, "scat: realloc fails\n");
+			exits("realloc");
+		}
+	}
+	memmove(orec, rec, nrec*sizeof(Record));
+	norec = nrec;
+}
+
+int
+eopen(char *s)
+{
+	char buf[128];
+	int f;
+
+	sprint(buf, "%s/%s.scat", dir, s);
+	f = open(buf, 0);
+	if(f<0){
+		fprint(2, "scat: can't open %s\n", buf);
+		exits("open");
+	}
+	return f;
+}
+
+
+void
+Eread(int f, char *name, void *addr, long n)
+{
+	if(read(f, addr, n) != n){	/* BUG! */
+		fprint(2, "scat: read error on %s\n", name);
+		exits("read");
+	}
+}
+
+char*
+skipbl(char *s)
+{
+	while(*s!=0 && (*s==' ' || *s=='\t'))
+		s++;
+	return s;
+}
+
+char*
+skipstr(char *s, char *t)
+{
+	while(*s && *s==*t)
+		s++, t++;
+	return skipbl(s);
+}
+
+/* produce little-endian long at address l */
+long
+Long(long *l)
+{
+	uchar *p;
+
+	p = (uchar*)l;
+	return (long)p[0]|((long)p[1]<<8)|((long)p[2]<<16)|((long)p[3]<<24);
+}
+
+/* produce little-endian long at address l */
+int
+Short(short *s)
+{
+	uchar *p;
+
+	p = (uchar*)s;
+	return p[0]|(p[1]<<8);
+}
+
+void
+nameopen(void)
+{
+	Biobuf b;
+	int i;
+	char *l, *p;
+
+	if(namedb == 0){
+		namedb = eopen("name");
+		Binit(&b, namedb, OREAD);
+		for(i=0; i<NName; i++){
+			l = Brdline(&b, '\n');
+			if(l == 0)
+				break;
+			p = strchr(l, '\t');
+			if(p == 0){
+		Badformat:
+				Bprint(&bout, "warning: name.scat bad format; line %d\n", i+1);
+				break;
+			}
+			*p++ = 0;
+			strcpy(name[i].name, l);
+			if(strncmp(p, "ngc", 3) == 0)
+				name[i].ngc = atoi(p+3);
+			else if(strncmp(p, "ic", 2) == 0)
+				name[i].ngc = atoi(p+2)+NNGC;
+			else if(strncmp(p, "sao", 3) == 0)
+				name[i].sao = atoi(p+3);
+			else if(strncmp(p, "abell", 5) == 0)
+				name[i].abell = atoi(p+5);
+			else
+				goto Badformat;
+		}
+		if(i == NName)
+			Bprint(&bout, "warning: too many names in name.scat (max %d); extra ignored\n", NName);
+		close(namedb);
+
+		bayerdb = eopen("bayer");
+		Eread(bayerdb, "bayer", bayer, sizeof bayer);
+		close(bayerdb);
+		for(i=0; i<NBayer; i++)
+			bayer[i].sao = Long(&bayer[i].sao);
+	}
+}
+
+void
+saoopen(void)
+{
+	if(saodb == 0){
+		nameopen();
+		saodb = eopen("sao");
+	}
+}
+
+void
+ngcopen(void)
+{
+	if(ngcdb == 0){
+		nameopen();
+		ngcdb = eopen("ngc2000");
+		ngctypedb = eopen("ngc2000type");
+		Eread(ngctypedb, "ngctype", ngctype, sizeof ngctype);
+		close(ngctypedb);
+	}
+}
+
+void
+abellopen(void)
+{
+	/* nothing extra to do with abell: it's directly indexed by number */
+	if(abelldb == 0)
+		abelldb = eopen("abell");
+}
+
+void
+patchopen(void)
+{
+	Biobuf *b;
+	long l, m;
+	char buf[100];
+
+	if(patchdb == 0){
+		patchdb = eopen("patch");
+		sprint(buf, "%s/patchindex.scat", dir);
+		b = Bopen(buf, OREAD);
+		if(b == 0){
+			fprint(2, "can't open %s\n", buf);
+			exits("open");
+		}
+		for(m=0,l=0; l<=Npatch; l++)
+			patchaddr[l] = m += Bgetc(b)*4;
+		Bterm(b);
+	}
+}
+
+void
+mopen(void)
+{
+	int i;
+
+	if(mindexdb == 0){
+		mindexdb = eopen("mindex");
+		Eread(mindexdb, "mindex", mindex, sizeof mindex);
+		close(mindexdb);
+		for(i=0; i<NMrec; i++)
+			mindex[i].ngc = Short(&mindex[i].ngc);
+	}
+}
+
+void
+constelopen(void)
+{
+	int i;
+
+	if(condb == 0){
+		condb = eopen("con");
+		conindexdb = eopen("conindex");
+		Eread(conindexdb, "conindex", conindex, sizeof conindex);
+		close(conindexdb);
+		for(i=0; i<Ncon+1; i++)
+			conindex[i] = Short((short*)&conindex[i]);
+	}
+}
+
+void
+lowercase(char *s)
+{
+	for(; *s; s++)
+		if('A'<=*s && *s<='Z')
+			*s += 'a'-'A';
+}
+
+int
+loadngc(long index)
+{
+	static int failed;
+	long j;
+
+	ngcopen();
+	j = (index-1)*sizeof(NGCrec);
+	grow();
+	cur->type = NGC;
+	cur->index = index;
+	seek(ngcdb, j, 0);
+	/* special case: NGC data may not be available */
+	if(read(ngcdb, &cur->ngc, sizeof(NGCrec)) != sizeof(NGCrec)){
+		if(!failed){
+			fprint(2, "scat: NGC database not available\n");
+			failed++;
+		}
+		cur->type = NONGC;
+		cur->ngc.ngc = 0;
+		cur->ngc.ra = 0;
+		cur->ngc.dec = 0;
+		cur->ngc.diam = 0;
+		cur->ngc.mag = 0;
+		return 0;
+	}
+	cur->ngc.ngc = Short(&cur->ngc.ngc);
+	cur->ngc.ra = Long(&cur->ngc.ra);
+	cur->ngc.dec = Long(&cur->ngc.dec);
+	cur->ngc.diam = Long(&cur->ngc.diam);
+	cur->ngc.mag = Short(&cur->ngc.mag);
+	return 1;
+}
+
+int
+loadabell(long index)
+{
+	long j;
+
+	abellopen();
+	j = index-1;
+	grow();
+	cur->type = Abell;
+	cur->index = index;
+	seek(abelldb, j*sizeof(Abellrec), 0);
+	Eread(abelldb, "abell", &cur->abell, sizeof(Abellrec));
+	cur->abell.abell = Short(&cur->abell.abell);
+	if(cur->abell.abell != index){
+		fprint(2, "bad format in abell catalog\n");
+		exits("abell");
+	}
+	cur->abell.ra = Long(&cur->abell.ra);
+	cur->abell.dec = Long(&cur->abell.dec);
+	cur->abell.glat = Long(&cur->abell.glat);
+	cur->abell.glong = Long(&cur->abell.glong);
+	cur->abell.rad = Long(&cur->abell.rad);
+	cur->abell.mag10 = Short(&cur->abell.mag10);
+	cur->abell.pop = Short(&cur->abell.pop);
+	cur->abell.dist = Short(&cur->abell.dist);
+	return 1;
+}
+
+int
+loadsao(int index)
+{
+	if(index<=0 || index>NSAO)
+		return 0;
+	saoopen();
+	grow();
+	cur->type = SAO;
+	cur->index = index;
+	seek(saodb, (index-1)*sizeof(SAOrec), 0);
+	Eread(saodb, "sao", &cur->sao, sizeof(SAOrec));
+	cur->sao.ra = Long(&cur->sao.ra);
+	cur->sao.dec = Long(&cur->sao.dec);
+	cur->sao.dra = Long(&cur->sao.dra);
+	cur->sao.ddec = Long(&cur->sao.ddec);
+	cur->sao.mag = Short(&cur->sao.mag);
+	cur->sao.mpg = Short(&cur->sao.mpg);
+	cur->sao.hd = Long(&cur->sao.hd);
+	return 1;
+}
+
+int
+loadplanet(int index, Record *r)
+{
+	if(index<0 || index>NPlanet || planet[index].name[0]=='\0')
+		return 0;
+	grow();
+	cur->type = Planet;
+	cur->index = index;
+	/* check whether to take new or existing record */
+	if(r == nil)
+		memmove(&cur->planet, &planet[index], sizeof(Planetrec));
+	else
+		memmove(&cur->planet, &r->planet, sizeof(Planetrec));
+	return 1;
+}
+
+int
+loadpatch(long index)
+{
+	int i;
+
+	patchopen();
+	if(index<=0 || index>Npatch)
+		return 0;
+	grow();
+	cur->type = Patch;
+	cur->index = index;
+	seek(patchdb, patchaddr[index-1], 0);
+	cur->patch.nkey = (patchaddr[index]-patchaddr[index-1])/4;
+	Eread(patchdb, "patch", cur->patch.key, cur->patch.nkey*4);
+	for(i=0; i<cur->patch.nkey; i++)
+		cur->patch.key[i] = Long(&cur->patch.key[i]);
+	return 1;
+}
+
+int
+loadtype(int t)
+{
+	int i;
+
+	ngcopen();
+	for(i=0; i<NNGCrec; i++)
+		if(t == (ngctype[i])){
+			grow();
+			cur->type = NGCN;
+			cur->index = i+1;
+		}
+	return 1;
+}
+
+void
+flatten(void)
+{
+	int i, j, notflat;
+	Record *or;
+	long key;
+
+    loop:
+	copy();
+	reset();
+	notflat = 0;
+	for(i=0,or=orec; i<norec; i++,or++){
+		switch(or->type){
+		default:
+			fprint(2, "bad type %d in flatten\n", or->type);
+			break;
+
+		case NONGC:
+			break;
+
+		case Planet:
+		case Abell:
+		case NGC:
+		case SAO:
+			grow();
+			memmove(cur, or, sizeof(Record));
+			break;
+
+		case NGCN:
+			if(loadngc(or->index))
+				notflat = 1;
+			break;
+
+		case NamedSAO:
+			loadsao(or->index);
+			notflat = 1;
+			break;
+
+		case NamedNGC:
+			if(loadngc(or->index))
+				notflat = 1;
+			break;
+
+		case NamedAbell:
+			loadabell(or->index);
+			notflat = 1;
+			break;
+
+		case PatchC:
+			loadpatch(or->index);
+			notflat = 1;
+			break;
+
+		case Patch:
+			for(j=1; j<or->patch.nkey; j++){
+				key = or->patch.key[j];
+				if((key&0x3F) == SAO)
+					loadsao((key>>8)&0xFFFFFF);
+				else if((key&0x3F) == Abell)
+					loadabell((key>>8)&0xFFFFFF);
+				else
+					loadngc((key>>16)&0xFFFF);
+			}
+			break;
+		}
+	}
+	if(notflat)
+		goto loop;
+}
+
+int
+ism(int index)
+{
+	int i;
+
+	for(i=0; i<NMrec; i++)
+		if(mindex[i].ngc == index)
+			return 1;
+	return 0;
+}
+
+char*
+alpha(char *s, char *t)
+{
+	int n;
+
+	n = strlen(t);
+	if(strncmp(s, t, n)==0 && (s[n]<'a' || 'z'<s[n]))
+		return skipbl(s+n);
+	return 0;
+	
+}
+
+char*
+text(char *s, char *t)
+{
+	int n;
+
+	n = strlen(t);
+	if(strncmp(s, t, n)==0 && (s[n]==0 || s[n]==' ' || s[n]=='\t'))
+		return skipbl(s+n);
+	return 0;
+	
+}
+
+int
+cull(char *s, int keep, int dobbox)
+{
+	int i, j, nobj, keepthis;
+	Record *or;
+	char *t;
+	int dogrtr, doless, dom, dosao, dongc, doabell;
+	int mgrtr, mless;
+	char obj[100];
+
+	memset(obj, 0, sizeof(obj));
+	nobj = 0;
+	dogrtr = 0;
+	doless = 0;
+	dom = 0;
+	dongc = 0;
+	dosao = 0;
+	doabell = 0;
+	mgrtr = mless= 0;
+	if(dobbox)
+		goto Cull;
+	for(;;){
+		if(s[0] == '>'){
+			dogrtr = 1;
+			mgrtr = 10 * strtod(s+1, &t);
+			if(mgrtr==0  && t==s+1){
+				fprint(2, "bad magnitude\n");
+				return 0;
+			}
+			s = skipbl(t);
+			continue;
+		}
+		if(s[0] == '<'){
+			doless = 1;
+			mless = 10 * strtod(s+1, &t);
+			if(mless==0  && t==s+1){
+				fprint(2, "bad magnitude\n");
+				return 0;
+			}
+			s = skipbl(t);
+			continue;
+		}
+		if(t = text(s, "m")){
+ 			dom = 1;
+			s = t;
+			continue;
+		}
+		if(t = text(s, "sao")){
+			dosao = 1;
+			s = t;
+			continue;
+		}
+		if(t = text(s, "ngc")){
+			dongc = 1;
+			s = t;
+			continue;
+		}
+		if(t = text(s, "abell")){
+			doabell = 1;
+			s = t;
+			continue;
+		}
+		for(i=0; names[i].name; i++)
+			if(t = alpha(s, names[i].name)){
+				if(nobj > 100){
+					fprint(2, "too many object types\n");
+					return 0;
+				}
+				obj[nobj++] = names[i].type;
+				s = t;
+				goto Continue;
+			}
+		break;
+	    Continue:;
+	}
+	if(*s){
+		fprint(2, "syntax error in object list\n");
+		return 0;
+	}
+
+    Cull:
+	flatten();
+	copy();
+	reset();
+	if(dom)
+		mopen();
+	if(dosao)
+		saoopen();
+	if(dongc || nobj)
+		ngcopen();
+	if(doabell)
+		abellopen();
+	for(i=0,or=orec; i<norec; i++,or++){
+		keepthis = !keep;
+		if(dobbox && inbbox(or->ngc.ra, or->ngc.dec))
+			keepthis = keep;
+		if(doless && or->ngc.mag <= mless)
+			keepthis = keep;
+		if(dogrtr && or->ngc.mag >= mgrtr)
+			keepthis = keep;
+		if(dom && (or->type==NGC && ism(or->ngc.ngc)))
+			keepthis = keep;
+		if(dongc && or->type==NGC)
+			keepthis = keep;
+		if(doabell && or->type==Abell)
+			keepthis = keep;
+		if(dosao && or->type==SAO)
+			keepthis = keep;
+		for(j=0; j<nobj; j++)
+			if(or->type==NGC && or->ngc.type==obj[j])
+				keepthis = keep;
+		if(keepthis){
+			grow();
+			memmove(cur, or, sizeof(Record));
+		}
+	}
+	return 1;
+}
+
+int
+compar(const void *va, const void *vb)
+{
+	Record *a=(Record*)va, *b=(Record*)vb;
+
+	if(a->type == b->type)
+		return a->index - b->index;
+	return a->type - b->type;
+}
+
+void
+sort(void)
+{
+	int i;
+	Record *r, *s;
+
+	if(nrec == 0)
+		return;
+	qsort(rec, nrec, sizeof(Record), compar);
+	r = rec+1;
+	s = rec;
+	for(i=1; i<nrec; i++,r++){
+		/* may have multiple instances of a planet in the scene */
+		if(r->type==s->type && r->index==s->index && r->type!=Planet)
+			continue;
+		memmove(++s, r, sizeof(Record));
+	}
+	nrec = (s+1)-rec;
+}
+
+char	greekbuf[128];
+
+char*
+togreek(char *s)
+{
+	char *t;
+	int i, n;
+	Rune r;
+
+	t = greekbuf;
+	while(*s){
+		for(i=1; i<=24; i++){
+			n = strlen(greek[i]);
+			if(strncmp(s, greek[i], n)==0 && (s[n]==' ' || s[n]=='\t')){
+				s += n;
+				t += runetochar(t, &greeklet[i]);
+				goto Cont;
+			}
+		}
+		n = chartorune(&r, s);
+		for(i=0; i<n; i++)
+			*t++ = *s++;
+    Cont:;
+	}
+	*t = 0;
+	return greekbuf;
+}
+
+char*
+fromgreek(char *s)
+{
+	char *t;
+	int i, n;
+	Rune r;
+
+	t = greekbuf;
+	while(*s){
+		n = chartorune(&r, s);
+		for(i=1; i<=24; i++){
+			if(r == greeklet[i]){
+				strcpy(t, greek[i]);
+				t += strlen(greek[i]);
+				s += n;
+				goto Cont;
+			}
+		}
+		for(i=0; i<n; i++)
+			*t++ = *s++;
+    Cont:;
+	}
+	*t = 0;
+	return greekbuf;
+}
+
+#ifdef OLD
+/*
+ * Old version
+ */
+int
+coords(int deg)
+{
+	int i;
+	int x, y;
+	Record *or;
+	long dec, ra, ndec, nra;
+	int rdeg;
+
+	flatten();
+	copy();
+	reset();
+	deg *= 2;
+	for(i=0,or=orec; i<norec; i++,or++){
+		if(or->type == Planet)	/* must keep it here */
+			loadplanet(or->index, or);
+		dec = or->ngc.dec/MILLIARCSEC;
+		ra = or->ngc.ra/MILLIARCSEC;
+		rdeg = deg/cos((dec*PI)/180);
+		for(y=-deg; y<=+deg; y++){
+			ndec = dec*2+y;
+			if(ndec/2>=90 || ndec/2<=-90)
+				continue;
+			/* fp errors hurt here, so we round 1' to the pole */
+			if(ndec >= 0)
+				ndec = ndec*500*60*60 + 60000;
+			else
+				ndec = ndec*500*60*60 - 60000;
+			for(x=-rdeg; x<=+rdeg; x++){
+				nra = ra*2+x;
+				if(nra/2 < 0)
+					nra += 360*2;
+				if(nra/2 >= 360)
+					nra -= 360*2;
+				/* fp errors hurt here, so we round up 1' */
+				nra = nra/2*MILLIARCSEC + 60000;
+				loadpatch(patcha(angle(nra), angle(ndec)));
+			}
+		}
+	}
+	sort();
+	return 1;
+}
+#endif
+
+/*
+ * New version attempts to match the boundaries of the plot better.
+ */
+int
+coords(int deg)
+{
+	int i;
+	int x, y, xx;
+	Record *or;
+	long min, circle;
+	double factor;
+
+	flatten();
+	circle = 360*MILLIARCSEC;
+	deg *= MILLIARCSEC;
+
+	/* find center */
+	folded = 0;
+	bbox(0, 0, 0);
+	/* now expand */
+	factor = cos(angle((decmax+decmin)/2));
+	if(factor < .2)
+		factor = .2;
+	factor = floor(1/factor);
+	folded = 0;
+	bbox(factor*deg, deg, 1);
+	Bprint(&bout, "%s to ", hms(angle(ramin)));
+	Bprint(&bout, "%s\n", hms(angle(ramax)));
+	Bprint(&bout, "%s to ", dms(angle(decmin)));
+	Bprint(&bout, "%s\n", dms(angle(decmax)));
+	copy();
+	reset();
+	for(i=0,or=orec; i<norec; i++,or++)
+		if(or->type == Planet)	/* must keep it here */
+			loadplanet(or->index, or);
+	min = ramin;
+	if(ramin > ramax)
+		min -= circle;
+	for(x=min; x<=ramax; x+=250*60*60){
+		xx = x;
+		if(xx < 0)
+			xx += circle;
+		for(y=decmin; y<=decmax; y+=250*60*60)
+			if(-circle/4 < y && y < circle/4)
+				loadpatch(patcha(angle(xx), angle(y)));
+	}
+	sort();
+	cull(nil, 1, 1);
+	return 1;
+}
+
+void
+pplate(char *flags)
+{
+	int i;
+	long c;
+	int na, rah, ram, d1, d2;
+	double r0;
+	int ra, dec;
+	long ramin, ramax, decmin, decmax;	/* all in degrees */
+	Record *r;
+	int folded;
+	Angle racenter, deccenter, rasize, decsize, a[4];
+	Picture *pic;
+
+	rasize = -1.0;
+	decsize = -1.0;
+	na = 0;
+	for(;;){
+		while(*flags==' ')
+			flags++;
+		if(('0'<=*flags && *flags<='9') || *flags=='+' || *flags=='-'){
+			if(na >= 3)
+				goto err;
+			a[na++] = getra(flags);
+			while(*flags && *flags!=' ')
+				flags++;
+			continue;
+		}
+		if(*flags){
+	err:
+			Bprint(&bout, "syntax error in plate\n");
+			return;
+		}
+		break;
+	}
+	switch(na){
+	case 0:
+		break;
+	case 1:
+		rasize = a[0];
+		decsize = rasize;
+		break;
+	case 2:
+		rasize = a[0];
+		decsize = a[1];
+		break;
+	case 3:
+	case 4:
+		racenter = a[0];
+		deccenter = a[1];
+		rasize = a[2];
+		if(na == 4)
+			decsize = a[3];
+		else
+			decsize = rasize;
+		if(rasize<0.0 || decsize<0.0){
+			Bprint(&bout, "negative sizes\n");
+			return;
+		}
+		goto done;
+	}
+	folded = 0;
+	/* convert to milliarcsec */
+	c = 1000*60*60;
+    Again:
+	if(nrec == 0){
+		Bprint(&bout, "empty\n");
+		return;
+	}
+	ramin = 0x7FFFFFFF;
+	ramax = -0x7FFFFFFF;
+	decmin = 0x7FFFFFFF;
+	decmax = -0x7FFFFFFF;
+	for(r=rec,i=0; i<nrec; i++,r++){
+		if(r->type == Patch){
+			radec(r->index, &rah, &ram, &dec);
+			ra = 15*rah+ram/4;
+			r0 = c/cos(RAD(dec));
+			ra *= c;
+			dec *= c;
+			if(dec == 0)
+				d1 = c, d2 = c;
+			else if(dec < 0)
+				d1 = c, d2 = 0;
+			else
+				d1 = 0, d2 = c;
+		}else if(r->type==SAO || r->type==NGC || r->type==Abell){
+			ra = r->ngc.ra;
+			dec = r->ngc.dec;
+			d1 = 0, d2 = 0, r0 = 0;
+		}else if(r->type==NGCN){
+			loadngc(r->index);
+			continue;
+		}else if(r->type==NamedSAO){
+			loadsao(r->index);
+			continue;
+		}else if(r->type==NamedNGC){
+			loadngc(r->index);
+			continue;
+		}else if(r->type==NamedAbell){
+			loadabell(r->index);
+			continue;
+		}else
+			continue;
+		if(dec+d2 > decmax)
+			decmax = dec+d2;
+		if(dec-d1 < decmin)
+			decmin = dec-d1;
+		if(folded){
+			ra -= 180*c;
+			if(ra < 0)
+				ra += 360*c;
+		}
+		if(ra+r0 > ramax)
+			ramax = ra+r0;
+		if(ra < ramin)
+			ramin = ra;
+	}
+	if(!folded && ramax-ramin>270*c){
+		folded = 1;
+		goto Again;
+	}
+	racenter = angle(ramin+(ramax-ramin)/2);
+	deccenter = angle(decmin+(decmax-decmin)/2);
+	if(rasize<0 || decsize<0){
+		rasize = angle(ramax-ramin)*cos(deccenter);
+		decsize = angle(decmax-decmin);
+	}
+    done:
+	if(DEG(rasize)>1.1 || DEG(decsize)>1.1){
+		Bprint(&bout, "plate too big: %s", ms(rasize));
+		Bprint(&bout, " x %s\n", ms(decsize));
+		Bprint(&bout, "trimming to 30'x30'\n");
+		rasize = RAD(0.5);
+		decsize = RAD(0.5);
+	}
+	Bprint(&bout, "%s %s ", hms(racenter), dms(deccenter));
+	Bprint(&bout, "%s", ms(rasize));
+	Bprint(&bout, " x %s\n", ms(decsize));
+	Bflush(&bout);
+	flatten();
+	pic = image(racenter, deccenter, rasize, decsize);
+	if(pic == 0)
+		return;
+	Bprint(&bout, "plate %s locn %d %d %d %d\n", pic->name, pic->minx, pic->miny, pic->maxx, pic->maxy);
+	Bflush(&bout);
+	displaypic(pic);
+}
+
+void
+lookup(char *s, int doreset)
+{
+	int i, j, k;
+	int rah, ram, deg;
+	char *starts, *inputline=s, *t, *u;
+	Record *r;
+	Rune c;
+	long n;
+	double x;
+	Angle ra;
+
+	lowercase(s);
+	s = skipbl(s);
+
+	if(*s == 0)
+		goto Print;
+
+	if(t = alpha(s, "flat")){
+		if(*t){
+			fprint(2, "flat takes no arguments\n");
+			return;
+		}
+		if(nrec == 0){
+			fprint(2, "no records\n");
+			return;
+		}
+		flatten();
+		goto Print;
+	}
+
+	if(t = alpha(s, "print")){
+		if(*t){
+			fprint(2, "print takes no arguments\n");
+			return;
+		}
+		for(i=0,r=rec; i<nrec; i++,r++)
+			prrec(r);
+		return;
+	}
+
+	if(t = alpha(s, "add")){
+		lookup(t, 0);
+		return;
+	}
+
+	if(t = alpha(s, "sao")){
+		n = strtoul(t, &u, 10);
+		if(n<=0 || n>NSAO)
+			goto NotFound;
+		t = skipbl(u);
+		if(*t){
+			fprint(2, "syntax error in sao\n");
+			return;
+		}
+		if(doreset)
+			reset();
+		if(!loadsao(n))
+			goto NotFound;
+		goto Print;
+	}
+
+	if(t = alpha(s, "ngc")){
+		n = strtoul(t, &u, 10);
+		if(n<=0 || n>NNGC)
+			goto NotFound;
+		t = skipbl(u);
+		if(*t){
+			fprint(2, "syntax error in ngc\n");
+			return;
+		}
+		if(doreset)
+			reset();
+		if(!loadngc(n))
+			goto NotFound;
+		goto Print;
+	}
+
+	if(t = alpha(s, "ic")){
+		n = strtoul(t, &u, 10);
+		if(n<=0 || n>NIC)
+			goto NotFound;
+		t = skipbl(u);
+		if(*t){
+			fprint(2, "syntax error in ic\n");
+			return;
+		}
+		if(doreset)
+			reset();
+		if(!loadngc(n+NNGC))
+			goto NotFound;
+		goto Print;
+	}
+
+	if(t = alpha(s, "abell")){
+		n = strtoul(t, &u, 10);
+		if(n<=0 || n>NAbell)
+			goto NotFound;
+		if(doreset)
+			reset();
+		if(!loadabell(n))
+			goto NotFound;
+		goto Print;
+	}
+
+	if(t = alpha(s, "m")){
+		n = strtoul(t, &u, 10);
+		if(n<=0 || n>NM)
+			goto NotFound;
+		mopen();
+		for(j=n-1; mindex[j].m<n; j++)
+			;
+		if(doreset)
+			reset();
+		while(mindex[j].m == n){
+			if(mindex[j].ngc){
+				grow();
+				cur->type = NGCN;
+				cur->index = mindex[j].ngc;
+			}
+			j++;
+		}
+		goto Print;
+	}
+
+	for(i=1; i<=Ncon; i++)
+		if(t = alpha(s, constel[i])){
+			if(*t){
+				fprint(2, "syntax error in constellation\n");
+				return;
+			}
+			constelopen();
+			seek(condb, 4L*conindex[i-1], 0);
+			j = conindex[i]-conindex[i-1];
+			Eread(condb, "con", con, 4*j);
+			if(doreset)
+				reset();
+			for(k=0; k<j; k++){
+				grow();
+				cur->type = PatchC;
+				cur->index = Long(&con[k]);
+			}
+			goto Print;
+		}
+
+	if(t = alpha(s, "expand")){
+		n = 0;
+		if(*t){
+			if(*t<'0' && '9'<*t){
+		Expanderr:
+				fprint(2, "syntax error in expand\n");
+				return;
+			}
+			n = strtoul(t, &u, 10);
+			t = skipbl(u);
+			if(*t)
+				goto Expanderr;
+		}
+		coords(n);
+		goto Print;
+	}
+
+	if(t = alpha(s, "plot")){
+		if(nrec == 0){
+			Bprint(&bout, "empty\n");
+			return;
+		}
+		plot(t);
+		return;
+	}
+
+	if(t = alpha(s, "astro")){
+		astro(t, 0);
+		return;
+	}
+
+	if(t = alpha(s, "plate")){
+		pplate(t);
+		return;
+	}
+
+	if(t = alpha(s, "gamma")){
+		while(*t==' ')
+			t++;
+		u = t;
+		x = strtod(t, &u);
+		if(u > t)
+			gam.gamma = x;
+		Bprint(&bout, "%.2f\n", gam.gamma);
+		return;
+	}
+
+	if(t = alpha(s, "keep")){
+		if(!cull(t, 1, 0))
+			return;
+		goto Print;
+	}
+
+	if(t = alpha(s, "drop")){
+		if(!cull(t, 0, 0))
+			return;
+		goto Print;
+	}
+
+	for(i=0; planet[i].name[0]; i++){
+		if(t = alpha(s, planet[i].name)){
+			if(doreset)
+				reset();
+			loadplanet(i, nil);
+			goto Print;
+		}
+	}
+
+	for(i=0; names[i].name; i++){
+		if(t = alpha(s, names[i].name)){
+			if(*t){
+				fprint(2, "syntax error in type\n");
+				return;
+			}
+			if(doreset)
+				reset();
+			loadtype(names[i].type);
+			goto Print;
+		}
+	}
+
+	switch(s[0]){
+	case '"':
+		starts = ++s;
+		while(*s != '"')
+			if(*s++ == 0){
+				fprint(2, "bad star name\n");
+				return;
+			}
+		*s = 0;
+		if(doreset)
+			reset();
+		j = nrec;
+		saoopen();
+		starts = fromgreek(starts);
+		for(i=0; i<NName; i++)
+			if(equal(starts, name[i].name)){
+				grow();
+				if(name[i].sao){
+					rec[j].type = NamedSAO;
+					rec[j].index = name[i].sao;
+				}
+				if(name[i].ngc){
+					rec[j].type = NamedNGC;
+					rec[j].index = name[i].ngc;
+				}
+				if(name[i].abell){
+					rec[j].type = NamedAbell;
+					rec[j].index = name[i].abell;
+				}
+				strcpy(rec[j].named.name, name[i].name);
+				j++;
+			}
+		if(parsename(starts))
+			for(i=0; i<NBayer; i++)
+				if(bayer[i].name[0]==parsed[0] &&
+				  (bayer[i].name[1]==parsed[1] || parsed[1]==0) &&
+				   bayer[i].name[2]==parsed[2]){
+					grow();
+					rec[j].type = NamedSAO;
+					rec[j].index = bayer[i].sao;
+					strncpy(rec[j].named.name, starts, sizeof(rec[j].named.name));
+					j++;
+				}
+		if(j == 0){
+			*s = '"';
+			goto NotFound;
+		}
+		break;
+
+	case '0': case '1': case '2': case '3': case '4':
+	case '5': case '6': case '7': case '8': case '9':
+		strtoul(s, &t, 10);
+		if(*t != 'h'){
+	BadCoords:
+			fprint(2, "bad coordinates %s\n", inputline);
+			break;
+		}
+		ra = DEG(getra(s));
+		while(*s && *s!=' ' && *s!='\t')
+			s++;
+		rah = ra/15;
+		ra = ra-rah*15;
+		ram = ra*4;
+		deg = strtol(s, &t, 10);
+		if(t == s)
+			goto BadCoords;
+		/* degree sign etc. is optional */
+		chartorune(&c, t);
+		if(c == 0xb0)
+			deg = DEG(getra(s));
+		if(doreset)
+			reset();
+		if(abs(deg)>=90 || rah>=24)
+			goto BadCoords;
+		if(!loadpatch(patch(rah, ram, deg)))
+			goto NotFound;
+		break;
+
+	default:
+		fprint(2, "unknown command %s\n", inputline);
+		return;
+	}
+
+    Print:
+	if(nrec == 0)
+		Bprint(&bout, "empty\n");
+	else if(nrec <= 2)
+		for(i=0; i<nrec; i++)
+			prrec(rec+i);
+	else
+		Bprint(&bout, "%ld items\n", nrec);
+	return;
+
+    NotFound:
+	fprint(2, "%s not found\n", inputline);
+	return;
+}
+
+char *ngctypes[] =
+{
+[Galaxy] 		"Gx",
+[PlanetaryN]	"Pl",
+[OpenCl]		"OC",
+[GlobularCl]	"Gb",
+[DiffuseN]		"Nb",
+[NebularCl]	"C+N",
+[Asterism]		"Ast",
+[Knot]		"Kt",
+[Triple]		"***",
+[Double]		"D*",
+[Single]		"*",
+[Uncertain]	"?",
+[Nonexistent]	"-",
+[Unknown]	" ",
+[PlateDefect]	"PD",
+};
+
+char*
+ngcstring(int d)
+{
+	if(d<Galaxy || d>PlateDefect)
+		return "can't happen";
+	return ngctypes[d];
+}
+
+short	descindex[NINDEX];
+
+void
+printnames(Record *r)
+{
+	int i, ok, done;
+
+	done = 0;
+	for(i=0; i<NName; i++){	/* stupid linear search! */
+		ok = 0;
+		if(r->type==SAO && r->index==name[i].sao)
+			ok = 1;
+		if(r->type==NGC && r->ngc.ngc==name[i].ngc)
+			ok = 1;
+		if(r->type==Abell && r->abell.abell==name[i].abell)
+			ok = 1;
+		if(ok){
+			if(done++ == 0)
+				Bprint(&bout, "\t");
+			Bprint(&bout, " \"%s\"", togreek(name[i].name));
+		}
+	}
+	if(done)
+		Bprint(&bout, "\n");
+}
+
+int
+equal(char *s1, char *s2)
+{
+	int c;
+
+	while(*s1){
+		if(*s1==' '){
+			while(*s1==' ')
+				s1++;
+			continue;
+		}
+		while(*s2==' ')
+			s2++;
+		c=*s2;
+		if('A'<=*s2 && *s2<='Z')
+			c^=' ';
+		if(*s1!=c)
+			return 0;
+		s1++, s2++;
+	}
+	return 1;
+}
+
+int
+parsename(char *s)
+{
+	char *blank;
+	int i;
+
+	blank = strchr(s, ' ');
+	if(blank==0 || strchr(blank+1, ' ') || strlen(blank+1)!=3)
+		return 0;
+	blank++;
+	parsed[0] = parsed[1] = parsed[2] = 0;
+	if('0'<=s[0] && s[0]<='9'){
+		i = atoi(s);
+		parsed[0] = i;
+		if(i > 100)
+			return 0;
+	}else{
+		for(i=1; i<=24; i++)
+			if(strncmp(greek[i], s, strlen(greek[i]))==0){
+				parsed[0]=100+i;
+				goto out;
+			}
+		return 0;
+	    out:
+		if('0'<=s[strlen(greek[i])] && s[strlen(greek[i])]<='9')
+			parsed[1]=s[strlen(greek[i])]-'0';
+	}
+	for(i=1; i<=88; i++)
+		if(strcmp(constel[i], blank)==0){
+			parsed[2] = i;
+			return 1;
+		}
+	return 0;
+}
+
+char*
+dist_grp(int dg)
+{
+	switch(dg){
+	default:
+		return "unknown";
+	case 1:
+		return "13.3-14.0";
+	case 2:
+		return "14.1-14.8";
+	case 3:
+		return "14.9-15.6";
+	case 4:
+		return "15.7-16.4";
+	case 5:
+		return "16.5-17.2";
+	case 6:
+		return "17.3-18.0";
+	case 7:
+		return ">18.0";
+	}
+}
+
+char*
+rich_grp(int dg)
+{
+	switch(dg){
+	default:
+		return "unknown";
+	case 0:
+		return "30-40";
+	case 1:
+		return "50-79";
+	case 2:
+		return "80-129";
+	case 3:
+		return "130-199";
+	case 4:
+		return "200-299";
+	case 5:
+		return ">=300";
+	}
+}
+
+char*
+nameof(Record *r)
+{
+	NGCrec *n;
+	SAOrec *s;
+	Abellrec *a;
+	static char buf[128];
+	int i;
+
+	switch(r->type){
+	default:
+		return nil;
+	case SAO:
+		s = &r->sao;
+		if(s->name[0] == 0)
+			return nil;
+		if(s->name[0] >= 100){
+			i = snprint(buf, sizeof buf, "%C", greeklet[s->name[0]-100]);
+			if(s->name[1])
+				i += snprint(buf+i, sizeof buf-i, "%d", s->name[1]);
+		}else
+			i = snprint(buf, sizeof buf, " %d", s->name[0]);
+		snprint(buf+i, sizeof buf-i, " %s", constel[(uchar)s->name[2]]);
+		break;
+	case NGC:
+		n = &r->ngc;
+		if(n->type >= Uncertain)
+			return nil;
+		if(n->ngc <= NNGC)
+			snprint(buf, sizeof buf, "NGC%4d ", n->ngc);
+		else
+			snprint(buf, sizeof buf, "IC%4d ", n->ngc-NNGC);
+		break;
+	case Abell:
+		a = &r->abell;
+		snprint(buf, sizeof buf, "Abell%4d", a->abell);
+		break;
+	}
+	return buf;
+}
+
+void
+prrec(Record *r)
+{
+	NGCrec *n;
+	SAOrec *s;
+	Abellrec *a;
+	Planetrec *p;
+	int i, rah, ram, dec, nn;
+	long key;
+
+	if(r) switch(r->type){
+	default:
+		fprint(2, "can't prrec type %d\n", r->type);
+		exits("type");
+
+	case Planet:
+		p = &r->planet;
+		Bprint(&bout, "%s", p->name);
+		Bprint(&bout, "\t%s %s",
+			hms(angle(p->ra)),
+			dms(angle(p->dec)));
+		Bprint(&bout, " %3.2f° %3.2f°",
+			p->az/(double)MILLIARCSEC, p->alt/(double)MILLIARCSEC);
+		Bprint(&bout, " %s",
+			ms(angle(p->semidiam)));
+		if(r->index <= 1)
+			Bprint(&bout, " %g", p->phase);
+		Bprint(&bout, "\n");
+		break;
+
+	case NGC:
+		n = &r->ngc;
+		if(n->ngc <= NNGC)
+			Bprint(&bout, "NGC%4d ", n->ngc);
+		else
+			Bprint(&bout, "IC%4d ", n->ngc-NNGC);
+		Bprint(&bout, "%s ", ngcstring(n->type));
+		if(n->mag == UNKNOWNMAG)
+			Bprint(&bout, "----");
+		else
+			Bprint(&bout, "%.1f%c", n->mag/10.0, n->magtype);
+		Bprint(&bout, "\t%s %s\t%c%.1f'\n",
+			hm(angle(n->ra)),
+			dm(angle(n->dec)),
+			n->diamlim,
+			DEG(angle(n->diam))*60.);
+		prdesc(n->desc, desctab, descindex);
+		printnames(r);
+		break;
+
+	case Abell:
+		a = &r->abell;
+		Bprint(&bout, "Abell%4d  %.1f %.2f° %dMpc", a->abell, a->mag10/10.0,
+			DEG(angle(a->rad)), a->dist);
+		Bprint(&bout, "\t%s %s\t%.2f %.2f\n",
+			hm(angle(a->ra)),
+			dm(angle(a->dec)),
+			DEG(angle(a->glat)),
+			DEG(angle(a->glong)));
+		Bprint(&bout, "\tdist grp: %s  rich grp: %s  %d galaxies/°²\n",
+			dist_grp(a->distgrp),
+			rich_grp(a->richgrp),
+			a->pop);
+		printnames(r);
+		break;
+
+	case SAO:
+		s = &r->sao;
+		Bprint(&bout, "SAO%6ld  ", r->index);
+		if(s->mag==UNKNOWNMAG)
+			Bprint(&bout, "---");
+		else
+			Bprint(&bout, "%.1f", s->mag/10.0);
+		if(s->mpg==UNKNOWNMAG)
+			Bprint(&bout, ",---");
+		else
+			Bprint(&bout, ",%.1f", s->mpg/10.0);
+		Bprint(&bout, "  %s %s  %.4fs %.3f\"",
+			hms(angle(s->ra)),
+			dms(angle(s->dec)),
+			DEG(angle(s->dra))*(4*60),
+			DEG(angle(s->ddec))*(60*60));
+		Bprint(&bout, "  %.3s %c %.2s %ld %d",
+			s->spec, s->code, s->compid, s->hd, s->hdcode);
+		if(s->name[0])
+			Bprint(&bout, " \"%s\"", nameof(r));
+		Bprint(&bout, "\n");
+		printnames(r);
+		break;
+
+	case Patch:
+		radec(r->index, &rah, &ram, &dec);
+		Bprint(&bout, "%dh%dm %d°", rah, ram, dec);
+		key = r->patch.key[0];
+		Bprint(&bout, " %s", constel[key&0xFF]);
+		if((key>>=8) & 0xFF)
+			Bprint(&bout, " %s", constel[key&0xFF]);
+		if((key>>=8) & 0xFF)
+			Bprint(&bout, " %s", constel[key&0xFF]);
+		if((key>>=8) & 0xFF)
+			Bprint(&bout, " %s", constel[key&0xFF]);
+		for(i=1; i<r->patch.nkey; i++){
+			key = r->patch.key[i];
+			switch(key&0x3F){
+			case SAO:
+				Bprint(&bout, " SAO%ld", (key>>8)&0xFFFFFF);
+				break;
+			case Abell:
+				Bprint(&bout, " Abell%ld", (key>>8)&0xFFFFFF);
+				break;
+			default:	/* NGC */
+				nn = (key>>16)&0xFFFF;
+				if(nn > NNGC)
+					Bprint(&bout, " IC%d", nn-NNGC);
+				else
+					Bprint(&bout, " NGC%d", nn);
+				Bprint(&bout, "(%s)", ngcstring(key&0x3F));
+				break;
+			}
+		}
+		Bprint(&bout, "\n");
+		break;
+
+	case NGCN:
+		if(r->index <= NNGC)
+			Bprint(&bout, "NGC%ld\n", r->index);
+		else
+			Bprint(&bout, "IC%ld\n", r->index-NNGC);
+		break;
+
+	case NamedSAO:
+		Bprint(&bout, "SAO%ld \"%s\"\n", r->index, togreek(r->named.name));
+		break;
+
+	case NamedNGC:
+		if(r->index <= NNGC)
+			Bprint(&bout, "NGC%ld \"%s\"\n", r->index, togreek(r->named.name));
+		else
+			Bprint(&bout, "IC%ld \"%s\"\n", r->index-NNGC, togreek(r->named.name));
+		break;
+
+	case NamedAbell:
+		Bprint(&bout, "Abell%ld \"%s\"\n", r->index, togreek(r->named.name));
+		break;
+
+	case PatchC:
+		radec(r->index, &rah, &ram, &dec);
+		Bprint(&bout, "%dh%dm %d\n", rah, ram, dec);
+		break;
+	}
+}
diff --git a/src/cmd/scat/sky.h b/src/cmd/scat/sky.h
new file mode 100644
index 00000000..420a2a9a
--- /dev/null
+++ b/src/cmd/scat/sky.h
@@ -0,0 +1,413 @@
+#define DIR	"#9/sky"
+/*
+ *	This code reflects many years of changes.  There remain residues
+ *		of prior implementations.
+ *
+ *	Keys:
+ *		32 bits long. High 26 bits are encoded as described below.
+ *		Low 6 bits are types:
+ *
+ *		Patch is ~ one square degree of sky.  It points to an otherwise
+ *			anonymous list of Catalog keys.  The 0th key is special:
+ *			it contains up to 4 constellation identifiers.
+ *		Catalogs (SAO,NGC,M,...) are:
+ *			31.........8|76|543210
+ *			  catalog # |BB|catalog name
+ *			BB is two bits of brightness:
+ *				00	-inf <  m <=  7
+ *				01	   7 <  m <= 10
+ *				10	  10 <  m <= 13
+ *				11	  13 <  m <  inf
+ *			The BB field is a dreg, and correct only for SAO and NGC.
+ *			IC(n) is just NGC(n+7840)
+ *		Others should be self-explanatory.
+ *	
+ *	Records:
+ *
+ *	Star is an SAOrec
+ *	Galaxy, PlanetaryN, OpenCl, GlobularCl, DiffuseN, etc., are NGCrecs.
+ *	Abell is an Abellrec
+ *	The Namedxxx records hold a name and a catalog entry; they result from
+ *		name lookups.
+ */
+
+typedef enum
+{
+	Planet,
+	Patch,
+	SAO,
+	NGC,
+	M,
+	Constel_deprecated,
+	Nonstar_deprecated,
+	NamedSAO,
+	NamedNGC,
+	NamedAbell,
+	Abell,
+	/* NGC types */
+	Galaxy,
+	PlanetaryN,
+	OpenCl,
+	GlobularCl,
+	DiffuseN,
+	NebularCl,
+	Asterism,
+	Knot,
+	Triple,
+	Double,
+	Single,
+	Uncertain,
+	Nonexistent,
+	Unknown,
+	PlateDefect,
+	/* internal */
+	NGCN,
+	PatchC,
+	NONGC,
+}Type;
+
+enum
+{
+	/*
+	 * parameters for plate
+	 */
+	Pppo1	= 0,
+	Pppo2,
+	Pppo3,
+	Pppo4,
+	Pppo5,
+	Pppo6,
+	Pamdx1,
+	Pamdx2,
+	Pamdx3,
+	Pamdx4,
+	Pamdx5,
+	Pamdx6,
+	Pamdx7,
+	Pamdx8,
+	Pamdx9,
+	Pamdx10,
+	Pamdx11,
+	Pamdx12,
+	Pamdx13,
+	Pamdx14,
+	Pamdx15,
+	Pamdx16,
+	Pamdx17,
+	Pamdx18,
+	Pamdx19,
+	Pamdx20,
+	Pamdy1,
+	Pamdy2,
+	Pamdy3,
+	Pamdy4,
+	Pamdy5,
+	Pamdy6,
+	Pamdy7,
+	Pamdy8,
+	Pamdy9,
+	Pamdy10,
+	Pamdy11,
+	Pamdy12,
+	Pamdy13,
+	Pamdy14,
+	Pamdy15,
+	Pamdy16,
+	Pamdy17,
+	Pamdy18,
+	Pamdy19,
+	Pamdy20,
+	Ppltscale,
+	Pxpixelsz,
+	Pypixelsz,
+	Ppltra,
+	Ppltrah,
+	Ppltram,
+	Ppltras,
+	Ppltdec,
+	Ppltdecd,
+	Ppltdecm,
+	Ppltdecs,
+	Pnparam,
+};
+
+#define	UNKNOWNMAG	32767
+#define	NPlanet			20
+
+typedef float	Angle;	/* in radians */
+typedef long	DAngle;	/* on disk: in units of milliarcsec */
+typedef short	Mag;	/* multiplied by 10 */
+typedef long	Key;	/* known to be 4 bytes, unfortunately */
+
+/*
+ * All integers are stored in little-endian order.
+ */
+typedef struct NGCrec NGCrec;
+struct NGCrec{
+	DAngle	ra;
+	DAngle	dec;
+	DAngle	dummy1;	/* compatibility with old RNGC version */
+	DAngle	diam;
+	Mag	mag;
+	short	ngc;	/* if >NNGC, IC number is ngc-NNGC */
+	char	diamlim;
+	char	type;
+	char	magtype;
+	char	dummy2;
+	char	desc[52];	/* 0-terminated Dreyer description */
+};
+
+typedef struct Abellrec Abellrec;
+struct Abellrec{
+	DAngle	ra;
+	DAngle	dec;
+	DAngle	glat;
+	DAngle	glong;
+	Mag	mag10;	/* mag of 10th brightest cluster member; in same place as ngc.mag*/
+	short	abell;
+	DAngle	rad;
+	short	pop;
+	short	dist;
+	char	distgrp;
+	char	richgrp;
+	char	flag;
+	char	pad;
+};
+
+typedef struct Planetrec Planetrec;
+struct Planetrec{
+	DAngle	ra;
+	DAngle	dec;
+	DAngle	az;
+	DAngle	alt;
+	DAngle	semidiam;
+	double	phase;
+	char		name[16];
+};
+
+/*
+ * Star names: 0,0==unused.  Numbers are name[0]=1,..,99.
+ * Greek letters are alpha=101, etc.
+ * Constellations are alphabetical order by abbreviation, and=1, etc.
+ */
+typedef struct SAOrec SAOrec;
+struct SAOrec{
+	DAngle	ra;
+	DAngle	dec;
+	DAngle	dra;
+	DAngle	ddec;
+	Mag	mag;		/* visual */
+	Mag	mpg;
+	char	spec[3];
+	char	code;
+	char	compid[2];
+	char	hdcode;
+	char	pad1;
+	long	hd;		/* HD catalog number */
+	char	name[3];	/* name[0]=alpha name[1]=2 name[3]=ori */
+	char	nname;		/* number of prose names */
+	/* 36 bytes to here */
+};
+
+typedef struct Mindexrec Mindexrec;
+struct Mindexrec{	/* code knows the bit patterns in here; this is a long */
+	char	m;		/* M number */
+	char	dummy;
+	short	ngc;
+};
+
+typedef struct Bayerec Bayerec;
+struct Bayerec{
+	long	sao;
+	char	name[3];
+	char	pad;
+};
+
+/*
+ * Internal form
+ */
+
+typedef struct Namedrec Namedrec;
+struct Namedrec{
+	char	name[36];
+};
+
+typedef struct Namerec Namerec;
+struct Namerec{
+	long	sao;
+	long	ngc;
+	long	abell;
+	char	name[36];	/* null terminated */
+};
+
+typedef struct Patchrec Patchrec;
+struct Patchrec{
+	int	nkey;
+	long	key[60];
+};
+
+typedef struct Record Record;
+struct Record{
+	Type	type;
+	long	index;
+	union{
+		SAOrec	sao;
+		NGCrec	ngc;
+		Abellrec	abell;
+		Namedrec	named;
+		Patchrec	patch;
+		Planetrec	planet;
+		/* PatchCrec is empty */
+	};
+};
+
+typedef struct Name Name;
+struct Name{
+	char	*name;
+	int	type;
+};
+
+typedef	struct	Plate	Plate;
+struct	Plate
+{
+	char	rgn[7];
+	char	disk;
+	Angle	ra;
+	Angle	dec;
+};
+
+typedef	struct	Header	Header;
+struct	Header
+{
+	float	param[Pnparam];
+	int	amdflag;
+
+	float	x;
+	float	y;
+	float	xi;
+	float	eta;
+};
+typedef	long	Pix;
+
+typedef struct	Img Img;
+struct	Img
+{
+	int	nx;
+	int	ny;	/* ny is the fast-varying dimension */
+	Pix	a[1];
+};
+
+#define	RAD(x)	((x)*PI_180)
+#define	DEG(x)	((x)/PI_180)
+#define	ARCSECONDS_PER_RADIAN	(DEG(1)*3600)
+#define	MILLIARCSEC	(1000*60*60)
+
+int	nplate;
+Plate	plate[2000];		/* needs to go to 2000 when the north comes */
+double	PI_180;
+double	TWOPI;
+double	LN2;
+int	debug;
+struct
+{
+	float	min;
+	float	max;
+	float	gamma;
+	float	absgamma;
+	float	mult1;
+	float	mult2;
+	int	neg;
+} gam;
+
+typedef struct Picture Picture;
+struct Picture
+{
+	int	minx;
+	int	miny;
+	int	maxx;
+	int	maxy;
+	char	name[16];
+	uchar	*data;
+};
+
+#ifndef _DRAW_H_
+typedef struct Image Image;
+#endif
+
+extern	double	PI_180;
+extern	double	TWOPI;
+extern	char	*progname;
+extern	char	*desctab[][2];
+extern	Name	names[];
+extern	Record	*rec;
+extern	long		nrec;
+extern	Planetrec	*planet;
+/* for bbox: */
+extern	int		folded;
+extern	DAngle	ramin;
+extern	DAngle	ramax;
+extern	DAngle	decmin;
+extern	DAngle	decmax;
+extern	Biobuf	bout;
+
+extern void saoopen(void);
+extern void ngcopen(void);
+extern void patchopen(void);
+extern void mopen(void);
+extern void constelopen(void);
+extern void lowercase(char*);
+extern void lookup(char*, int);
+extern int typetab(int);
+extern char*ngcstring(int);
+extern char*skip(int, char*);
+extern void prrec(Record*);
+extern int equal(char*, char*);
+extern int parsename(char*);
+extern void radec(int, int*, int*, int*);
+extern int btag(short);
+extern long patcha(Angle, Angle);
+extern long patch(int, int, int);
+extern char*hms(Angle);
+extern char*dms(Angle);
+extern char*ms(Angle);
+extern char*hm(Angle);
+extern char*dm(Angle);
+extern char*deg(Angle);
+extern char*hm5(Angle);
+extern long dangle(Angle);
+extern Angle angle(DAngle);
+extern void prdesc(char*, char*(*)[2], short*);
+extern double	xsqrt(double);
+extern Angle	dist(Angle, Angle, Angle, Angle);
+extern Header*	getheader(char*);
+extern char*	getword(char*, char*);
+extern void	amdinv(Header*, Angle, Angle, float, float);
+extern void	ppoinv(Header*, Angle, Angle);
+extern void	xypos(Header*, Angle, Angle, float, float);
+extern void	traneqstd(Header*, Angle, Angle);
+extern Angle	getra(char*);
+extern Angle	getdec(char*);
+extern void	getplates(void);
+extern Img*	dssread(char*);
+extern void	hinv(Pix*, int, int);
+extern int	input_bit(Biobuf*);
+extern int	input_nbits(Biobuf*, int);
+extern int	input_huffman(Biobuf*);
+extern	int	input_nybble(Biobuf*);
+extern void	qtree_decode(Biobuf*, Pix*, int, int, int, int);
+extern void	start_inputing_bits(void);
+extern Picture*	image(Angle, Angle, Angle, Angle);
+extern char*	dssmount(int);
+extern int	dogamma(Pix);
+extern void	displaypic(Picture*);
+extern void	displayimage(Image*);
+extern void	plot(char*);
+extern void	astro(char*, int);
+extern char*	alpha(char*, char*);
+extern char*	skipbl(char*);
+extern void	flatten(void);
+extern int		bbox(long, long, int);
+extern int		inbbox(DAngle, DAngle);
+extern char*	nameof(Record*);
+
+#define	NINDEX	400
diff --git a/src/cmd/scat/strings.c b/src/cmd/scat/strings.c
new file mode 100644
index 00000000..e60246c7
--- /dev/null
+++ b/src/cmd/scat/strings.c
@@ -0,0 +1,52 @@
+char *greek[]={ 0,	/* 1-indexed */
+	"alpha", "beta", "gamma", "delta", "epsilon", "zeta", "eta", "theta",
+	"iota", "kappa", "lambda", "mu", "nu", "xsi", "omicron", "pi", "rho",
+	"sigma", "tau", "upsilon", "phi", "chi", "psi", "omega",
+};
+
+Rune greeklet[]={ 0,
+	0x3b1, 0x3b2, 0x3b3, 0x3b4, 0x3b5, 0x3b6, 0x3b7, 0x3b8, 0x3b9, 0x3ba, 0x3bb,
+	0x3bc, 0x3bd, 0x3be, 0x3bf, 0x3c0, 0x3c1, 0x3c3, 0x3c4, 0x3c5, 0x3c6, 0x3c7,
+	0x3c8, 0x3c9,
+};
+
+char *constel[]={ 0,	/* 1-indexed */
+	"and", "ant", "aps", "aql", "aqr", "ara", "ari", "aur", "boo", "cae",
+	"cam", "cap", "car", "cas", "cen", "cep", "cet", "cha", "cir", "cma",
+	"cmi", "cnc", "col", "com", "cra", "crb", "crt", "cru", "crv", "cvn",
+	"cyg", "del", "dor", "dra", "equ", "eri", "for", "gem", "gru", "her",
+	"hor", "hya", "hyi", "ind", "lac", "leo", "lep", "lib", "lmi", "lup",
+	"lyn", "lyr", "men", "mic", "mon", "mus", "nor", "oct", "oph", "ori",
+	"pav", "peg", "per", "phe", "pic", "psa", "psc", "pup", "pyx", "ret",
+	"scl", "sco", "sct", "ser", "sex", "sge", "sgr", "tau", "tel", "tra",
+	"tri", "tuc", "uma", "umi", "vel", "vir", "vol", "vul",
+};
+Name names[]={
+	"gx",	Galaxy,
+	"pl",	PlanetaryN,
+	"oc",	OpenCl,
+	"gb",	GlobularCl,
+	"nb",	DiffuseN,
+	"c+n",NebularCl,
+	"ast",	Asterism,
+	"kt",	Knot,
+	"***",	Triple,
+	"d*",	Double,
+	"*",	Single,
+	"pd",	PlateDefect,
+	"galaxy",	Galaxy,
+	"planetary",	PlanetaryN,
+	"opencluster",	OpenCl,
+	"globularcluster",	GlobularCl,
+	"nebula",	DiffuseN,
+	"nebularcluster",NebularCl,
+	"asterism",	Asterism,
+	"knot",	Knot,
+	"triple",	Triple,
+	"double",	Double,
+	"single",	Single,
+	"nonexistent",	Nonexistent,
+	"unknown",	Unknown,
+	"platedefect",	PlateDefect,
+	0,		0,
+};
diff --git a/src/cmd/scat/util.c b/src/cmd/scat/util.c
new file mode 100644
index 00000000..33dca378
--- /dev/null
+++ b/src/cmd/scat/util.c
@@ -0,0 +1,368 @@
+#include <u.h>
+#include <libc.h>
+#include <bio.h>
+#include "sky.h"
+
+double	PI_180	= 0.0174532925199432957692369;
+double	TWOPI	= 6.2831853071795864769252867665590057683943387987502;
+double	LN2	= 0.69314718055994530941723212145817656807550013436025;
+static double angledangle=(180./PI)*MILLIARCSEC;
+
+#define rint scatrint
+
+int
+rint(char *p, int n)
+{
+	int i=0;
+
+	while(*p==' ' && n)
+		p++, --n;
+	while(n--)
+		i=i*10+*p++-'0';
+	return i;
+}
+
+DAngle
+dangle(Angle angle)
+{
+	return angle*angledangle;
+}
+
+Angle
+angle(DAngle dangle)
+{
+	return dangle/angledangle;
+}
+
+double
+rfloat(char *p, int n)
+{
+	double i, d=0;
+
+	while(*p==' ' && n)
+		p++, --n;
+	if(*p == '+')
+		return rfloat(p+1, n-1);
+	if(*p == '-')
+		return -rfloat(p+1, n-1);
+	while(*p == ' ' && n)
+		p++, --n;
+	if(n == 0)
+		return 0.0;
+	while(n-- && *p!='.')
+		d = d*10+*p++-'0';
+	if(n <= 0)
+		return d;
+	p++;
+	i = 1;
+	while(n--)
+		d+=(*p++-'0')/(i*=10.);
+	return d;
+}
+
+int
+sign(int c)
+{
+	if(c=='-')
+		return -1;
+	return 1;
+}
+
+char*
+hms(Angle a)
+{
+	static char buf[20];
+	double x;
+	int h, m, s, ts;
+
+	x=DEG(a)/15;
+	x += 0.5/36000.;	/* round up half of 0.1 sec */
+	h = floor(x);
+	x -= h;
+	x *= 60;
+	m = floor(x);
+	x -= m;
+	x *= 60;
+	s = floor(x);
+	x -= s;
+	ts = 10*x;
+	sprint(buf, "%dh%.2dm%.2d.%ds", h, m, s, ts);
+	return buf;
+}
+
+char*
+dms(Angle a)
+{
+	static char buf[20];
+	double x;
+	int sign, d, m, s, ts;
+
+	x = DEG(a);
+	sign='+';
+	if(a<0){
+		sign='-';
+		x=-x;
+	}
+	x += 0.5/36000.;	/* round up half of 0.1 arcsecond */
+	d = floor(x);
+	x -= d;
+	x *= 60;
+	m = floor(x);
+	x -= m;
+	x *= 60;
+	s = floor(x);
+	x -= s;
+	ts = floor(10*x);
+	sprint(buf, "%c%d°%.2d'%.2d.%d\"", sign, d, m, s, ts);
+	return buf;
+}
+
+char*
+ms(Angle a)
+{
+	static char buf[20];
+	double x;
+	int d, m, s, ts;
+
+	x = DEG(a);
+	x += 0.5/36000.;	/* round up half of 0.1 arcsecond */
+	d = floor(x);
+	x -= d;
+	x *= 60;
+	m = floor(x);
+	x -= m;
+	x *= 60;
+	s = floor(x);
+	x -= s;
+	ts = floor(10*x);
+	if(d != 0)
+		sprint(buf, "%d°%.2d'%.2d.%d\"", d, m, s, ts);
+	else
+		sprint(buf, "%.2d'%.2d.%d\"", m, s, ts);
+	return buf;
+}
+
+char*
+hm(Angle a)
+{
+	static char buf[20];
+	double x;
+	int h, m, n;
+
+	x = DEG(a)/15;
+	x += 0.5/600.;	/* round up half of tenth of minute */
+	h = floor(x);
+	x -= h;
+	x *= 60;
+	m = floor(x);
+	x -= m;
+	x *= 10;
+	n = floor(x);
+	sprint(buf, "%dh%.2d.%1dm", h, m, n);
+	return buf;
+}
+
+char*
+hm5(Angle a)
+{
+	static char buf[20];
+	double x;
+	int h, m;
+
+	x = DEG(a)/15;
+	x += 2.5/60.;	/* round up 2.5m */
+	h = floor(x);
+	x -= h;
+	x *= 60;
+	m = floor(x);
+	m -= m % 5;
+	sprint(buf, "%dh%.2dm", h, m);
+	return buf;
+}
+
+char*
+dm(Angle a)
+{
+	static char buf[20];
+	double x;
+	int sign, d, m, n;
+
+	x = DEG(a);
+	sign='+';
+	if(a<0){
+		sign='-';
+		x=-x;
+	}
+	x += 0.5/600.;	/* round up half of tenth of arcminute */
+	d = floor(x);
+	x -= d;
+	x *= 60;
+	m = floor(x);
+	x -= m;
+	x *= 10;
+	n = floor(x);
+	sprint(buf, "%c%d°%.2d.%.1d'", sign, d, m, n);
+	return buf;
+}
+
+char*
+deg(Angle a)
+{
+	static char buf[20];
+	double x;
+	int sign, d;
+
+	x = DEG(a);
+	sign='+';
+	if(a<0){
+		sign='-';
+		x=-x;
+	}
+	x += 0.5;	/* round up half degree */
+	d = floor(x);
+	sprint(buf, "%c%d°", sign, d);
+	return buf;
+}
+
+char*
+getword(char *ou, char *in)
+{
+	int c;
+
+	for(;;) {
+		c = *in++;
+		if(c == ' ' || c == '\t')
+			continue;
+		if(c == 0)
+			return 0;
+		break;
+	}
+
+	if(c == '\'')
+		for(;;) {
+			if(c >= 'A' && c <= 'Z')
+				c += 'a' - 'A';
+			*ou++ = c;
+			c = *in++;
+			if(c == 0)
+				return 0;
+			if(c == '\'') {
+				*ou = 0;
+				return in-1;
+			}
+		}
+	for(;;) {
+		if(c >= 'A' && c <= 'Z')
+			c += 'a' - 'A';
+		*ou++ = c;
+		c = *in++;
+		if(c == ' ' || c == '\t' || c == 0) {
+			*ou = 0;
+			return in-1;
+		}
+	}
+}
+
+/*
+ * Read formatted angle.  Must contain no embedded blanks
+ */
+Angle
+getra(char *p)
+{
+	Rune r;
+	char *q;
+	Angle f, d;
+	int neg;
+
+	neg = 0;
+	d = 0;
+	while(*p == ' ')
+		p++;
+	for(;;) {
+		if(*p == ' ' || *p=='\0')
+			goto Return;
+		if(*p == '-') {
+			neg = 1;
+			p++;
+		}
+		if(*p == '+') {
+			neg = 0;
+			p++;
+		}
+		q = p;
+		f = strtod(p, &q);
+		if(q > p) {
+			p = q;
+		}
+		p += chartorune(&r, p);
+		switch(r) {
+		default:
+		Return:
+			if(neg)
+				d = -d;
+			return RAD(d);
+		case 'h':
+			d += f*15;
+			break;
+		case 'm':
+			d += f/4;
+			break;
+		case 's':
+			d += f/240;
+			break;
+		case 0xB0:	/* ° */
+			d += f;
+			break;
+		case '\'':
+			d += f/60;
+			break;
+		case '\"':
+			d += f/3600;
+			break;
+		}
+	}
+	return 0;
+}
+
+double
+xsqrt(double a)
+{
+
+	if(a < 0)
+		return 0;
+	return sqrt(a);
+}
+
+Angle
+dist(Angle ra1, Angle dec1, Angle ra2, Angle dec2)
+{
+	double a;
+
+	a = sin(dec1) * sin(dec2) +
+		cos(dec1) * cos(dec2) *
+		cos(ra1 - ra2);
+	a = atan2(xsqrt(1 - a*a), a);
+	if(a < 0)
+		a = -a;
+	return a;
+}
+
+int
+dogamma(Pix c)
+{
+	float f;
+
+	f = c - gam.min;
+	if(f < 1)
+		f = 1;
+
+	if(gam.absgamma == 1)
+		c = f * gam.mult2;
+	else
+		c = exp(log(f*gam.mult1) * gam.absgamma) * 255;
+	if(c > 255)
+		c = 255;
+	if(gam.neg)
+		c = 255-c;
+	return c;
+}