Subversion Repositories planix.SVN

Rev

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

#include <u.h>
#include <libc.h>
#include <bio.h>
#include <draw.h>
#include <event.h>
#include <ctype.h>
#include "map.h"
#undef  RAD
#undef  TWOPI
#include "sky.h"

        /* convert to milliarcsec */
static long     c = MILLIARCSEC;        /* 1 degree */
static long     m5 = 1250*60*60;        /* 5 minutes of ra */

DAngle  ramin;
DAngle  ramax;
DAngle  decmin;
DAngle  decmax;
int             folded;

Image   *grey;
Image   *alphagrey;
Image   *green;
Image   *lightblue;
Image   *lightgrey;
Image   *ocstipple;
Image   *suncolor;
Image   *mooncolor;
Image   *shadowcolor;
Image   *mercurycolor;
Image   *venuscolor;
Image   *marscolor;
Image   *jupitercolor;
Image   *saturncolor;
Image   *uranuscolor;
Image   *neptunecolor;
Image   *plutocolor;
Image   *cometcolor;

Planetrec       *planet;

long    mapx0, mapy0;
long    mapra, mapdec;
double  mylat, mylon, mysid;
double  mapscale;
double  maps;
int (*projection)(struct place*, double*, double*);

char *fontname = "/lib/font/bit/lucida/unicode.6.font";

/* types Coord and Loc correspond to types in map(3) thus:
   Coord == struct coord;
   Loc == struct place, except longitudes are measured
     positive east for Loc and west for struct place
*/

typedef struct Xyz Xyz;
typedef struct coord Coord;
typedef struct Loc Loc;

struct Xyz{
        double x,y,z;
};

struct Loc{
        Coord nlat, elon;
};

Xyz north = { 0, 0, 1 };

int
plotopen(void)
{
        if(display != nil)
                return 1;
        display = initdisplay(nil, nil, drawerror);
        if(display == nil){
                fprint(2, "initdisplay failed: %r\n");
                return -1;
        }
        grey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x777777FF);
        lightgrey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0xAAAAAAFF);
        alphagrey = allocimage(display, Rect(0, 0, 2, 2), RGBA32, 1, 0x77777777);
        green = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x00AA00FF);
        lightblue = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x009EEEFF);
        ocstipple = allocimage(display, Rect(0, 0, 2, 2), CMAP8, 1, 0xAAAAAAFF);
        draw(ocstipple, Rect(0, 0, 1, 1), display->black, nil, ZP);
        draw(ocstipple, Rect(1, 1, 2, 2), display->black, nil, ZP);

        suncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFF77FF);
        mooncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAAAFF);
        shadowcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x00000055);
        mercurycolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFAAAAFF);
        venuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);
        marscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFF5555FF);
        jupitercolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFFAAFF);
        saturncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAFFAAFF);
        uranuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77DDDDFF);
        neptunecolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77FF77FF);
        plutocolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x7777FFFF);
        cometcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);
        font = openfont(display, fontname);
        if(font == nil)
                fprint(2, "warning: no font %s: %r\n", fontname);
        return 1;
}

/*
 * Function heavens() for setting up observer-eye-view
 * sky charts, plus subsidiary functions.
 * Written by Doug McIlroy.
 */

/* heavens(nlato, wlono, nlatc, wlonc) sets up a map(3)-style
   coordinate change (by calling orient()) and returns a
   projection function heavens for calculating a star map
   centered on nlatc,wlonc and rotated so it will appear
   upright as seen by a ground observer at nlato,wlono.
   all coordinates (north latitude and west longitude)
   are in degrees.
*/

Coord
mkcoord(double degrees)
{
        Coord c;

        c.l = degrees*PI/180;
        c.s = sin(c.l);
        c.c = cos(c.l);
        return c;
}

Xyz
ptov(struct place p)
{
        Xyz v;

        v.z = p.nlat.s;
        v.x = p.nlat.c * p.wlon.c;
        v.y = -p.nlat.c * p.wlon.s;
        return v;
}

double
dot(Xyz a, Xyz b)
{
        return a.x*b.x + a.y*b.y + a.z*b.z;
}

Xyz
cross(Xyz a, Xyz b)
{
        Xyz v;

        v.x = a.y*b.z - a.z*b.y;
        v.y = a.z*b.x - a.x*b.z;
        v.z = a.x*b.y - a.y*b.x;
        return v;
}

double
len(Xyz a)
{
        return sqrt(dot(a, a));
}

/* An azimuth vector from a to b is a tangent to
   the sphere at a pointing along the (shorter)
   great-circle direction to b.  It lies in the
   plane of the vectors a and b, is perpendicular
   to a, and has a positive compoent along b,
   provided a and b span a 2-space.  Otherwise
   it is null.  (aXb)Xa, where X denotes cross
   product, is such a vector. */

Xyz
azvec(Xyz a, Xyz b)
{
        return cross(cross(a,b), a);
}

/* Find the azimuth of point b as seen
   from point a, given that a is distinct
   from either pole and from b */

double
azimuth(Xyz a, Xyz b)
{
        Xyz aton = azvec(a,north);
        Xyz atob = azvec(a,b);
        double lenaton = len(aton);
        double lenatob = len(atob);
        double az = acos(dot(aton,atob)/(lenaton*lenatob));

        if(dot(b,cross(north,a)) < 0)
                az = -az;
        return az;
}

/* Find the rotation (3rd argument of orient() in the
   map projection package) for the map described
   below.  The return is radians; it must be converted
   to degrees for orient().
*/

double
rot(struct place center, struct place zenith)
{
        Xyz cen = ptov(center);
        Xyz zen = ptov(zenith);

        if(cen.z > 1-FUZZ)      /* center at N pole */
                return PI + zenith.wlon.l;
        if(cen.z < FUZZ-1)              /* at S pole */
                return -zenith.wlon.l;
        if(fabs(dot(cen,zen)) > 1-FUZZ) /* at zenith */
                return 0;
        return azimuth(cen, zen);
}

int (*
heavens(double zlatdeg, double zlondeg, double clatdeg, double clondeg))(Loc*, double*, double*)
{
        struct place center;
        struct place zenith;

        center.nlat = mkcoord(clatdeg);
        center.wlon = mkcoord(clondeg);
        zenith.nlat = mkcoord(zlatdeg);
        zenith.wlon = mkcoord(zlondeg);
        projection = stereographic();
        orient(clatdeg, clondeg, rot(center, zenith)*180/PI);
        return projection;
}

int
maptoxy(long ra, long dec, double *x, double *y)
{
        double lat, lon;
        struct place pl;

        lat = angle(dec);
        lon = angle(ra);
        pl.nlat.l = lat;
        pl.nlat.s = sin(lat);
        pl.nlat.c = cos(lat);
        pl.wlon.l = lon;
        pl.wlon.s = sin(lon);
        pl.wlon.c = cos(lon);
        normalize(&pl);
        return projection(&pl, x, y);
}

/* end of 'heavens' section */

int
setmap(long ramin, long ramax, long decmin, long decmax, Rectangle r, int zenithup)
{
        int c;
        double minx, maxx, miny, maxy;

        c = 1000*60*60;
        mapra = ramax/2+ramin/2;
        mapdec = decmax/2+decmin/2;

        if(zenithup)
                projection = heavens(mylat, mysid, (double)mapdec/c, (double)mapra/c);
        else{
                orient((double)mapdec/c, (double)mapra/c, 0.);
                projection = stereographic();
        }
        mapx0 = (r.max.x+r.min.x)/2;
        mapy0 = (r.max.y+r.min.y)/2;
        maps = ((double)Dy(r))/(double)(decmax-decmin);
        if(maptoxy(ramin, decmin, &minx, &miny) < 0)
                return -1;
        if(maptoxy(ramax, decmax, &maxx, &maxy) < 0)
                return -1;
        /*
         * It's tricky to get the scale right.  This noble attempt scales
         * to fit the lengths of the sides of the rectangle.
         */
        mapscale = Dx(r)/(minx-maxx);
        if(mapscale > Dy(r)/(maxy-miny))
                mapscale = Dy(r)/(maxy-miny);
        /*
         * near the pole It's not a rectangle, though, so this scales
         * to fit the corners of the trapezoid (a triangle at the pole).
         */
        mapscale *= (cos(angle(mapdec))+1.)/2;
        if(maxy  < miny){
                /* upside down, between zenith and pole */
                fprint(2, "reverse plot\n");
                mapscale = -mapscale;
        }
        return 1;
}

Point
map(long ra, long dec)
{
        double x, y;
        Point p;

        if(maptoxy(ra, dec, &x, &y) > 0){
                p.x = mapscale*x + mapx0;
                p.y = mapscale*-y + mapy0;
        }else{
                p.x = -100;
                p.y = -100;
        }
        return p;
}

int
dsize(int mag)  /* mag is 10*magnitude; return disc size */
{
        double d;

        mag += 25;      /* make mags all positive; sirius is -1.6m */
        d = (130-mag)/10;
        /* if plate scale is huge, adjust */
        if(maps < 100.0/MILLIARCSEC)
                d *= .71;
        if(maps < 50.0/MILLIARCSEC)
                d *= .71;
        return d;
}

void
drawname(Image *scr, Image *col, char *s, int ra, int dec)
{
        Point p;

        if(font == nil)
                return;
        p = addpt(map(ra, dec), Pt(4, -1));     /* font has huge ascent */
        string(scr, p, col, ZP, font, s);
}

int
npixels(DAngle diam)
{
        Point p0, p1;

        diam = DEG(angle(diam)*MILLIARCSEC);    /* diam in milliarcsec */
        diam /= 2;      /* radius */
        /* convert to plate scale */
        /* BUG: need +100 because we sometimes crash in library if we plot the exact center */
        p0 = map(mapra+100, mapdec);
        p1 = map(mapra+100+diam, mapdec);
        return hypot(p0.x-p1.x, p0.y-p1.y);
}

void
drawdisc(Image *scr, DAngle semidiam, int ring, Image *color, Point pt, char *s)
{
        int d;

        d = npixels(semidiam*2);
        if(d < 5)
                d = 5;
        fillellipse(scr, pt, d, d, color, ZP);
        if(ring == 1)
                ellipse(scr, pt, 5*d/2, d/2, 0, color, ZP);
        if(ring == 2)
                ellipse(scr, pt, d, d, 0, grey, ZP);
        if(s){
                d = stringwidth(font, s);
                pt.x -= d/2;
                pt.y -= font->height/2 + 1;
                string(scr, pt, display->black, ZP, font, s);
        }
}

void
drawplanet(Image *scr, Planetrec *p, Point pt)
{
        if(strcmp(p->name, "sun") == 0){
                drawdisc(scr, p->semidiam, 0, suncolor, pt, nil);
                return;
        }
        if(strcmp(p->name, "moon") == 0){
                drawdisc(scr, p->semidiam, 0, mooncolor, pt, nil);
                return;
        }
        if(strcmp(p->name, "shadow") == 0){
                drawdisc(scr, p->semidiam, 2, shadowcolor, pt, nil);
                return;
        }
        if(strcmp(p->name, "mercury") == 0){
                drawdisc(scr, p->semidiam, 0, mercurycolor, pt, "m");
                return;
        }
        if(strcmp(p->name, "venus") == 0){
                drawdisc(scr, p->semidiam, 0, venuscolor, pt, "v");
                return;
        }
        if(strcmp(p->name, "mars") == 0){
                drawdisc(scr, p->semidiam, 0, marscolor, pt, "M");
                return;
        }
        if(strcmp(p->name, "jupiter") == 0){
                drawdisc(scr, p->semidiam, 0, jupitercolor, pt, "J");
                return;
        }
        if(strcmp(p->name, "saturn") == 0){
                drawdisc(scr, p->semidiam, 1, saturncolor, pt, "S");
                
                return;
        }
        if(strcmp(p->name, "uranus") == 0){
                drawdisc(scr, p->semidiam, 0, uranuscolor, pt, "U");
                
                return;
        }
        if(strcmp(p->name, "neptune") == 0){
                drawdisc(scr, p->semidiam, 0, neptunecolor, pt, "N");
                
                return;
        }
        if(strcmp(p->name, "pluto") == 0){
                drawdisc(scr, p->semidiam, 0, plutocolor, pt, "P");
                
                return;
        }
        if(strcmp(p->name, "comet") == 0){
                drawdisc(scr, p->semidiam, 0, cometcolor, pt, "C");
                return;
        }

        pt.x -= stringwidth(font, p->name)/2;
        pt.y -= font->height/2;
        string(scr, pt, grey, ZP, font, p->name);
}

void
tolast(char *name)
{
        int i, nlast;
        Record *r, rr;

        /* stop early to simplify inner loop adjustment */
        nlast = 0;
        for(i=0,r=rec; i<nrec-nlast; i++,r++)
                if(r->type == Planet)
                        if(name==nil || strcmp(r->planet.name, name)==0){
                                rr = *r;
                                memmove(rec+i, rec+i+1, (nrec-i-1)*sizeof(Record));
                                rec[nrec-1] = rr;
                                --i;
                                --r;
                                nlast++;
                        }
}

int
bbox(long extrara, long extradec, int quantize)
{
        int i;
        Record *r;
        int ra, dec;
        int rah, ram, d1, d2;
        double r0;

        ramin = 0x7FFFFFFF;
        ramax = -0x7FFFFFFF;
        decmin = 0x7FFFFFFF;
        decmax = -0x7FFFFFFF;

        for(i=0,r=rec; i<nrec; i++,r++){
                if(r->type == Patch){
                        radec(r->index, &rah, &ram, &dec);
                        ra = 15*rah+ram/4;
                        r0 = c/cos(dec*PI/180);
                        ra *= c;
                        dec *= c;
                        if(dec == 0)
                                d1 = c, d2 = c;
                        else if(dec < 0)
                                d1 = c, d2 = 0;
                        else
                                d1 = 0, d2 = c;
                }else if(r->type==SAO || r->type==NGC || r->type==Planet || r->type==Abell){
                        ra = r->ngc.ra;
                        dec = r->ngc.dec;
                        d1 = 0, d2 = 0, r0 = 0;
                }else
                        continue;
                if(dec+d2+extradec > decmax)
                        decmax = dec+d2+extradec;
                if(dec-d1-extradec < decmin)
                        decmin = dec-d1-extradec;
                if(folded){
                        ra -= 180*c;
                        if(ra < 0)
                                ra += 360*c;
                }
                if(ra+r0+extrara > ramax)
                        ramax = ra+r0+extrara;
                if(ra-extrara < ramin)
                        ramin = ra-extrara;
        }

        if(decmax > 90*c)
                decmax = 90*c;
        if(decmin < -90*c)
                decmin = -90*c;
        if(ramin < 0)
                ramin += 360*c;
        if(ramax >= 360*c)
                ramax -= 360*c;

        if(quantize){
                /* quantize to degree boundaries */
                ramin -= ramin%m5;
                if(ramax%m5 != 0)
                        ramax += m5-(ramax%m5);
                if(decmin > 0)
                        decmin -= decmin%c;
                else
                        decmin -= c - (-decmin)%c;
                if(decmax > 0){
                        if(decmax%c != 0)
                                decmax += c-(decmax%c);
                }else if(decmax < 0){
                        if(decmax%c != 0)
                                decmax += ((-decmax)%c);
                }
        }

        if(folded){
                if(ramax-ramin > 270*c){
                        fprint(2, "ra range too wide %ldĀ°\n", (ramax-ramin)/c);
                        return -1;
                }
        }else if(ramax-ramin > 270*c){
                folded = 1;
                return bbox(0, 0, quantize);
        }

        return 1;
}

int
inbbox(DAngle ra, DAngle dec)
{
        int min;

        if(dec < decmin)
                return 0;
        if(dec > decmax)
                return 0;
        min = ramin;
        if(ramin > ramax){      /* straddling 0h0m */
                min -= 360*c;
                if(ra > 180*c)
                        ra -= 360*c;
        }
        if(ra < min)
                return 0;
        if(ra > ramax)
                return 0;
        return 1;
}

int
gridra(long mapdec)
{
        mapdec = abs(mapdec)+c/2;
        mapdec /= c;
        if(mapdec < 60)
                return m5;
        if(mapdec < 80)
                return 2*m5;
        if(mapdec < 85)
                return 4*m5;
        return 8*m5;
}

#define GREY    (nogrey? display->white : grey)

void
plot(char *flags)
{
        int i, j, k;
        char *t;
        long x, y;
        int ra, dec;
        int m;
        Point p, pts[10];
        Record *r;
        Rectangle rect, r1;
        int dx, dy, nogrid, textlevel, nogrey, zenithup;
        Image *scr;
        char *name, buf[32];
        double v;

        if(plotopen() < 0)
                return;
        nogrid = 0;
        nogrey = 0;
        textlevel = 1;
        dx = 512;
        dy = 512;
        zenithup = 0;
        for(;;){
                if(t = alpha(flags, "nogrid")){
                        nogrid = 1;
                        flags = t;
                        continue;
                }
                if((t = alpha(flags, "zenith")) || (t = alpha(flags, "zenithup")) ){
                        zenithup = 1;
                        flags = t;
                        continue;
                }
                if((t = alpha(flags, "notext")) || (t = alpha(flags, "nolabel")) ){
                        textlevel = 0;
                        flags = t;
                        continue;
                }
                if((t = alpha(flags, "alltext")) || (t = alpha(flags, "alllabel")) ){
                        textlevel = 2;
                        flags = t;
                        continue;
                }
                if(t = alpha(flags, "dx")){
                        dx = strtol(t, &t, 0);
                        if(dx < 100){
                                fprint(2, "dx %d too small (min 100) in plot\n", dx);
                                return;
                        }
                        flags = skipbl(t);
                        continue;
                }
                if(t = alpha(flags, "dy")){
                        dy = strtol(t, &t, 0);
                        if(dy < 100){
                                fprint(2, "dy %d too small (min 100) in plot\n", dy);
                                return;
                        }
                        flags = skipbl(t);
                        continue;
                }
                if((t = alpha(flags, "nogrey")) || (t = alpha(flags, "nogray"))){
                        nogrey = 1;
                        flags = skipbl(t);
                        continue;
                }
                if(*flags){
                        fprint(2, "syntax error in plot\n");
                        return;
                }
                break;
        }
        flatten();
        folded = 0;

        if(bbox(0, 0, 1) < 0)
                return;
        if(ramax-ramin<100 || decmax-decmin<100){
                fprint(2, "plot too small\n");
                return;
        }

        scr = allocimage(display, Rect(0, 0, dx, dy), RGB24, 0, DBlack);
        if(scr == nil){
                fprint(2, "can't allocate image: %r\n");
                return;
        }
        rect = scr->r;
        rect.min.x += 16;
        rect = insetrect(rect, 40);
        if(setmap(ramin, ramax, decmin, decmax, rect, zenithup) < 0){
                fprint(2, "can't set up map coordinates\n");
                return;
        }
        if(!nogrid){
                for(x=ramin; ; ){
                        for(j=0; j<nelem(pts); j++){
                                /* use double to avoid overflow */
                                v = (double)j / (double)(nelem(pts)-1);
                                v = decmin + v*(decmax-decmin);
                                pts[j] = map(x, v);
                        }
                        bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
                        ra = x;
                        if(folded){
                                ra -= 180*c;
                                if(ra < 0)
                                        ra += 360*c;
                        }
                        p = pts[0];
                        p.x -= stringwidth(font, hm5(angle(ra)))/2;
                        string(scr, p, GREY, ZP, font, hm5(angle(ra)));
                        p = pts[nelem(pts)-1];
                        p.x -= stringwidth(font, hm5(angle(ra)))/2;
                        p.y -= font->height;
                        string(scr, p, GREY, ZP, font, hm5(angle(ra)));
                        if(x == ramax)
                                break;
                        x += gridra(mapdec);
                        if(x > ramax)
                                x = ramax;
                }
                for(y=decmin; y<=decmax; y+=c){
                        for(j=0; j<nelem(pts); j++){
                                /* use double to avoid overflow */
                                v = (double)j / (double)(nelem(pts)-1);
                                v = ramin + v*(ramax-ramin);
                                pts[j] = map(v, y);
                        }
                        bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
                        p = pts[0];
                        p.x += 3;
                        p.y -= font->height/2;
                        string(scr, p, GREY, ZP, font, deg(angle(y)));
                        p = pts[nelem(pts)-1];
                        p.x -= 3+stringwidth(font, deg(angle(y)));
                        p.y -= font->height/2;
                        string(scr, p, GREY, ZP, font, deg(angle(y)));
                }
        }
        /* reorder to get planets in front of stars */
        tolast(nil);
        tolast("moon");         /* moon is in front of everything... */
        tolast("shadow");       /* ... except the shadow */

        for(i=0,r=rec; i<nrec; i++,r++){
                dec = r->ngc.dec;
                ra = r->ngc.ra;
                if(folded){
                        ra -= 180*c;
                        if(ra < 0)
                                ra += 360*c;
                }
                if(textlevel){
                        name = nameof(r);
                        if(name==nil && textlevel>1 && r->type==SAO){
                                snprint(buf, sizeof buf, "SAO%ld", r->index);
                                name = buf;
                        }
                        if(name)
                                drawname(scr, nogrey? display->white : alphagrey, name, ra, dec);
                }
                if(r->type == Planet){
                        drawplanet(scr, &r->planet, map(ra, dec));
                        continue;
                }
                if(r->type == SAO){
                        m = r->sao.mag;
                        if(m == UNKNOWNMAG)
                                m = r->sao.mpg;
                        if(m == UNKNOWNMAG)
                                continue;
                        m = dsize(m);
                        if(m < 3)
                                fillellipse(scr, map(ra, dec), m, m, nogrey? display->white : lightgrey, ZP);
                        else{
                                ellipse(scr, map(ra, dec), m+1, m+1, 0, display->black, ZP);
                                fillellipse(scr, map(ra, dec), m, m, display->white, ZP);
                        }
                        continue;
                }
                if(r->type == Abell){
                        ellipse(scr, addpt(map(ra, dec), Pt(-3, 2)), 2, 1, 0, lightblue, ZP);
                        ellipse(scr, addpt(map(ra, dec), Pt(3, 2)), 2, 1, 0, lightblue, ZP);
                        ellipse(scr, addpt(map(ra, dec), Pt(0, -2)), 1, 2, 0, lightblue, ZP);
                        continue;
                }
                switch(r->ngc.type){
                case Galaxy:
                        j = npixels(r->ngc.diam);
                        if(j < 4)
                                j = 4;
                        if(j > 10)
                                k = j/3;
                        else
                                k = j/2;
                        ellipse(scr, map(ra, dec), j, k, 0, lightblue, ZP);
                        break;

                case PlanetaryN:
                        p = map(ra, dec);
                        j = npixels(r->ngc.diam);
                        if(j < 3)
                                j = 3;
                        ellipse(scr, p, j, j, 0, green, ZP);
                        line(scr, Pt(p.x, p.y+j+1), Pt(p.x, p.y+j+4),
                                Endsquare, Endsquare, 0, green, ZP);
                        line(scr, Pt(p.x, p.y-(j+1)), Pt(p.x, p.y-(j+4)),
                                Endsquare, Endsquare, 0, green, ZP);
                        line(scr, Pt(p.x+j+1, p.y), Pt(p.x+j+4, p.y),
                                Endsquare, Endsquare, 0, green, ZP);
                        line(scr, Pt(p.x-(j+1), p.y), Pt(p.x-(j+4), p.y),
                                Endsquare, Endsquare, 0, green, ZP);
                        break;

                case DiffuseN:
                case NebularCl:
                        p = map(ra, dec);
                        j = npixels(r->ngc.diam);
                        if(j < 4)
                                j = 4;
                        r1.min = Pt(p.x-j, p.y-j);
                        r1.max = Pt(p.x+j, p.y+j);
                        if(r->ngc.type != DiffuseN)
                                draw(scr, r1, ocstipple, ocstipple, ZP);
                        line(scr, Pt(p.x-j, p.y-j), Pt(p.x+j, p.y-j),
                                Endsquare, Endsquare, 0, green, ZP);
                        line(scr, Pt(p.x-j, p.y+j), Pt(p.x+j, p.y+j),
                                Endsquare, Endsquare, 0, green, ZP);
                        line(scr, Pt(p.x-j, p.y-j), Pt(p.x-j, p.y+j),
                                Endsquare, Endsquare, 0, green, ZP);
                        line(scr, Pt(p.x+j, p.y-j), Pt(p.x+j, p.y+j),
                                Endsquare, Endsquare, 0, green, ZP);
                        break;

                case OpenCl:
                        p = map(ra, dec);
                        j = npixels(r->ngc.diam);
                        if(j < 4)
                                j = 4;
                        fillellipse(scr, p, j, j, ocstipple, ZP);
                        break;

                case GlobularCl:
                        j = npixels(r->ngc.diam);
                        if(j < 4)
                                j = 4;
                        p = map(ra, dec);
                        ellipse(scr, p, j, j, 0, lightgrey, ZP);
                        line(scr, Pt(p.x-(j-1), p.y), Pt(p.x+j, p.y),
                                Endsquare, Endsquare, 0, lightgrey, ZP);
                        line(scr, Pt(p.x, p.y-(j-1)), Pt(p.x, p.y+j),
                                Endsquare, Endsquare, 0, lightgrey, ZP);
                        break;

                }
        }
        flushimage(display, 1);
        displayimage(scr);
}

int
runcommand(char *command, int p[2])
{
        switch(rfork(RFPROC|RFFDG|RFNOWAIT)){
        case -1:
                return -1;
        default:
                break;
        case 0:
                dup(p[1], 1);
                close(p[0]);
                close(p[1]);
                execl("/bin/rc", "rc", "-c", command, nil);
                fprint(2, "can't exec %s: %r", command);
                exits(nil);
        }
        return 1;
}

void
parseplanet(char *line, Planetrec *p)
{
        char *fld[6];
        int i, nfld;
        char *s;

        if(line[0] == '\0')
                return;
        line[10] = '\0';        /* terminate name */
        s = strrchr(line, ' ');
        if(s == nil)
                s = line;
        else
                s++;
        strcpy(p->name, s);
        for(i=0; s[i]!='\0'; i++)
                p->name[i] = tolower(s[i]);
        p->name[i] = '\0';
        nfld = getfields(line+11, fld, nelem(fld), 1, " ");
        p->ra = dangle(getra(fld[0]));
        p->dec = dangle(getra(fld[1]));
        p->az = atof(fld[2])*MILLIARCSEC;
        p->alt = atof(fld[3])*MILLIARCSEC;
        p->semidiam = atof(fld[4])*1000;
        if(nfld > 5)
                p->phase = atof(fld[5]);
        else
                p->phase = 0;
}

void
astro(char *flags, int initial)
{
        int p[2];
        int i, n, np;
        char cmd[256], buf[4096], *lines[20], *fld[10];

        snprint(cmd, sizeof cmd, "/bin/astro -p %s", flags);
        if(pipe(p) < 0){
                fprint(2, "can't pipe: %r\n");
                return;
        }
        if(runcommand(cmd, p) < 0){
                close(p[0]);
                close(p[1]);
                fprint(2, "can't run astro: %r");
                return;
        }
        close(p[1]);
        n = readn(p[0], buf, sizeof buf-1);
        if(n <= 0){
                fprint(2, "no data from astro\n");
                return;
        }
        if(!initial)
                Bwrite(&bout, buf, n);
        buf[n] = '\0';
        np = getfields(buf, lines, nelem(lines), 0, "\n");
        if(np <= 1){
                fprint(2, "astro: not enough output\n");
                return;
        }
        Bprint(&bout, "%s\n", lines[0]);
        Bflush(&bout);
        /* get latitude and longitude */
        if(getfields(lines[0], fld, nelem(fld), 1, " ") < 8)
                fprint(2, "astro: can't read longitude: too few fields\n");
        else{
                mysid = getra(fld[5])*180./PI;
                mylat = getra(fld[6])*180./PI;
                mylon = getra(fld[7])*180./PI;
        }
        /*
         * Each time we run astro, we generate a new planet list
         * so multiple appearances of a planet may exist as we plot
         * its motion over time.
         */
        planet = malloc(NPlanet*sizeof planet[0]);
        if(planet == nil){
                fprint(2, "astro: malloc failed: %r\n");
                exits("malloc");
        }
        memset(planet, 0, NPlanet*sizeof planet[0]);
        for(i=1; i<np; i++)
                parseplanet(lines[i], &planet[i-1]);
}