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
plotopen(void)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
mkcoord(double degrees)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
ptov(struct place p)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
dot(Xyz a,Xyz b)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
cross(Xyz a,Xyz b)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
len(Xyz a)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
azvec(Xyz a,Xyz b)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
azimuth(Xyz a,Xyz b)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
rot(struct place center,struct place zenith)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 (*
heavens(double zlatdeg,double zlondeg,double clatdeg,double clondeg)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
maptoxy(int32 ra,int32 dec,double * x,double * y)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
setmap(int32 ramin,int32 ramax,int32 decmin,int32 decmax,Rectangle r,int zenithup)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
map(int32 ra,int32 dec)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
dsize(int mag)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
drawname(Image * scr,Image * col,char * s,int ra,int dec)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
npixels(DAngle diam)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
drawdisc(Image * scr,DAngle semidiam,int ring,Image * color,Point pt,char * s)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
drawplanet(Image * scr,Planetrec * p,Point pt)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
tolast(char * name)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
bbox(int32 extrara,int32 extradec,int quantize)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
inbbox(DAngle ra,DAngle dec)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
gridra(int32 mapdec)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
plot(char * flags)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
runcommand(char * command,int p[2])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
parseplanet(char * line,Planetrec * p)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
astro(char * flags,int initial)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 }
953