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

1 files changed, 952 insertions, 0 deletions

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]);
+}