1 #include <u.h>
2 #include <libc.h>
3 #include <bio.h>
4 #include <draw.h>
5 #include <event.h>
6 #include <ctype.h>
7 #include "map.h"
8 #undef	RAD
9 #undef	TWOPI
10 #include "sky.h"
11 
12 	/* convert to milliarcsec */
13 static int32	c = MILLIARCSEC;	/* 1 degree */
14 static int32	m5 = 1250*60*60;	/* 5 minutes of ra */
15 
16 DAngle	ramin;
17 DAngle	ramax;
18 DAngle	decmin;
19 DAngle	decmax;
20 int		folded;
21 
22 Image	*grey;
23 Image	*alphagrey;
24 Image	*green;
25 Image	*lightblue;
26 Image	*lightgrey;
27 Image	*ocstipple;
28 Image	*suncolor;
29 Image	*mooncolor;
30 Image	*shadowcolor;
31 Image	*mercurycolor;
32 Image	*venuscolor;
33 Image	*marscolor;
34 Image	*jupitercolor;
35 Image	*saturncolor;
36 Image	*uranuscolor;
37 Image	*neptunecolor;
38 Image	*plutocolor;
39 Image	*cometcolor;
40 
41 Planetrec	*planet;
42 
43 int32	mapx0, mapy0;
44 int32	mapra, mapdec;
45 double	mylat, mylon, mysid;
46 double	mapscale;
47 double	maps;
48 int (*projection)(struct place*, double*, double*);
49 
50 char *fontname = "/lib/font/bit/luc/unicode.6.font";
51 
52 /* types Coord and Loc correspond to types in map(3) thus:
53    Coord == struct coord;
54    Loc == struct place, except longitudes are measured
55      positive east for Loc and west for struct place
56 */
57 
58 typedef struct Xyz Xyz;
59 typedef struct coord Coord;
60 typedef struct Loc Loc;
61 
62 struct Xyz{
63 	double x,y,z;
64 };
65 
66 struct Loc{
67 	Coord nlat, elon;
68 };
69 
70 Xyz north = { 0, 0, 1 };
71 
72 int
73 plotopen(void)
74 {
75 	if(display != nil)
76 		return 1;
77 	if(initdraw(drawerror, nil, nil) < 0){
78 		fprint(2, "initdisplay failed: %r\n");
79 		return -1;
80 	}
81 	grey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x777777FF);
82 	lightgrey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0xAAAAAAFF);
83 	alphagrey = allocimage(display, Rect(0, 0, 2, 2), RGBA32, 1, 0x77777777);
84 	green = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x00AA00FF);
85 	lightblue = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x009EEEFF);
86 	ocstipple = allocimage(display, Rect(0, 0, 2, 2), CMAP8, 1, 0xAAAAAAFF);
87 	draw(ocstipple, Rect(0, 0, 1, 1), display->black, nil, ZP);
88 	draw(ocstipple, Rect(1, 1, 2, 2), display->black, nil, ZP);
89 
90 	suncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFF77FF);
91 	mooncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAAAFF);
92 	shadowcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x00000055);
93 	mercurycolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFAAAAFF);
94 	venuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);
95 	marscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFF5555FF);
96 	jupitercolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFFAAFF);
97 	saturncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAFFAAFF);
98 	uranuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77DDDDFF);
99 	neptunecolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77FF77FF);
100 	plutocolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x7777FFFF);
101 	cometcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);
102 	font = openfont(display, fontname);
103 	if(font == nil)
104 		fprint(2, "warning: no font %s: %r\n", fontname);
105 	return 1;
106 }
107 
108 /*
109  * Function heavens() for setting up observer-eye-view
110  * sky charts, plus subsidiary functions.
111  * Written by Doug McIlroy.
112  */
113 
114 /* heavens(nlato, wlono, nlatc, wlonc) sets up a map(3)-style
115    coordinate change (by calling orient()) and returns a
116    projection function heavens for calculating a star map
117    centered on nlatc,wlonc and rotated so it will appear
118    upright as seen by a ground observer at nlato,wlono.
119    all coordinates (north latitude and west longitude)
120    are in degrees.
121 */
122 
123 Coord
124 mkcoord(double degrees)
125 {
126 	Coord c;
127 
128 	c.l = degrees*PI/180;
129 	c.s = sin(c.l);
130 	c.c = cos(c.l);
131 	return c;
132 }
133 
134 Xyz
135 ptov(struct place p)
136 {
137 	Xyz v;
138 
139 	v.z = p.nlat.s;
140 	v.x = p.nlat.c * p.wlon.c;
141 	v.y = -p.nlat.c * p.wlon.s;
142 	return v;
143 }
144 
145 double
146 dot(Xyz a, Xyz b)
147 {
148 	return a.x*b.x + a.y*b.y + a.z*b.z;
149 }
150 
151 Xyz
152 cross(Xyz a, Xyz b)
153 {
154 	Xyz v;
155 
156 	v.x = a.y*b.z - a.z*b.y;
157 	v.y = a.z*b.x - a.x*b.z;
158 	v.z = a.x*b.y - a.y*b.x;
159 	return v;
160 }
161 
162 double
163 len(Xyz a)
164 {
165 	return sqrt(dot(a, a));
166 }
167 
168 /* An azimuth vector from a to b is a tangent to
169    the sphere at a pointing along the (shorter)
170    great-circle direction to b.  It lies in the
171    plane of the vectors a and b, is perpendicular
172    to a, and has a positive compoent along b,
173    provided a and b span a 2-space.  Otherwise
174    it is null.  (aXb)Xa, where X denotes cross
175    product, is such a vector. */
176 
177 Xyz
178 azvec(Xyz a, Xyz b)
179 {
180 	return cross(cross(a,b), a);
181 }
182 
183 /* Find the azimuth of point b as seen
184    from point a, given that a is distinct
185    from either pole and from b */
186 
187 double
188 azimuth(Xyz a, Xyz b)
189 {
190 	Xyz aton = azvec(a,north);
191 	Xyz atob = azvec(a,b);
192 	double lenaton = len(aton);
193 	double lenatob = len(atob);
194 	double az = acos(dot(aton,atob)/(lenaton*lenatob));
195 
196 	if(dot(b,cross(north,a)) < 0)
197 		az = -az;
198 	return az;
199 }
200 
201 /* Find the rotation (3rd argument of orient() in the
202    map projection package) for the map described
203    below.  The return is radians; it must be converted
204    to degrees for orient().
205 */
206 
207 double
208 rot(struct place center, struct place zenith)
209 {
210 	Xyz cen = ptov(center);
211 	Xyz zen = ptov(zenith);
212 
213 	if(cen.z > 1-FUZZ) 	/* center at N pole */
214 		return PI + zenith.wlon.l;
215 	if(cen.z < FUZZ-1)		/* at S pole */
216 		return -zenith.wlon.l;
217 	if(fabs(dot(cen,zen)) > 1-FUZZ)	/* at zenith */
218 		return 0;
219 	return azimuth(cen, zen);
220 }
221 
222 int (*
223 heavens(double zlatdeg, double zlondeg, double clatdeg, double clondeg))(struct place*, double*, double*)
224 {
225 	struct place center;
226 	struct place zenith;
227 
228 	center.nlat = mkcoord(clatdeg);
229 	center.wlon = mkcoord(clondeg);
230 	zenith.nlat = mkcoord(zlatdeg);
231 	zenith.wlon = mkcoord(zlondeg);
232 	projection = stereographic();
233 	orient(clatdeg, clondeg, rot(center, zenith)*180/PI);
234 	return projection;
235 }
236 
237 int
238 maptoxy(int32 ra, int32 dec, double *x, double *y)
239 {
240 	double lat, lon;
241 	struct place pl;
242 
243 	lat = angle(dec);
244 	lon = angle(ra);
245 	pl.nlat.l = lat;
246 	pl.nlat.s = sin(lat);
247 	pl.nlat.c = cos(lat);
248 	pl.wlon.l = lon;
249 	pl.wlon.s = sin(lon);
250 	pl.wlon.c = cos(lon);
251 	normalize(&pl);
252 	return projection(&pl, x, y);
253 }
254 
255 /* end of 'heavens' section */
256 
257 int
258 setmap(int32 ramin, int32 ramax, int32 decmin, int32 decmax, Rectangle r, int zenithup)
259 {
260 	int c;
261 	double minx, maxx, miny, maxy;
262 
263 	c = 1000*60*60;
264 	mapra = ramax/2+ramin/2;
265 	mapdec = decmax/2+decmin/2;
266 
267 	if(zenithup)
268 		projection = heavens(mylat, mysid, (double)mapdec/c, (double)mapra/c);
269 	else{
270 		orient((double)mapdec/c, (double)mapra/c, 0.);
271 		projection = stereographic();
272 	}
273 	mapx0 = (r.max.x+r.min.x)/2;
274 	mapy0 = (r.max.y+r.min.y)/2;
275 	maps = ((double)Dy(r))/(double)(decmax-decmin);
276 	if(maptoxy(ramin, decmin, &minx, &miny) < 0)
277 		return -1;
278 	if(maptoxy(ramax, decmax, &maxx, &maxy) < 0)
279 		return -1;
280 	/*
281 	 * It's tricky to get the scale right.  This noble attempt scales
282 	 * to fit the lengths of the sides of the rectangle.
283 	 */
284 	mapscale = Dx(r)/(minx-maxx);
285 	if(mapscale > Dy(r)/(maxy-miny))
286 		mapscale = Dy(r)/(maxy-miny);
287 	/*
288 	 * near the pole It's not a rectangle, though, so this scales
289 	 * to fit the corners of the trapezoid (a triangle at the pole).
290 	 */
291 	mapscale *= (cos(angle(mapdec))+1.)/2;
292 	if(maxy  < miny){
293 		/* upside down, between zenith and pole */
294 		fprint(2, "reverse plot\n");
295 		mapscale = -mapscale;
296 	}
297 	return 1;
298 }
299 
300 Point
301 map(int32 ra, int32 dec)
302 {
303 	double x, y;
304 	Point p;
305 
306 	if(maptoxy(ra, dec, &x, &y) > 0){
307 		p.x = mapscale*x + mapx0;
308 		p.y = mapscale*-y + mapy0;
309 	}else{
310 		p.x = -100;
311 		p.y = -100;
312 	}
313 	return p;
314 }
315 
316 int
317 dsize(int mag)	/* mag is 10*magnitude; return disc size */
318 {
319 	double d;
320 
321 	mag += 25;	/* make mags all positive; sirius is -1.6m */
322 	d = (130-mag)/10;
323 	/* if plate scale is huge, adjust */
324 	if(maps < 100.0/MILLIARCSEC)
325 		d *= .71;
326 	if(maps < 50.0/MILLIARCSEC)
327 		d *= .71;
328 	return d;
329 }
330 
331 void
332 drawname(Image *scr, Image *col, char *s, int ra, int dec)
333 {
334 	Point p;
335 
336 	if(font == nil)
337 		return;
338 	p = addpt(map(ra, dec), Pt(4, -1));	/* font has huge ascent */
339 	string(scr, p, col, ZP, font, s);
340 }
341 
342 int
343 npixels(DAngle diam)
344 {
345 	Point p0, p1;
346 
347 	diam = DEG(angle(diam)*MILLIARCSEC);	/* diam in milliarcsec */
348 	diam /= 2;	/* radius */
349 	/* convert to plate scale */
350 	/* BUG: need +100 because we sometimes crash in library if we plot the exact center */
351 	p0 = map(mapra+100, mapdec);
352 	p1 = map(mapra+100+diam, mapdec);
353 	return hypot(p0.x-p1.x, p0.y-p1.y);
354 }
355 
356 void
357 drawdisc(Image *scr, DAngle semidiam, int ring, Image *color, Point pt, char *s)
358 {
359 	int d;
360 
361 	d = npixels(semidiam*2);
362 	if(d < 5)
363 		d = 5;
364 	fillellipse(scr, pt, d, d, color, ZP);
365 	if(ring == 1)
366 		ellipse(scr, pt, 5*d/2, d/2, 0, color, ZP);
367 	if(ring == 2)
368 		ellipse(scr, pt, d, d, 0, grey, ZP);
369 	if(s){
370 		d = stringwidth(font, s);
371 		pt.x -= d/2;
372 		pt.y -= font->height/2 + 1;
373 		string(scr, pt, display->black, ZP, font, s);
374 	}
375 }
376 
377 void
378 drawplanet(Image *scr, Planetrec *p, Point pt)
379 {
380 	if(strcmp(p->name, "sun") == 0){
381 		drawdisc(scr, p->semidiam, 0, suncolor, pt, nil);
382 		return;
383 	}
384 	if(strcmp(p->name, "moon") == 0){
385 		drawdisc(scr, p->semidiam, 0, mooncolor, pt, nil);
386 		return;
387 	}
388 	if(strcmp(p->name, "shadow") == 0){
389 		drawdisc(scr, p->semidiam, 2, shadowcolor, pt, nil);
390 		return;
391 	}
392 	if(strcmp(p->name, "mercury") == 0){
393 		drawdisc(scr, p->semidiam, 0, mercurycolor, pt, "m");
394 		return;
395 	}
396 	if(strcmp(p->name, "venus") == 0){
397 		drawdisc(scr, p->semidiam, 0, venuscolor, pt, "v");
398 		return;
399 	}
400 	if(strcmp(p->name, "mars") == 0){
401 		drawdisc(scr, p->semidiam, 0, marscolor, pt, "M");
402 		return;
403 	}
404 	if(strcmp(p->name, "jupiter") == 0){
405 		drawdisc(scr, p->semidiam, 0, jupitercolor, pt, "J");
406 		return;
407 	}
408 	if(strcmp(p->name, "saturn") == 0){
409 		drawdisc(scr, p->semidiam, 1, saturncolor, pt, "S");
410 
411 		return;
412 	}
413 	if(strcmp(p->name, "uranus") == 0){
414 		drawdisc(scr, p->semidiam, 0, uranuscolor, pt, "U");
415 
416 		return;
417 	}
418 	if(strcmp(p->name, "neptune") == 0){
419 		drawdisc(scr, p->semidiam, 0, neptunecolor, pt, "N");
420 
421 		return;
422 	}
423 	if(strcmp(p->name, "pluto") == 0){
424 		drawdisc(scr, p->semidiam, 0, plutocolor, pt, "P");
425 
426 		return;
427 	}
428 	if(strcmp(p->name, "comet") == 0){
429 		drawdisc(scr, p->semidiam, 0, cometcolor, pt, "C");
430 		return;
431 	}
432 
433 	pt.x -= stringwidth(font, p->name)/2;
434 	pt.y -= font->height/2;
435 	string(scr, pt, grey, ZP, font, p->name);
436 }
437 
438 void
439 tolast(char *name)
440 {
441 	int i, nlast;
442 	Record *r, rr;
443 
444 	/* stop early to simplify inner loop adjustment */
445 	nlast = 0;
446 	for(i=0,r=rec; i<nrec-nlast; i++,r++)
447 		if(r->type == Planet)
448 			if(name==nil || strcmp(r->u.planet.name, name)==0){
449 				rr = *r;
450 				memmove(rec+i, rec+i+1, (nrec-i-1)*sizeof(Record));
451 				rec[nrec-1] = rr;
452 				--i;
453 				--r;
454 				nlast++;
455 			}
456 }
457 
458 int
459 bbox(int32 extrara, int32 extradec, int quantize)
460 {
461 	int i;
462 	Record *r;
463 	int ra, dec;
464 	int rah, ram, d1, d2;
465 	double r0;
466 
467 	ramin = 0x7FFFFFFF;
468 	ramax = -0x7FFFFFFF;
469 	decmin = 0x7FFFFFFF;
470 	decmax = -0x7FFFFFFF;
471 
472 	for(i=0,r=rec; i<nrec; i++,r++){
473 		if(r->type == Patch){
474 			radec(r->index, &rah, &ram, &dec);
475 			ra = 15*rah+ram/4;
476 			r0 = c/cos(dec*PI/180);
477 			ra *= c;
478 			dec *= c;
479 			if(dec == 0)
480 				d1 = c, d2 = c;
481 			else if(dec < 0)
482 				d1 = c, d2 = 0;
483 			else
484 				d1 = 0, d2 = c;
485 		}else if(r->type==SAO || r->type==NGC || r->type==Planet || r->type==Abell){
486 			ra = r->u.ngc.ra;
487 			dec = r->u.ngc.dec;
488 			d1 = 0, d2 = 0, r0 = 0;
489 		}else
490 			continue;
491 		if(dec+d2+extradec > decmax)
492 			decmax = dec+d2+extradec;
493 		if(dec-d1-extradec < decmin)
494 			decmin = dec-d1-extradec;
495 		if(folded){
496 			ra -= 180*c;
497 			if(ra < 0)
498 				ra += 360*c;
499 		}
500 		if(ra+r0+extrara > ramax)
501 			ramax = ra+r0+extrara;
502 		if(ra-extrara < ramin)
503 			ramin = ra-extrara;
504 	}
505 
506 	if(decmax > 90*c)
507 		decmax = 90*c;
508 	if(decmin < -90*c)
509 		decmin = -90*c;
510 	if(ramin < 0)
511 		ramin += 360*c;
512 	if(ramax >= 360*c)
513 		ramax -= 360*c;
514 
515 	if(quantize){
516 		/* quantize to degree boundaries */
517 		ramin -= ramin%m5;
518 		if(ramax%m5 != 0)
519 			ramax += m5-(ramax%m5);
520 		if(decmin > 0)
521 			decmin -= decmin%c;
522 		else
523 			decmin -= c - (-decmin)%c;
524 		if(decmax > 0){
525 			if(decmax%c != 0)
526 				decmax += c-(decmax%c);
527 		}else if(decmax < 0){
528 			if(decmax%c != 0)
529 				decmax += ((-decmax)%c);
530 		}
531 	}
532 
533 	if(folded){
534 		if(ramax-ramin > 270*c){
535 			fprint(2, "ra range too wide %ld°\n", (ramax-ramin)/c);
536 			return -1;
537 		}
538 	}else if(ramax-ramin > 270*c){
539 		folded = 1;
540 		return bbox(0, 0, quantize);
541 	}
542 
543 	return 1;
544 }
545 
546 int
547 inbbox(DAngle ra, DAngle dec)
548 {
549 	int min;
550 
551 	if(dec < decmin)
552 		return 0;
553 	if(dec > decmax)
554 		return 0;
555 	min = ramin;
556 	if(ramin > ramax){	/* straddling 0h0m */
557 		min -= 360*c;
558 		if(ra > 180*c)
559 			ra -= 360*c;
560 	}
561 	if(ra < min)
562 		return 0;
563 	if(ra > ramax)
564 		return 0;
565 	return 1;
566 }
567 
568 int
569 gridra(int32 mapdec)
570 {
571 	mapdec = abs(mapdec)+c/2;
572 	mapdec /= c;
573 	if(mapdec < 60)
574 		return m5;
575 	if(mapdec < 80)
576 		return 2*m5;
577 	if(mapdec < 85)
578 		return 4*m5;
579 	return 8*m5;
580 }
581 
582 #define	GREY	(nogrey? display->white : grey)
583 
584 void
585 plot(char *flags)
586 {
587 	int i, j, k;
588 	char *t;
589 	int32 x, y;
590 	int ra, dec;
591 	int m;
592 	Point p, pts[10];
593 	Record *r;
594 	Rectangle rect, r1;
595 	int dx, dy, nogrid, textlevel, nogrey, zenithup;
596 	Image *scr;
597 	char *name, buf[32];
598 	double v;
599 
600 	if(plotopen() < 0)
601 		return;
602 	nogrid = 0;
603 	nogrey = 0;
604 	textlevel = 1;
605 	dx = 512;
606 	dy = 512;
607 	zenithup = 0;
608 	for(;;){
609 		if(t = alpha(flags, "nogrid")){
610 			nogrid = 1;
611 			flags = t;
612 			continue;
613 		}
614 		if((t = alpha(flags, "zenith")) || (t = alpha(flags, "zenithup")) ){
615 			zenithup = 1;
616 			flags = t;
617 			continue;
618 		}
619 		if((t = alpha(flags, "notext")) || (t = alpha(flags, "nolabel")) ){
620 			textlevel = 0;
621 			flags = t;
622 			continue;
623 		}
624 		if((t = alpha(flags, "alltext")) || (t = alpha(flags, "alllabel")) ){
625 			textlevel = 2;
626 			flags = t;
627 			continue;
628 		}
629 		if(t = alpha(flags, "dx")){
630 			dx = strtol(t, &t, 0);
631 			if(dx < 100){
632 				fprint(2, "dx %d too small (min 100) in plot\n", dx);
633 				return;
634 			}
635 			flags = skipbl(t);
636 			continue;
637 		}
638 		if(t = alpha(flags, "dy")){
639 			dy = strtol(t, &t, 0);
640 			if(dy < 100){
641 				fprint(2, "dy %d too small (min 100) in plot\n", dy);
642 				return;
643 			}
644 			flags = skipbl(t);
645 			continue;
646 		}
647 		if((t = alpha(flags, "nogrey")) || (t = alpha(flags, "nogray"))){
648 			nogrey = 1;
649 			flags = skipbl(t);
650 			continue;
651 		}
652 		if(*flags){
653 			fprint(2, "syntax error in plot\n");
654 			return;
655 		}
656 		break;
657 	}
658 	flatten();
659 	folded = 0;
660 
661 	if(bbox(0, 0, 1) < 0)
662 		return;
663 	if(ramax-ramin<100 || decmax-decmin<100){
664 		fprint(2, "plot too small\n");
665 		return;
666 	}
667 
668 	scr = allocimage(display, Rect(0, 0, dx, dy), RGB24, 0, DBlack);
669 	if(scr == nil){
670 		fprint(2, "can't allocate image: %r\n");
671 		return;
672 	}
673 	rect = scr->r;
674 	rect.min.x += 16;
675 	rect = insetrect(rect, 40);
676 	if(setmap(ramin, ramax, decmin, decmax, rect, zenithup) < 0){
677 		fprint(2, "can't set up map coordinates\n");
678 		return;
679 	}
680 	if(!nogrid){
681 		for(x=ramin; ; ){
682 			for(j=0; j<nelem(pts); j++){
683 				/* use double to avoid overflow */
684 				v = (double)j / (double)(nelem(pts)-1);
685 				v = decmin + v*(decmax-decmin);
686 				pts[j] = map(x, v);
687 			}
688 			bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
689 			ra = x;
690 			if(folded){
691 				ra -= 180*c;
692 				if(ra < 0)
693 					ra += 360*c;
694 			}
695 			p = pts[0];
696 			p.x -= stringwidth(font, hm5(angle(ra)))/2;
697 			string(scr, p, GREY, ZP, font, hm5(angle(ra)));
698 			p = pts[nelem(pts)-1];
699 			p.x -= stringwidth(font, hm5(angle(ra)))/2;
700 			p.y -= font->height;
701 			string(scr, p, GREY, ZP, font, hm5(angle(ra)));
702 			if(x == ramax)
703 				break;
704 			x += gridra(mapdec);
705 			if(x > ramax)
706 				x = ramax;
707 		}
708 		for(y=decmin; y<=decmax; y+=c){
709 			for(j=0; j<nelem(pts); j++){
710 				/* use double to avoid overflow */
711 				v = (double)j / (double)(nelem(pts)-1);
712 				v = ramin + v*(ramax-ramin);
713 				pts[j] = map(v, y);
714 			}
715 			bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
716 			p = pts[0];
717 			p.x += 3;
718 			p.y -= font->height/2;
719 			string(scr, p, GREY, ZP, font, deg(angle(y)));
720 			p = pts[nelem(pts)-1];
721 			p.x -= 3+stringwidth(font, deg(angle(y)));
722 			p.y -= font->height/2;
723 			string(scr, p, GREY, ZP, font, deg(angle(y)));
724 		}
725 	}
726 	/* reorder to get planets in front of stars */
727 	tolast(nil);
728 	tolast("moon");		/* moon is in front of everything... */
729 	tolast("shadow");	/* ... except the shadow */
730 
731 	for(i=0,r=rec; i<nrec; i++,r++){
732 		dec = r->u.ngc.dec;
733 		ra = r->u.ngc.ra;
734 		if(folded){
735 			ra -= 180*c;
736 			if(ra < 0)
737 				ra += 360*c;
738 		}
739 		if(textlevel){
740 			name = nameof(r);
741 			if(name==nil && textlevel>1 && r->type==SAO){
742 				snprint(buf, sizeof buf, "SAO%ld", r->index);
743 				name = buf;
744 			}
745 			if(name)
746 				drawname(scr, nogrey? display->white : alphagrey, name, ra, dec);
747 		}
748 		if(r->type == Planet){
749 			drawplanet(scr, &r->u.planet, map(ra, dec));
750 			continue;
751 		}
752 		if(r->type == SAO){
753 			m = r->u.sao.mag;
754 			if(m == UNKNOWNMAG)
755 				m = r->u.sao.mpg;
756 			if(m == UNKNOWNMAG)
757 				continue;
758 			m = dsize(m);
759 			if(m < 3)
760 				fillellipse(scr, map(ra, dec), m, m, nogrey? display->white : lightgrey, ZP);
761 			else{
762 				ellipse(scr, map(ra, dec), m+1, m+1, 0, display->black, ZP);
763 				fillellipse(scr, map(ra, dec), m, m, display->white, ZP);
764 			}
765 			continue;
766 		}
767 		if(r->type == Abell){
768 			ellipse(scr, addpt(map(ra, dec), Pt(-3, 2)), 2, 1, 0, lightblue, ZP);
769 			ellipse(scr, addpt(map(ra, dec), Pt(3, 2)), 2, 1, 0, lightblue, ZP);
770 			ellipse(scr, addpt(map(ra, dec), Pt(0, -2)), 1, 2, 0, lightblue, ZP);
771 			continue;
772 		}
773 		switch(r->u.ngc.type){
774 		case Galaxy:
775 			j = npixels(r->u.ngc.diam);
776 			if(j < 4)
777 				j = 4;
778 			if(j > 10)
779 				k = j/3;
780 			else
781 				k = j/2;
782 			ellipse(scr, map(ra, dec), j, k, 0, lightblue, ZP);
783 			break;
784 
785 		case PlanetaryN:
786 			p = map(ra, dec);
787 			j = npixels(r->u.ngc.diam);
788 			if(j < 3)
789 				j = 3;
790 			ellipse(scr, p, j, j, 0, green, ZP);
791 			line(scr, Pt(p.x, p.y+j+1), Pt(p.x, p.y+j+4),
792 				Endsquare, Endsquare, 0, green, ZP);
793 			line(scr, Pt(p.x, p.y-(j+1)), Pt(p.x, p.y-(j+4)),
794 				Endsquare, Endsquare, 0, green, ZP);
795 			line(scr, Pt(p.x+j+1, p.y), Pt(p.x+j+4, p.y),
796 				Endsquare, Endsquare, 0, green, ZP);
797 			line(scr, Pt(p.x-(j+1), p.y), Pt(p.x-(j+4), p.y),
798 				Endsquare, Endsquare, 0, green, ZP);
799 			break;
800 
801 		case DiffuseN:
802 		case NebularCl:
803 			p = map(ra, dec);
804 			j = npixels(r->u.ngc.diam);
805 			if(j < 4)
806 				j = 4;
807 			r1.min = Pt(p.x-j, p.y-j);
808 			r1.max = Pt(p.x+j, p.y+j);
809 			if(r->u.ngc.type != DiffuseN)
810 				draw(scr, r1, ocstipple, ocstipple, ZP);
811 			line(scr, Pt(p.x-j, p.y-j), Pt(p.x+j, p.y-j),
812 				Endsquare, Endsquare, 0, green, ZP);
813 			line(scr, Pt(p.x-j, p.y+j), Pt(p.x+j, p.y+j),
814 				Endsquare, Endsquare, 0, green, ZP);
815 			line(scr, Pt(p.x-j, p.y-j), Pt(p.x-j, p.y+j),
816 				Endsquare, Endsquare, 0, green, ZP);
817 			line(scr, Pt(p.x+j, p.y-j), Pt(p.x+j, p.y+j),
818 				Endsquare, Endsquare, 0, green, ZP);
819 			break;
820 
821 		case OpenCl:
822 			p = map(ra, dec);
823 			j = npixels(r->u.ngc.diam);
824 			if(j < 4)
825 				j = 4;
826 			fillellipse(scr, p, j, j, ocstipple, ZP);
827 			break;
828 
829 		case GlobularCl:
830 			j = npixels(r->u.ngc.diam);
831 			if(j < 4)
832 				j = 4;
833 			p = map(ra, dec);
834 			ellipse(scr, p, j, j, 0, lightgrey, ZP);
835 			line(scr, Pt(p.x-(j-1), p.y), Pt(p.x+j, p.y),
836 				Endsquare, Endsquare, 0, lightgrey, ZP);
837 			line(scr, Pt(p.x, p.y-(j-1)), Pt(p.x, p.y+j),
838 				Endsquare, Endsquare, 0, lightgrey, ZP);
839 			break;
840 
841 		}
842 	}
843 	flushimage(display, 1);
844 	displayimage(scr);
845 }
846 
847 int
848 runcommand(char *command, int p[2])
849 {
850 	switch(rfork(RFPROC|RFFDG|RFNOWAIT)){
851 	case -1:
852 		return -1;
853 	default:
854 		break;
855 	case 0:
856 		dup(p[1], 1);
857 		close(p[0]);
858 		close(p[1]);
859 		execlp("rc", "rc", "-c", command, nil);
860 		fprint(2, "can't exec %s: %r", command);
861 		exits(nil);
862 	}
863 	return 1;
864 }
865 
866 void
867 parseplanet(char *line, Planetrec *p)
868 {
869 	char *fld[6];
870 	int i, nfld;
871 	char *s;
872 
873 	if(line[0] == '\0')
874 		return;
875 	line[10] = '\0';	/* terminate name */
876 	s = strrchr(line, ' ');
877 	if(s == nil)
878 		s = line;
879 	else
880 		s++;
881 	strcpy(p->name, s);
882 	for(i=0; s[i]!='\0'; i++)
883 		p->name[i] = tolower(s[i]);
884 	p->name[i] = '\0';
885 	nfld = getfields(line+11, fld, nelem(fld), 1, " ");
886 	p->ra = dangle(getra(fld[0]));
887 	p->dec = dangle(getra(fld[1]));
888 	p->az = atof(fld[2])*MILLIARCSEC;
889 	p->alt = atof(fld[3])*MILLIARCSEC;
890 	p->semidiam = atof(fld[4])*1000;
891 	if(nfld > 5)
892 		p->phase = atof(fld[5]);
893 	else
894 		p->phase = 0;
895 }
896 
897 void
898 astro(char *flags, int initial)
899 {
900 	int p[2];
901 	int i, n, np;
902 	char cmd[256], buf[4096], *lines[20], *fld[10];
903 
904 	snprint(cmd, sizeof cmd, "astro -p %s", flags);
905 	if(pipe(p) < 0){
906 		fprint(2, "can't pipe: %r\n");
907 		return;
908 	}
909 	if(runcommand(cmd, p) < 0){
910 		close(p[0]);
911 		close(p[1]);
912 		fprint(2, "can't run astro: %r");
913 		return;
914 	}
915 	close(p[1]);
916 	n = readn(p[0], buf, sizeof buf-1);
917 	if(n <= 0){
918 		fprint(2, "no data from astro\n");
919 		return;
920 	}
921 	if(!initial)
922 		Bwrite(&bout, buf, n);
923 	buf[n] = '\0';
924 	np = getfields(buf, lines, nelem(lines), 0, "\n");
925 	if(np <= 1){
926 		fprint(2, "astro: not enough output\n");
927 		return;
928 	}
929 	Bprint(&bout, "%s\n", lines[0]);
930 	Bflush(&bout);
931 	/* get latitude and longitude */
932 	if(getfields(lines[0], fld, nelem(fld), 1, " ") < 8)
933 		fprint(2, "astro: can't read longitude: too few fields\n");
934 	else{
935 		mysid = getra(fld[5])*180./PI;
936 		mylat = getra(fld[6])*180./PI;
937 		mylon = getra(fld[7])*180./PI;
938 	}
939 	/*
940 	 * Each time we run astro, we generate a new planet list
941 	 * so multiple appearances of a planet may exist as we plot
942 	 * its motion over time.
943 	 */
944 	planet = malloc(NPlanet*sizeof planet[0]);
945 	if(planet == nil){
946 		fprint(2, "astro: malloc failed: %r\n");
947 		exits("malloc");
948 	}
949 	memset(planet, 0, NPlanet*sizeof planet[0]);
950 	for(i=1; i<np; i++)
951 		parseplanet(lines[i], &planet[i-1]);
952 }