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 } 954