sys/src/cmd/scat/plot.c

   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 }