12 /* convert to milliarcsec */
13 static long c = MILLIARCSEC; /* 1 degree */
14 static long m5 = 1250*60*60; /* 5 minutes of ra */
45 double mylat, mylon, mysid;
48 int (*projection)(struct place*, double*, double*);
50 char *fontname = "/lib/font/bit/luc/unicode.6.font";
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
58 typedef struct Xyz Xyz;
59 typedef struct coord Coord;
60 typedef struct Loc Loc;
70 Xyz north = { 0, 0, 1 };
77 if(initdraw(drawerror, nil, nil) < 0){
78 fprint(2, "initdisplay failed: %r\n");
81 grey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x777777FF);
82 lightgrey = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0xAAAAAAFF);
83 alphagrey = allocimage(display, Rect(0, 0, 2, 2), RGBA32, 1, 0x77777777);
84 green = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x00AA00FF);
85 lightblue = allocimage(display, Rect(0, 0, 1, 1), CMAP8, 1, 0x009EEEFF);
86 ocstipple = allocimage(display, Rect(0, 0, 2, 2), CMAP8, 1, 0xAAAAAAFF);
87 draw(ocstipple, Rect(0, 0, 1, 1), display->black, nil, ZP);
88 draw(ocstipple, Rect(1, 1, 2, 2), display->black, nil, ZP);
90 suncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFF77FF);
91 mooncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAAAFF);
92 shadowcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x00000055);
93 mercurycolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFAAAAFF);
94 venuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);
95 marscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFF5555FF);
96 jupitercolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xFFFFAAFF);
97 saturncolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAFFAAFF);
98 uranuscolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77DDDDFF);
99 neptunecolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x77FF77FF);
100 plutocolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0x7777FFFF);
101 cometcolor = allocimage(display, Rect(0, 0, 1, 1), RGBA32, 1, 0xAAAAFFFF);
102 font = openfont(display, fontname);
104 fprint(2, "warning: no font %s: %r\n", fontname);
109 * Function heavens() for setting up observer-eye-view
110 * sky charts, plus subsidiary functions.
111 * Written by Doug McIlroy.
114 /* heavens(nlato, wlono, nlatc, wlonc) sets up a map(3)-style
115 coordinate change (by calling orient()) and returns a
116 projection function heavens for calculating a star map
117 centered on nlatc,wlonc and rotated so it will appear
118 upright as seen by a ground observer at nlato,wlono.
119 all coordinates (north latitude and west longitude)
124 mkcoord(double degrees)
128 c.l = degrees*PI/180;
140 v.x = p.nlat.c * p.wlon.c;
141 v.y = -p.nlat.c * p.wlon.s;
148 return a.x*b.x + a.y*b.y + a.z*b.z;
156 v.x = a.y*b.z - a.z*b.y;
157 v.y = a.z*b.x - a.x*b.z;
158 v.z = a.x*b.y - a.y*b.x;
165 return sqrt(dot(a, a));
168 /* An azimuth vector from a to b is a tangent to
169 the sphere at a pointing along the (shorter)
170 great-circle direction to b. It lies in the
171 plane of the vectors a and b, is perpendicular
172 to a, and has a positive compoent along b,
173 provided a and b span a 2-space. Otherwise
174 it is null. (aXb)Xa, where X denotes cross
175 product, is such a vector. */
180 return cross(cross(a,b), a);
183 /* Find the azimuth of point b as seen
184 from point a, given that a is distinct
185 from either pole and from b */
188 azimuth(Xyz a, Xyz b)
190 Xyz aton = azvec(a,north);
191 Xyz atob = azvec(a,b);
192 double lenaton = len(aton);
193 double lenatob = len(atob);
194 double az = acos(dot(aton,atob)/(lenaton*lenatob));
196 if(dot(b,cross(north,a)) < 0)
201 /* Find the rotation (3rd argument of orient() in the
202 map projection package) for the map described
203 below. The return is radians; it must be converted
204 to degrees for orient().
208 rot(struct place center, struct place zenith)
210 Xyz cen = ptov(center);
211 Xyz zen = ptov(zenith);
213 if(cen.z > 1-FUZZ) /* center at N pole */
214 return PI + zenith.wlon.l;
215 if(cen.z < FUZZ-1) /* at S pole */
216 return -zenith.wlon.l;
217 if(fabs(dot(cen,zen)) > 1-FUZZ) /* at zenith */
219 return azimuth(cen, zen);
223 heavens(double zlatdeg, double zlondeg, double clatdeg, double clondeg))(struct place*, double*, double*)
228 center.nlat = mkcoord(clatdeg);
229 center.wlon = mkcoord(clondeg);
230 zenith.nlat = mkcoord(zlatdeg);
231 zenith.wlon = mkcoord(zlondeg);
232 projection = stereographic();
233 orient(clatdeg, clondeg, rot(center, zenith)*180/PI);
238 maptoxy(long ra, long dec, double *x, double *y)
246 pl.nlat.s = sin(lat);
247 pl.nlat.c = cos(lat);
249 pl.wlon.s = sin(lon);
250 pl.wlon.c = cos(lon);
252 return projection(&pl, x, y);
255 /* end of 'heavens' section */
258 setmap(long ramin, long ramax, long decmin, long decmax, Rectangle r, int zenithup)
261 double minx, maxx, miny, maxy;
264 mapra = ramax/2+ramin/2;
265 mapdec = decmax/2+decmin/2;
268 projection = heavens(mylat, mysid, (double)mapdec/c, (double)mapra/c);
270 orient((double)mapdec/c, (double)mapra/c, 0.);
271 projection = stereographic();
273 mapx0 = (r.max.x+r.min.x)/2;
274 mapy0 = (r.max.y+r.min.y)/2;
275 maps = ((double)Dy(r))/(double)(decmax-decmin);
276 if(maptoxy(ramin, decmin, &minx, &miny) < 0)
278 if(maptoxy(ramax, decmax, &maxx, &maxy) < 0)
281 * It's tricky to get the scale right. This noble attempt scales
282 * to fit the lengths of the sides of the rectangle.
284 mapscale = Dx(r)/(minx-maxx);
285 if(mapscale > Dy(r)/(maxy-miny))
286 mapscale = Dy(r)/(maxy-miny);
288 * near the pole It's not a rectangle, though, so this scales
289 * to fit the corners of the trapezoid (a triangle at the pole).
291 mapscale *= (cos(angle(mapdec))+1.)/2;
293 /* upside down, between zenith and pole */
294 fprint(2, "reverse plot\n");
295 mapscale = -mapscale;
301 map(long ra, long dec)
306 if(maptoxy(ra, dec, &x, &y) > 0){
307 p.x = mapscale*x + mapx0;
308 p.y = mapscale*-y + mapy0;
317 dsize(int mag) /* mag is 10*magnitude; return disc size */
321 mag += 25; /* make mags all positive; sirius is -1.6m */
323 /* if plate scale is huge, adjust */
324 if(maps < 100.0/MILLIARCSEC)
326 if(maps < 50.0/MILLIARCSEC)
332 drawname(Image *scr, Image *col, char *s, int ra, int dec)
338 p = addpt(map(ra, dec), Pt(4, -1)); /* font has huge ascent */
339 string(scr, p, col, ZP, font, s);
347 diam = DEG(angle(diam)*MILLIARCSEC); /* diam in milliarcsec */
348 diam /= 2; /* radius */
349 /* convert to plate scale */
350 /* BUG: need +100 because we sometimes crash in library if we plot the exact center */
351 p0 = map(mapra+100, mapdec);
352 p1 = map(mapra+100+diam, mapdec);
353 return hypot(p0.x-p1.x, p0.y-p1.y);
357 drawdisc(Image *scr, DAngle semidiam, int ring, Image *color, Point pt, char *s)
361 d = npixels(semidiam*2);
364 fillellipse(scr, pt, d, d, color, ZP);
366 ellipse(scr, pt, 5*d/2, d/2, 0, color, ZP);
368 ellipse(scr, pt, d, d, 0, grey, ZP);
370 d = stringwidth(font, s);
372 pt.y -= font->height/2 + 1;
373 string(scr, pt, display->black, ZP, font, s);
378 drawplanet(Image *scr, Planetrec *p, Point pt)
380 if(strcmp(p->name, "sun") == 0){
381 drawdisc(scr, p->semidiam, 0, suncolor, pt, nil);
384 if(strcmp(p->name, "moon") == 0){
385 drawdisc(scr, p->semidiam, 0, mooncolor, pt, nil);
388 if(strcmp(p->name, "shadow") == 0){
389 drawdisc(scr, p->semidiam, 2, shadowcolor, pt, nil);
392 if(strcmp(p->name, "mercury") == 0){
393 drawdisc(scr, p->semidiam, 0, mercurycolor, pt, "m");
396 if(strcmp(p->name, "venus") == 0){
397 drawdisc(scr, p->semidiam, 0, venuscolor, pt, "v");
400 if(strcmp(p->name, "mars") == 0){
401 drawdisc(scr, p->semidiam, 0, marscolor, pt, "M");
404 if(strcmp(p->name, "jupiter") == 0){
405 drawdisc(scr, p->semidiam, 0, jupitercolor, pt, "J");
408 if(strcmp(p->name, "saturn") == 0){
409 drawdisc(scr, p->semidiam, 1, saturncolor, pt, "S");
413 if(strcmp(p->name, "uranus") == 0){
414 drawdisc(scr, p->semidiam, 0, uranuscolor, pt, "U");
418 if(strcmp(p->name, "neptune") == 0){
419 drawdisc(scr, p->semidiam, 0, neptunecolor, pt, "N");
423 if(strcmp(p->name, "pluto") == 0){
424 drawdisc(scr, p->semidiam, 0, plutocolor, pt, "P");
428 if(strcmp(p->name, "comet") == 0){
429 drawdisc(scr, p->semidiam, 0, cometcolor, pt, "C");
433 pt.x -= stringwidth(font, p->name)/2;
434 pt.y -= font->height/2;
435 string(scr, pt, grey, ZP, font, p->name);
444 /* stop early to simplify inner loop adjustment */
446 for(i=0,r=rec; i<nrec-nlast; i++,r++)
447 if(r->type == Planet)
448 if(name==nil || strcmp(r->u.planet.name, name)==0){
450 memmove(rec+i, rec+i+1, (nrec-i-1)*sizeof(Record));
459 bbox(long extrara, long extradec, int quantize)
464 int rah, ram, d1, d2;
470 decmax = -0x7FFFFFFF;
472 for(i=0,r=rec; i<nrec; i++,r++){
473 if(r->type == Patch){
474 radec(r->index, &rah, &ram, &dec);
476 r0 = c/cos(dec*PI/180);
485 }else if(r->type==SAO || r->type==NGC || r->type==Planet || r->type==Abell){
488 d1 = 0, d2 = 0, r0 = 0;
491 if(dec+d2+extradec > decmax)
492 decmax = dec+d2+extradec;
493 if(dec-d1-extradec < decmin)
494 decmin = dec-d1-extradec;
500 if(ra+r0+extrara > ramax)
501 ramax = ra+r0+extrara;
502 if(ra-extrara < ramin)
516 /* quantize to degree boundaries */
519 ramax += m5-(ramax%m5);
523 decmin -= c - (-decmin)%c;
526 decmax += c-(decmax%c);
527 }else if(decmax < 0){
529 decmax += ((-decmax)%c);
534 if(ramax-ramin > 270*c){
535 fprint(2, "ra range too wide %ldĀ°\n", (ramax-ramin)/c);
538 }else if(ramax-ramin > 270*c){
540 return bbox(0, 0, quantize);
547 inbbox(DAngle ra, DAngle dec)
556 if(ramin > ramax){ /* straddling 0h0m */
571 mapdec = abs(mapdec)+c/2;
582 #define GREY (nogrey? display->white : grey)
595 int dx, dy, nogrid, textlevel, nogrey, zenithup;
609 if(t = alpha(flags, "nogrid")){
614 if((t = alpha(flags, "zenith")) || (t = alpha(flags, "zenithup")) ){
619 if((t = alpha(flags, "notext")) || (t = alpha(flags, "nolabel")) ){
624 if((t = alpha(flags, "alltext")) || (t = alpha(flags, "alllabel")) ){
629 if(t = alpha(flags, "dx")){
630 dx = strtol(t, &t, 0);
632 fprint(2, "dx %d too small (min 100) in plot\n", dx);
638 if(t = alpha(flags, "dy")){
639 dy = strtol(t, &t, 0);
641 fprint(2, "dy %d too small (min 100) in plot\n", dy);
647 if((t = alpha(flags, "nogrey")) || (t = alpha(flags, "nogray"))){
653 fprint(2, "syntax error in plot\n");
661 if(bbox(0, 0, 1) < 0)
663 if(ramax-ramin<100 || decmax-decmin<100){
664 fprint(2, "plot too small\n");
668 scr = allocimage(display, Rect(0, 0, dx, dy), RGB24, 0, DBlack);
670 fprint(2, "can't allocate image: %r\n");
675 rect = insetrect(rect, 40);
676 if(setmap(ramin, ramax, decmin, decmax, rect, zenithup) < 0){
677 fprint(2, "can't set up map coordinates\n");
682 for(j=0; j<nelem(pts); j++){
683 /* use double to avoid overflow */
684 v = (double)j / (double)(nelem(pts)-1);
685 v = decmin + v*(decmax-decmin);
688 bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
696 p.x -= stringwidth(font, hm5(angle(ra)))/2;
697 string(scr, p, GREY, ZP, font, hm5(angle(ra)));
698 p = pts[nelem(pts)-1];
699 p.x -= stringwidth(font, hm5(angle(ra)))/2;
701 string(scr, p, GREY, ZP, font, hm5(angle(ra)));
708 for(y=decmin; y<=decmax; y+=c){
709 for(j=0; j<nelem(pts); j++){
710 /* use double to avoid overflow */
711 v = (double)j / (double)(nelem(pts)-1);
712 v = ramin + v*(ramax-ramin);
715 bezspline(scr, pts, nelem(pts), Endsquare, Endsquare, 0, GREY, ZP);
718 p.y -= font->height/2;
719 string(scr, p, GREY, ZP, font, deg(angle(y)));
720 p = pts[nelem(pts)-1];
721 p.x -= 3+stringwidth(font, deg(angle(y)));
722 p.y -= font->height/2;
723 string(scr, p, GREY, ZP, font, deg(angle(y)));
726 /* reorder to get planets in front of stars */
728 tolast("moon"); /* moon is in front of everything... */
729 tolast("shadow"); /* ... except the shadow */
731 for(i=0,r=rec; i<nrec; i++,r++){
741 if(name==nil && textlevel>1 && r->type==SAO){
742 snprint(buf, sizeof buf, "SAO%ld", r->index);
746 drawname(scr, nogrey? display->white : alphagrey, name, ra, dec);
748 if(r->type == Planet){
749 drawplanet(scr, &r->u.planet, map(ra, dec));
760 fillellipse(scr, map(ra, dec), m, m, nogrey? display->white : lightgrey, ZP);
762 ellipse(scr, map(ra, dec), m+1, m+1, 0, display->black, ZP);
763 fillellipse(scr, map(ra, dec), m, m, display->white, ZP);
767 if(r->type == Abell){
768 ellipse(scr, addpt(map(ra, dec), Pt(-3, 2)), 2, 1, 0, lightblue, ZP);
769 ellipse(scr, addpt(map(ra, dec), Pt(3, 2)), 2, 1, 0, lightblue, ZP);
770 ellipse(scr, addpt(map(ra, dec), Pt(0, -2)), 1, 2, 0, lightblue, ZP);
773 switch(r->u.ngc.type){
775 j = npixels(r->u.ngc.diam);
782 ellipse(scr, map(ra, dec), j, k, 0, lightblue, ZP);
787 j = npixels(r->u.ngc.diam);
790 ellipse(scr, p, j, j, 0, green, ZP);
791 line(scr, Pt(p.x, p.y+j+1), Pt(p.x, p.y+j+4),
792 Endsquare, Endsquare, 0, green, ZP);
793 line(scr, Pt(p.x, p.y-(j+1)), Pt(p.x, p.y-(j+4)),
794 Endsquare, Endsquare, 0, green, ZP);
795 line(scr, Pt(p.x+j+1, p.y), Pt(p.x+j+4, p.y),
796 Endsquare, Endsquare, 0, green, ZP);
797 line(scr, Pt(p.x-(j+1), p.y), Pt(p.x-(j+4), p.y),
798 Endsquare, Endsquare, 0, green, ZP);
804 j = npixels(r->u.ngc.diam);
807 r1.min = Pt(p.x-j, p.y-j);
808 r1.max = Pt(p.x+j, p.y+j);
809 if(r->u.ngc.type != DiffuseN)
810 draw(scr, r1, ocstipple, ocstipple, ZP);
811 line(scr, Pt(p.x-j, p.y-j), Pt(p.x+j, p.y-j),
812 Endsquare, Endsquare, 0, green, ZP);
813 line(scr, Pt(p.x-j, p.y+j), Pt(p.x+j, p.y+j),
814 Endsquare, Endsquare, 0, green, ZP);
815 line(scr, Pt(p.x-j, p.y-j), Pt(p.x-j, p.y+j),
816 Endsquare, Endsquare, 0, green, ZP);
817 line(scr, Pt(p.x+j, p.y-j), Pt(p.x+j, p.y+j),
818 Endsquare, Endsquare, 0, green, ZP);
823 j = npixels(r->u.ngc.diam);
826 fillellipse(scr, p, j, j, ocstipple, ZP);
830 j = npixels(r->u.ngc.diam);
834 ellipse(scr, p, j, j, 0, lightgrey, ZP);
835 line(scr, Pt(p.x-(j-1), p.y), Pt(p.x+j, p.y),
836 Endsquare, Endsquare, 0, lightgrey, ZP);
837 line(scr, Pt(p.x, p.y-(j-1)), Pt(p.x, p.y+j),
838 Endsquare, Endsquare, 0, lightgrey, ZP);
843 flushimage(display, 1);
848 runcommand(char *command, int p[2])
850 switch(rfork(RFPROC|RFFDG|RFNOWAIT)){
859 execlp("rc", "rc", "-c", command, nil);
860 fprint(2, "can't exec %s: %r", command);
867 parseplanet(char *line, Planetrec *p)
875 line[10] = '\0'; /* terminate name */
876 s = strrchr(line, ' ');
882 for(i=0; s[i]!='\0'; i++)
883 p->name[i] = tolower(s[i]);
885 nfld = getfields(line+11, fld, nelem(fld), 1, " ");
886 p->ra = dangle(getra(fld[0]));
887 p->dec = dangle(getra(fld[1]));
888 p->az = atof(fld[2])*MILLIARCSEC;
889 p->alt = atof(fld[3])*MILLIARCSEC;
890 p->semidiam = atof(fld[4])*1000;
892 p->phase = atof(fld[5]);
898 astro(char *flags, int initial)
902 char cmd[256], buf[4096], *lines[20], *fld[10];
904 snprint(cmd, sizeof cmd, "astro -p %s", flags);
906 fprint(2, "can't pipe: %r\n");
909 if(runcommand(cmd, p) < 0){
912 fprint(2, "can't run astro: %r");
916 n = readn(p[0], buf, sizeof buf-1);
918 fprint(2, "no data from astro\n");
922 Bwrite(&bout, buf, n);
924 np = getfields(buf, lines, nelem(lines), 0, "\n");
926 fprint(2, "astro: not enough output\n");
929 Bprint(&bout, "%s\n", lines[0]);
931 /* get latitude and longitude */
932 if(getfields(lines[0], fld, nelem(fld), 1, " ") < 8)
933 fprint(2, "astro: can't read longitude: too few fields\n");
935 mysid = getra(fld[5])*180./PI;
936 mylat = getra(fld[6])*180./PI;
937 mylon = getra(fld[7])*180./PI;
940 * Each time we run astro, we generate a new planet list
941 * so multiple appearances of a planet may exist as we plot
942 * its motion over time.
944 planet = malloc(NPlanet*sizeof planet[0]);
946 fprint(2, "astro: malloc failed: %r\n");
949 memset(planet, 0, NPlanet*sizeof planet[0]);
951 parseplanet(lines[i], &planet[i-1]);