harvey/sys/src/cmd/scat/plot.c

#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;
	display = initdisplay(nil, nil, drawerror);
	if(display == nil){
		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))(Loc*, 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]);
		execl("/bin/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, "/bin/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]);
}