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