Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

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