123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953 |
- #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]);
- }
|