/* Astronomical Sky Simulator "Skyplot" by David Coffin 1987/11 Wrote similar program in CP/M Basic, text output only, based on program in Oct 87 issue of Sky & Telescope 1988/02/25 Parents bought me a Leading Edge Model D with Hercules Graphics Card 1989/03/24 Wrote Skyplot in Turbo Pascal with graphics 1993/05/14 Rewrote it in Turbo C 1993/11 Bought myself a Tandy 4850 with SVGA, didn't port Skyplot 1995/09 Upgraded from DOS 5.0 to Linux, using XEphem when needed 2018/06/26 Ported to Linux/X11 Changes from the Turbo Pascal/Turbo C code: - Latitude and longitude are decimal instead of DMS. - Two-digit years are prefixed with "20" instead of "19". - Show Mercury in daytime, in addition to Venus and Jupiter. - Allow mouse clicks to quickly identify objects. - Use a faster, standard formula to calculate Keplerian orbits. - All questions may be skipped for default location and current date & time. - Dark side of the moon is drawn with black pixels to simulate eclipses. - Larger star glyphs designed for 3:4 PAR are replaced with solid circles. - Filled in stars south of -50 to magnitude 4.50 like the rest of the sky. - There used to be a text screen to change settings, all poked into video RAM. Could be re-implemented with vt100 codes but not worth the bother. - Added animation for today's much faster computers. */ #define _GNU_SOURCE #include #include #include #include #include #include #include #include #include #include #if !defined(uchar) #define uchar unsigned char #endif #if !defined(ushort) #define ushort unsigned short #endif /* My custom monochrome frame-buffer code. Originally written for the Hercules Graphics Card (720x348, 3:4 pixel aspect ratio) 720x464 achieves the same display aspect ratio with square pixels. */ #define COLOR 0xffbf00 // amber #define WIDE 720 #define HIGH 464 #define PAD 32 #if (WIDE % 8) #error WIDE must be a multiple of 8!! #endif /* Simple monochrome frame buffer: page 0 -- active drawing area page 1 -- background of non-moving objects page 2 -- last image sent to X server */ uchar fbuf[3][HIGH][WIDE/8]; #define INS(x,y) ((unsigned)(x) < WIDE && (unsigned)(y) < HIGH) #define POINTON(page,x,y) fbuf[page][(int)(y)][(int)(x)>>3] |= 1 << ((int)(x) & 7) #define POINTOFF(page,x,y) fbuf[page][(int)(y)][(int)(x)>>3] &= ~(1 << ((int)(x) & 7)) #define TOGGLE(page,x,y) fbuf[page][(int)(y)][(int)(x)>>3] ^= 1 << ((int)(x) & 7) #define GETPOINT(page,x,y) (fbuf[page][(int)(y)][(int)(x)>>3] >> ((int)(x) & 7) & 1) #define SQ(x) ((x)*(x)) #define INSKY(x,y) \ (SQ((x << 1 | 1) - WIDE) + SQ((y << 1 | 1) - (HIGH*2 - WIDE)) < SQ(WIDE+2) && y >= 0) #define SWAP(a,b) { a=a+b; b=a-b; a=a-b; } #define SIZE(alt) 1/SQ(cos(M_PI_4-(alt)/2)) void putpixel (int page, int x, int y, int val) { if (!INS(x,y)) return; if (val) POINTON(page,x,y); else POINTOFF(page,x,y); } void putskypixel (int page, int x, int y, int val) { if (!INSKY(x,y)) return; if (val) POINTON(page,x,y); else POINTOFF(page,x,y); } /* When pixel coords are expressed in floating point, each pixel (x,y) is a unit square tile centered at (x+0.5, y+0.5) and the boundaries of the screen are at 0, 0, WIDE, and HIGH. */ // Avoid the problem that (int) incorrectly maps the interval (-1,0) to 0. #define FLOOR(var,con) ((int)(var + (1+con)) - 1) // Occasionally leaves a one-pixel gap where a line of slope slightly // greater than one meets a line of slope slightly less than one. // I don't think this can be solved without retaining state between calls. void line (float x1, float y1, float x2, float y2) { float f, slope; int i, j, end; if (fabs(x1-x2) > fabs(y1-y2)) { slope = (y1-y2) / (x1-x2); if (x1 > x2) { SWAP (x1,x2); SWAP (y1,y2); } i = FLOOR (x1, 0.5); if (i < 0) i = 0; f = y1 + (i+0.5 - x1)*slope; end = FLOOR (x2, 0.5); if (end > WIDE) end = WIDE; for (; i < end; i++) { if ((unsigned)(j = FLOOR(f,0)) < HIGH) POINTON (0, i, j); f += slope; } } else { slope = (x1-x2) / (y1-y2); if (y1 > y2) { SWAP (x1,x2); SWAP (y1,y2); } i = FLOOR (y1, 0.5); if (i < 0) i = 0; f = x1 + (i+0.5 - y1)*slope; end = FLOOR (y2, 0.5); if (end > HIGH) end = HIGH; for (; i < end; i++) { if ((unsigned)(j = FLOOR(f,0)) < WIDE) POINTON (0, j, i); f += slope; } } } #undef FLOOR double ssqrt (double x) // safe square root { if (x <= 0) return 0; return sqrt(x); } void circle (float x, float y, float rad, int filled) { double rad2 = SQ(rad), clo, chi; int ys, yf, j, i; if ((ys = y - rad) < 0) ys = 0; if ((yf = y + rad + 1) > HIGH) yf = HIGH; for (j = ys; j < yf; j++) { clo = ssqrt (rad2 - SQ(y-j)); chi = ssqrt (rad2 - SQ(y-j-1)); if (clo > chi) SWAP(clo,chi); if (filled) for (i = x-chi; i < x+chi; i++) putskypixel(0,i,j,1); else { i = x-chi; do { putskypixel(0,i,j,1); } while (++i < (int)(x-clo)); i = ceil(x+chi)-1; do { putskypixel(0,i,j,1); } while (i-- > ceil(x+clo)); } } } void rectangle (int xs, int ys, int xf, int yf) { int i; for (i=xs; i <= xf; i++) { putpixel (0, i, ys, 1); putpixel (0, i, yf, 1); } for (i=ys+1; i < yf; i++) { putpixel (0, xs, i, 1); putpixel (0, xf, i, 1); } } void clearrect (int xs, int ys, int xf, int yf) { uchar row[WIDE/8]; int *clip[] = { &xs, &ys, &xf, &yf }; int clim[] = { WIDE-1,HIGH-1,WIDE-1,HIGH-1 }; int xsb, xfb, x, y; if (xs > xf) SWAP(xs,xf); if (ys > yf) SWAP(ys,yf); for (x=0; x < 4; x++) { if (*clip[x] < 0) *clip[x] = 0; if (*clip[x] > clim[x]) *clip[x] = clim[x]; } memset(row,255,sizeof row); xsb = xs >> 3; xfb = xf >> 3; for (x = xsb; x <= xfb; x++) row[x]=0; row[xsb] |= ~(254 << (xs & 7)); row[xfb] |= 254 << (xf & 7); for (y = ys; y <= yf; y++) for (x = xsb; x <= xfb; x++) fbuf[0][y][x] &= row[x]; } void outtextxy (int x, int y, char *st) { /* Taken from the 8x13 X11 font */ static const uchar font[95][11] = { { 0,0,0,0,0,0,0,0,0,0,0 }, { 8,8,8,8,8,8,8,0,8,0,0 }, { 36,36,36,0,0,0,0,0,0,0,0 }, { 0,36,36,126,36,126,36,36,0,0,0 }, { 8,60,10,10,28,40,40,30,8,0,0 }, { 68,74,36,16,16,8,36,84,34,0,0 }, { 0,0,12,18,18,12,82,34,92,0,0 }, { 8,8,8,0,0,0,0,0,0,0,0 }, { 32,16,16,8,8,8,16,16,32,0,0 }, { 4,8,8,16,16,16,8,8,4,0,0 }, { 36,24,126,24,36,0,0,0,0,0,0 }, { 0,0,8,8,62,8,8,0,0,0,0 }, { 0,0,0,0,0,0,0,28,12,2,0 }, { 0,0,0,0,62,0,0,0,0,0,0 }, { 0,0,0,0,0,0,0,8,28,8,0 }, { 64,64,32,16,8,4,2,1,1,0,0 }, { 24,36,66,66,66,66,66,36,24,0,0 }, { 8,12,10,8,8,8,8,8,62,0,0 }, { 60,66,66,64,32,24,4,2,126,0,0 }, { 126,64,32,16,56,64,64,66,60,0,0 }, { 32,48,40,36,34,34,126,32,32,0,0 }, { 126,2,2,58,70,64,64,66,60,0,0 }, { 56,4,2,2,58,70,66,66,60,0,0 }, { 126,64,32,16,16,8,8,4,4,0,0 }, { 60,66,66,66,60,66,66,66,60,0,0 }, { 60,66,66,98,92,64,64,32,28,0,0 }, { 0,0,8,28,8,0,0,8,28,8,0 }, { 0,0,8,28,8,0,0,28,12,2,0 }, { 64,32,16,8,4,8,16,32,64,0,0 }, { 0,0,0,126,0,0,126,0,0,0,0 }, { 2,4,8,16,32,16,8,4,2,0,0 }, { 60,66,66,64,32,16,16,0,16,0,0 }, { 60,66,66,114,74,106,82,2,60,0,0 }, { 24,36,66,66,66,126,66,66,66,0,0 }, { 30,34,66,34,30,34,66,34,30,0,0 }, { 60,66,2,2,2,2,2,66,60,0,0 }, { 30,34,66,66,66,66,66,34,30,0,0 }, { 126,2,2,2,30,2,2,2,126,0,0 }, { 126,2,2,2,30,2,2,2,2,0,0 }, { 60,66,2,2,2,114,66,98,92,0,0 }, { 66,66,66,66,126,66,66,66,66,0,0 }, { 62,8,8,8,8,8,8,8,62,0,0 }, { 120,32,32,32,32,32,32,34,28,0,0 }, { 66,34,18,10,6,10,18,34,66,0,0 }, { 2,2,2,2,2,2,2,2,126,0,0 }, { 65,65,99,85,73,73,65,65,65,0,0 }, { 66,66,70,74,82,98,66,66,66,0,0 }, { 60,66,66,66,66,66,66,66,60,0,0 }, { 62,66,66,66,62,2,2,2,2,0,0 }, { 60,66,66,66,66,66,74,82,60,64,0 }, { 62,66,66,66,62,10,18,34,66,0,0 }, { 60,66,2,2,60,64,64,66,60,0,0 }, { 127,8,8,8,8,8,8,8,8,0,0 }, { 66,66,66,66,66,66,66,66,60,0,0 }, { 65,65,34,34,34,20,20,20,8,0,0 }, { 65,65,65,65,73,73,73,85,34,0,0 }, { 65,65,34,20,8,20,34,65,65,0,0 }, { 65,65,34,20,8,8,8,8,8,0,0 }, { 126,64,32,16,8,4,2,2,126,0,0 }, { 60,4,4,4,4,4,4,4,60,0,0 }, { 1,1,2,4,8,16,32,64,64,0,0 }, { 30,16,16,16,16,16,16,16,30,0,0 }, { 8,20,34,0,0,0,0,0,0,0,0 }, { 0,0,0,0,0,0,0,0,0,127,0 }, { 8,16,0,0,0,0,0,0,0,0,0 }, { 0,0,0,60,64,124,66,98,92,0,0 }, { 2,2,2,58,70,66,66,70,58,0,0 }, { 0,0,0,60,66,2,2,66,60,0,0 }, { 64,64,64,92,98,66,66,98,92,0,0 }, { 0,0,0,60,66,126,2,66,60,0,0 }, { 56,68,4,4,62,4,4,4,4,0,0 }, { 0,0,0,92,34,34,28,2,60,66,60 }, { 2,2,2,58,70,66,66,66,66,0,0 }, { 0,8,0,12,8,8,8,8,62,0,0 }, { 0,32,0,48,32,32,32,32,34,34,28 }, { 2,2,2,34,18,14,18,34,66,0,0 }, { 12,8,8,8,8,8,8,8,62,0,0 }, { 0,0,0,55,73,73,73,73,65,0,0 }, { 0,0,0,58,70,66,66,66,66,0,0 }, { 0,0,0,60,66,66,66,66,60,0,0 }, { 0,0,0,58,70,66,70,58,2,2,2 }, { 0,0,0,92,98,66,98,92,64,64,64 }, { 0,0,0,58,68,4,4,4,4,0,0 }, { 0,0,0,60,66,12,48,66,60,0,0 }, { 0,4,4,62,4,4,4,68,56,0,0 }, { 0,0,0,34,34,34,34,34,92,0,0 }, { 0,0,0,34,34,34,20,20,8,0,0 }, { 0,0,0,65,65,73,73,85,34,0,0 }, { 0,0,0,66,36,24,24,36,66,0,0 }, { 0,0,0,66,66,66,98,92,64,66,60 }, { 0,0,0,126,32,16,8,4,126,0,0 }, { 112,8,8,16,12,16,8,8,112,0,0 }, { 8,8,8,8,8,8,8,8,8,0,0 }, { 14,16,16,8,48,8,16,16,14,0,0 }, { 36,42,18,0,0,0,0,0,0,0,0 }, }; uchar c, i; if (y < 0 || y+11 > HIGH) return; x >>= 3; while ((c = *st++)) { if ((unsigned)(c -= 32) < 95) for (i=0; i < 11; i++) fbuf[0][y+i][x] = font[c][i]; x++; } } struct { char name[20]; int mag; double ra, dec; float alt, azim, xplot, yplot; } object[960], userob[20]; long long seconds; // since 0:00 UTC March 1, 0 AD (aka 1 BC) double sidtime, lat, lon, direc, moondis; int tzone=0, nvis=0, nobj=0, nuser=0, idir, animate=0; Display *dis; int screen, exposed; Window win; XImage *ximage; GC gc; #define BOX WIDE/90 #define DATE (seconds / 86400.0 + 61) // old Skyplot format #define TILT 0.4090509 void set_date (int year, int month, int day, int hour, int minute) { if ((month -= 3) < 0) { month += 12; year--; } seconds = year*365 + year/4 + year/400 - year/100; seconds += (month*306 + 5) / 10 + day - 1; seconds = ((seconds*24 + hour)*60 + minute)*60 - tzone; } void get_date (int *year, int *month, int *day, int *hour, int *minute) { long long i; int cent; i = (seconds + tzone) / 60; *minute = i % 60; i /= 60; *hour = i % 24; i /= 24; cent = (i*4 + 3) / 146097; i += cent - cent/4; *year = (i*4 + 3) / 1461; i = (i - *year * 1461 / 4) * 10 + 5; if ((*month = i / 306 + 3) > 12) { *month -= 12; (*year)++; } *day = i % 306 / 10 + 1; } void rotate (double *x, double *y, double angle) { double cosine, sine, tx; cosine=cos(angle); sine = sin(angle); tx = *x*cosine - *y*sine; *y = *y*cosine + *x*sine; *x = tx; } #define O object[nobj] void loadstars() { FILE *fp; static char name[][8]= { "MOON","SUN","VENUS","JUPITER","MERCURY","MARS", "SATURN","URANUS","NEPTUNE","PLUTO" }; int rahr, ramn, decdg, decmn, i; long long epoch; double precor, ra, dec, vx, vy, vz; char line[64], echar; for (nobj=0; nobj < 10; nobj++) /* Name the sun and planets */ strcpy(O.name,name[nobj]); if (!(fp = fopen("starcat.txt","r"))) { perror("starcat.txt"); exit(1); } printf("Loading star catalog and converting to %6.1lf coordinates...\n",DATE/365.2425); for (; nobj < 922; nobj++) { /* Read the stars */ fgets (line, sizeof line, fp); memcpy (O.name, line, 16); O.name[16] = echar = 0; i = sscanf (line+16,"%d%d%d%d%d %c", &rahr, &ramn, &decdg, &decmn, &O.mag, &echar); if (decdg < 0) decmn = -decmn; epoch = 61530883200 + 1577880000 * (echar == 'G'); /* G = J2000, else J1950 */ precor = (seconds - epoch) * (5028.796195 * M_PI / (3600.0 * 180 * 36525 * 86400)); ra = (rahr*600 + ramn) * (M_PI/7200); dec = (decdg*60 + decmn) * (M_PI/10800); vx = 1 + (vy = vz = 0); rotate (&vx, &vz, dec); rotate (&vx, &vy, ra ); rotate (&vy, &vz, -TILT); rotate (&vx, &vy, precor); rotate (&vy, &vz, TILT); O.dec = asin(vz); O.ra = atan2(-vy,-vx) + M_PI; } fclose(fp); for (i=0; i < nuser; i++) /* Add user-defined objects */ object[nobj++]=userob[i]; } #undef O void prepscreen() { char mon[][10]= { "JANUARY","FEBRUARY","MARCH","APRIL","MAY","JUNE", "JULY","AUGUST","SEPTEMBER","OCTOBER","NOVEMBER","DECEMBER" }; char wkday[][10]= { "WEDNESDAY","THURSDAY","FRIDAY","SATURDAY","SUNDAY","MONDAY","TUESDAY" }; int sidhr, i; float angle, dist, sidmn, x, y, xp, yp; char str[50]; float v0[] = /* Compass letter N */ { -352,3.15452,13,0,15.9,2.53087,13,0,0,0 }; float v1[] = /* Letter E */ { -352,1.55787,9.1,M_PI,13,-M_PI_2,9.1,0,-11.26,2.5213,7,0,0 }; float v2[] = /* Letter S */ { -352,0.01925,9.1,-M_PI_2,6.5,M_PI,9.1,M_PI_2,6.5,M_PI,9.1,-M_PI_2,0 }; float v3[] = /* Letter W */ { -339,-1.5486,13.57,4.4209,9.05,1.8623,9.05,4.4209,13.57,1.8623,0 }; float *vect[4]; int year, month, day, hour, minute, dow; double d; exposed = 1; direc = idir * (M_PI/180); memset (fbuf[0], 0, sizeof fbuf[0]); get_date (&year, &month, &day, &hour, &minute); dow = (seconds + tzone) / 86400 % 7; vect[0]=v0,vect[1]=v1,vect[2]=v2,vect[3]=v3; sprintf (str,"%2d:%02d:%02lld UTC%+d:%02d", hour, minute, seconds%60, tzone/3600, abs(tzone)%3600/60); outtextxy (0,HIGH-40,str); sprintf (str,"%0.4f %0.4f %d",lat*180/M_PI, lon*180/M_PI, idir); outtextxy (0,HIGH-26,str); sprintf (str,"%d %s %d, %s", year, mon[month-1], day, wkday[dow]); outtextxy (0,HIGH-12,str); sprintf (str,"SUN %1.1lf",object[1].alt*(180/M_PI)); outtextxy (WIDE-72,HIGH-54,str); sprintf (str,"MOON %.1f DAYS", /* Approximate phase of moon */ (float) fmod (DATE + 6.05964122, 29.530588)); outtextxy (WIDE-112,HIGH-40,str); sidmn = modf (sidtime*(12/M_PI)+48,&d)*60; sidhr = (int) d % 24; sprintf (str,"SIDEREAL TIME %dh%4.1lfm",sidhr,sidmn); outtextxy (WIDE-176,HIGH-26,str); sprintf (str,"JULIAN DAY %1.3lf",DATE+1721058.5); outtextxy (WIDE-176,HIGH-12,str); circle (WIDE/2, HIGH-WIDE/2, WIDE/2, 0); for (i=0; i<4; i++) /* Draw compass headings */ { xp = WIDE/2; yp = HIGH - xp; while (*vect[i]!=0) { dist = *vect[i]++ * (WIDE/720.0); angle = direc + *vect[i]++; x = xp + fabs(dist) * sin(angle); y = yp - fabs(dist) * cos(angle); if (dist > 0) line (x, y, xp, yp); xp = x; yp = y; } } } #define S object[1] void drawsun() { float rad, rad2; int i; if (S.alt < 0) return; if (S.yplot < -20) return; rad = (WIDE/180.0) * SIZE(S.alt); rad2 = rad * 2.67; for (i=0; i < 8; i++) line (S.xplot, S.yplot, S.xplot + rad2*sin(i*M_PI_4), S.yplot + rad2*cos(i*M_PI_4)); circle (S.xplot, S.yplot, rad, 1); } #define M object[0] void drawmoon() { double sx, sy, sz, dir, rad, rad2, xp, yp; int x0, x1, y0, y1, x, y; if (M.alt < -0.1) return; if (M.yplot < -20) return; sx=0; sy=0; sz=1; /* Which way is the sun? */ rotate (&sz, &sy, S.alt); rotate (&sz, &sx, S.azim-M.azim); rotate (&sy, &sz, M.alt); dir = atan2(-sy,sx) - M.azim + direc; /* 0 = right, M_PI_2 = up */ rad = (WIDE/3.0) / moondis * SIZE(M.alt); rad2 = rad*rad; x0 = M.xplot - rad - 1; x1 = M.xplot + rad + 2; y0 = M.yplot - rad - 1; y1 = M.yplot + rad + 2; for (y = y0; y < y1; y++) for (x = x0; x < x1; x++) { xp = x - M.xplot; yp = y - M.yplot; rotate (&xp, &yp, -dir); if (xp*xp+yp*yp < rad2) putskypixel (0, x, y, xp > sz*sqrt(rad2-yp*yp)); } circle (WIDE/2, HIGH-WIDE/2, WIDE/2, 0); // may have damaged this } #undef M #undef S #define O object[ob] void plotstars() { int ob, x, y; float r; for (ob = 0; ob < nvis; ob++) { r = (WIDE/2-1) * tan (M_PI_4 - O.alt/2); // Stereographic projection O.xplot = WIDE/2 + r*sin(O.azim-direc); O.yplot = HIGH-WIDE/2 + r*cos(O.azim-direc); } drawsun(); drawmoon(); for (ob = 2; ob < nvis; ob++) { if (O.alt < 0) continue; if (O.yplot < -10.0) continue; if (O.mag == 2000) { x = O.xplot; y = O.yplot; line (x-2, y-2, x+3, y+3); line (x-2, y+3, x+3, y-2); } else if (nvis < nobj || O.mag > 350) putpixel (0, O.xplot, O.yplot, 1); else if (O.mag > 250) { putpixel (0, O.xplot-0.5, O.yplot, 1); putpixel (0, O.xplot+0.5, O.yplot, 1); } else if (O.mag > 150) { putpixel (0, O.xplot-0.5, O.yplot-0.5, 1); putpixel (0, O.xplot+0.5, O.yplot-0.5, 1); putpixel (0, O.xplot-0.5, O.yplot+0.5, 1); putpixel (0, O.xplot+0.5, O.yplot+0.5, 1); } else circle (O.xplot, O.yplot, (150-O.mag)/400.0 + 1, 1); } memcpy (fbuf[1],fbuf[0],sizeof *fbuf); } #undef O #define U userob[nuser] void userobject() { int rahr, decdg, decmn; float ramn; char line[64]; do { printf ("\nWhat is the object\'s name? "); fgets (U.name, 16, stdin); printf ("What is its Right Ascension? "); fgets (line, sizeof line, stdin); sscanf (line, "%d%f", &rahr, &ramn); printf ("What is its Declination? "); fgets (line, sizeof line, stdin); sscanf (line, "%d%d", &decdg, &decmn); if (decdg < 0) decmn = -abs(decmn); U.ra = (rahr*15 + ramn*0.25) * (M_PI/180); U.dec = (decdg + decmn/60.0) * (M_PI/180); U.mag = 2000; nuser++; printf ("Another object? "); fgets (line, sizeof line, stdin); } while (toupper(line[0]) == 'Y'); } #undef U /* Calculate the moon's position */ /* Sky & Telescope, July 1989, p. 78 */ #define O object[0] void mooncalc() { static double cor[] = { 0.374897, 0.259091, 0.827362, 0.347343, 0.993126, 0.606434 }, per[] = { 27.55454988, 27.21222075, 29.53058886, -6798.37799, 365.25964, 27.321582268 }; static int amp[]= { 19779,19779,8200,3257,1092,666,-664,-331,-304,-240,226,-108,-79,78, -10828,-1880,-1479,181,-147,-105,-75, 10478,-4105,-2130,-1779,1774,987,-338,-309,-190,-144,-144,-113,-92,-94 }; static char coeff[][5]= { { 0,1,0,1,0 },{ 0,1,0,1,0 },{ 0,1,0,0,0 },{ 1,-1,0,-1,0 }, { 1,1,0,1,0 },{ 1,-1,0,0,0 },{ 1,1,-2,1,0 },{ 0,1,-2,1,0 }, { 0,1,-2,0,0 },{ 1,-1,-2,-1,0 },{ 1,1,0,0,0 },{ 1,1,-2,0,0 }, { 0,1,0,-1,0 },{ 0,1,2,1,0 }, { 1,0,0,0,0 },{ 1,0,-2,0,0 },{ 0,0,2,0,0 },{ 2,0,-2,0,0 }, { 2,0,0,0,0 },{ 0,0,2,0,-1 },{ 1,0,-2,0,1 }, { 1,0,0,0,0 },{ 0,2,0,2,0 },{ 1,0,-2,0,0 },{ 0,2,0,1,0 }, { 0,0,0,1,0 },{ 0,0,2,0,0 },{ 1,-2,0,-2,0 },{ 0,0,0,0,1 }, { 0,2,0,0,0 },{ 1,0,0,1,0 },{ 0,-2,0,-1,0 },{ 1,2,0,2,0 }, { 1,0,-2,0,1 },{ 2,0,-2,0,0 } }; double phase[6], v,u,w,x; int i,j; for (i=0; i<6; i++) phase[i]=fmod((DATE-730486.5)/per[i]+cor[i],1)*(2*M_PI); for (v=0,i=0; i<14; i++) { for (x=0,j=0; j<5; j++) x+=coeff[i][j]*phase[j]; v+=amp[i]*sin(x); } v/=100000; for (u=100000; i<21; i++) { for (x=0,j=0; j<5; j++) x+=coeff[i][j]*phase[j]; u+=amp[i]*cos(x); } u/=100000; for (w=0; i<35; i++) { for (x=0,j=0; j<5; j++) x+=coeff[i][j]*phase[j]; w+=amp[i]*sin(x); } w/=100000; O.ra=asin(w/sqrt(u-v*v))+phase[5]; if (O.ra<0) O.ra+=(2*M_PI); if (O.ra>(2*M_PI)) O.ra-=(2*M_PI); O.dec=asin(v/sqrt(u)); moondis=60.40974*sqrt(u); /* Moon's distance in Earth radii */ O.mag=2000; } #undef O double bright (double th) { /* Calculate brightness change due to phase */ double st,s2t; /* Returns value between 0 and 1 */ th=fabs(th); while (th>M_PI) th-=(2*M_PI); th=fabs(th); st=sin(th); s2t=sin(2*th); return SQ((cos(th)*(M_PI-th+0.5*s2t)+st*st*st-0.29*s2t)/M_PI); } /* Taken from http://www.met.rdg.ac.uk/~ross/Astronomy/Planets.html and various Wikipedia pages. */ #define O object[p+1] void planetcalc () { static const double element[9][6] = { // J2000 epoch { 1.00000011, 0.01671022, 0.00005,-11.26064,102.94719,100.46435 }, // Earth { 0.72333199, 0.00677323, 3.39471, 76.68069,131.53298,181.97973 }, // Venus { 5.20336301, 0.04839266, 1.30530,100.55615, 14.75385, 34.40438 }, // Jupiter { 0.38709893, 0.20563069, 7.00487, 48.33167, 77.45645,252.25084 }, // Mercury { 1.52366231, 0.09341233, 1.85061, 49.57854,336.04084,355.45332 }, // Mars { 9.53707032, 0.05415060, 2.48446,113.71504, 92.43194, 49.94432 }, // Saturn { 19.19126393, 0.04716771, 0.76986, 74.22988,170.96424,313.23218 }, // Uranus { 30.06896348, 0.00858587, 1.76917,131.72169, 44.97135,304.88003 }, // Neptune { 39.48168677, 0.24880766,17.14175,110.30347,224.06676,238.92881 }, // Pluto }; static const double rate[9][6] = { { -0.00000005,-0.00003804,-46.94,-18228.25, 1198.28,129597740.63 }, { 0.00000092,-0.00004938, -2.86, -996.89, -108.80,210664136.06 }, { 0.00060737,-0.00012880, -4.15, 1217.17, 839.93, 10925078.35 }, { 0.00000066, 0.00002527,-23.51, -446.30, 573.57,538101628.29 }, { -0.00007221, 0.00011902,-25.47, -1020.19, 1560.78, 68905103.78 }, { -0.00301530,-0.00036762, 6.11, -1591.05,-1948.89, 4401052.95 }, { 0.00152025,-0.00019150, -2.09, -1681.40, 1312.56, 1542547.79 }, { -0.00125196, 0.0000251, -3.64, -151.25, -844.43, 786449.21 }, { -0.00076912, 0.00006465, 11.07, -37.33, -132.25, 522747.90 }, }; static const int amag[] = { -351,-438,-940,-25,-150,-885,-696,-668,-93 }; double elem[6], pos[3], epos[3], century, err; double a_mean, a_ecc, a_true, sdis, slon, edis, elon; int p, i; // Algol in eclipse? object[72].mag = (fmod((DATE-0.029212) * 0.348759245, 1) < 0.029) ? 340 : 215; century = (seconds - 63108763200) / (36525 * 86400.0); for (p=0; p < 9; p++) { for (i=0; i < 2; i++) elem[i] = element[p][i] + century*rate[p][i]; for (i=2; i < 6; i++) elem[i] = (element[p][i] + century*rate[p][i]/3600) * (M_PI/180); a_mean = fmod (elem[5]-elem[4], 2*M_PI); if (a_mean < 0) a_mean += 2*M_PI; a_ecc = a_mean; do { err = a_mean - a_ecc + elem[1]*sin(a_ecc); a_ecc += err / (1 - elem[1]*cos(a_ecc)); } while (fabs(err) > 1e-12); a_true = 2*atan(sqrt((1+elem[1])/(1-elem[1]) * SQ(tan(a_ecc/2)))); if (a_mean > M_PI) a_true = -a_true; pos[0] = sdis = elem[0] * (1 - elem[1]*cos(a_ecc)); pos[1] = pos[2] = 0; rotate (pos+0, pos+1, a_true + elem[4] - elem[3]); rotate (pos+1, pos+2, elem[2]); rotate (pos+0, pos+1, elem[3] + century*(5028.796195*M_PI/(180*3600))); if (p == 0) { for (i=0; i < 3; i++) epos[i] = pos[i] *= -1; O.mag = -2680; } else { slon = atan2 (pos[1],pos[0]); for (edis = i=0; i < 3; i++) { pos[i] += epos[i]; edis += SQ(pos[i]); } edis = sqrt (edis); elon = atan2 (pos[1],pos[0]); O.mag = floor (amag[p] + log (SQ(sdis*edis) / bright(slon-elon)) * 108.57362 + 0.5); } rotate (pos+1, pos+2, TILT); O.ra = atan2 (-pos[1], -pos[0]) + M_PI; O.dec = atan2 (pos[2], sqrt(SQ(pos[0]) + SQ(pos[1]))); #if 0 int rasec = O.ra/M_PI * (12*36000); int decsec = O.dec/M_PI * (180*3600); printf ("%-10s %2d:%02d:%04.1f %3d:%02d:%02d %9.6f\n", O.name, rasec/36000, rasec/600%60, rasec%600/10.0, decsec/3600, abs(decsec)/60%60, abs(decsec)%60, edis); #endif switch (p) { case 1: // Venus' atmosphere if (O.mag > -380 && edis < 1) O.mag=-380; break; case 5: // Saturn's rings O.mag -= floor(log(1 + fabs(sin(elon))) * 108.57362 + 0.5); } } } #undef O #define O object[b] void compaltaz() { int b; double px, py, pz; sidtime = modf (seconds/86164.090530832886, &px) * (2*M_PI) + 2.76617713242554 + lon; for (b=0, nvis = nobj; b < nvis; b++) { px = py = 0; pz = 1; /* px=South, py=West, pz=Zenith */ rotate (&px, &pz, O.dec); rotate (&pz, &py, sidtime-O.ra); rotate (&pz, &px, lat); O.azim = atan2(py,px) + M_PI; O.alt = asin(pz); if (b==0) O.alt -= atan2 (cos(O.alt), moondis-sin(O.alt)); if (b==1 && O.alt > -0.1) nvis = 5; } } void identify (int b) { char str[50]; double rahr,ramn,decdg,decmn; rectangle (0,0,234,46); clearrect (1,1,233,45); if (b < 0) { outtextxy(40,19,"NO OBJECT SEEN"); return; } if (O.mag == 2000) strcpy (str, O.name); else sprintf (str, "%s MAG %5.2f", O.name, O.mag*0.01); outtextxy(16,5,str); if (O.ra < 0) O.ra += 2*M_PI; ramn = modf (O.ra*(12/M_PI), &rahr)*60; decmn = modf (fabs(O.dec)*(180/M_PI)+0.008333, &decdg)*60; sprintf(str,"R.A. %dh %1.1lfm DEC %c%d %d\'", (int) rahr, ramn, O.dec < 0 ? '-':32, (int) decdg, (int) decmn); outtextxy(16,19,str); sprintf(str,"ALT %1.1lf AZIM %1.1lf",O.alt*(180/M_PI),O.azim*(180/M_PI)); outtextxy(16,33,str); } void display_x() { unsigned *old, *new, i, b, mask, pix, c, last_c=-1; if (exposed) { XSetForeground (dis, gc, COLOR); XPutImage (dis, win, gc, ximage, 0, 0, PAD, PAD, WIDE, HIGH); memcpy (fbuf[2], fbuf[0], sizeof *fbuf); exposed = 0; } else { old = (unsigned *) fbuf[2]; new = (unsigned *) fbuf[0]; for (i=0; i < HIGH*WIDE/32; i++) if ((mask = new[i] ^ old[i])) { for (b=0; b < 32; b++) if (mask & (1 << b)) { pix = i << 5 | b; c = new[i] >> b & 1; if (c != last_c) XSetForeground (dis, gc, COLOR * (last_c = c)); XDrawPoint (dis, win, gc, pix%WIDE+PAD, pix/WIDE+PAD); } old[i] = new[i]; } } } void event_loop() { XEvent xev; char key; static int boxx=(WIDE-BOX)/2, boxy=HIGH-(WIDE+BOX)/2; static int shift=0, control=0, anim=0; static const int anim_speed[] = { 0, 10, 60, 360, 86164, 86400 }; static const signed char dx[] = { -1,0,1,0,-1,0,1,0,-1,0,1 }; static const signed char dy[] = { -1,-1,-1,0,0,0,0,0,1,1,1 }; float dist, best; int b, found; dis = XOpenDisplay(0); screen = DefaultScreen (dis); win = XCreateSimpleWindow (dis, DefaultRootWindow(dis),0,0, WIDE+PAD*2, HIGH+PAD*2, 5, 0, 0); XSetStandardProperties (dis, win, "SKYPLOT","SKYPLOT",None,0,0,0); XSelectInput (dis, win, ExposureMask|KeyPressMask|KeyReleaseMask|ButtonPressMask); gc = XCreateGC (dis, win, 0,0); XClearWindow (dis, win); XMapRaised (dis, win); ximage = XCreateImage (dis, DefaultVisual(dis,screen), 1, XYBitmap, 0, (char *)fbuf[0][0], WIDE, HIGH, 8, WIDE/8); recalc: seconds += animate; if (animate) usleep(100000); mooncalc(); planetcalc(); compaltaz(); redraw: prepscreen(); plotstars(); while (1) { display_x(); if (animate) { if (!XCheckWindowEvent (dis, win, ExposureMask|KeyPressMask|KeyReleaseMask|ButtonPressMask, &xev)) goto recalc; } else XWindowEvent (dis, win, ExposureMask|KeyPressMask|KeyReleaseMask|ButtonPressMask, &xev); if (xev.type == Expose) { exposed = 1; continue; } if (xev.type == ButtonPress) { if (!INSKY (xev.xkey.x-PAD, xev.xkey.y-PAD)) continue; for (best=FLT_MAX, b=0; b < nvis; b++) { dist = SQ(xev.xkey.x-PAD-O.xplot-0.5) + SQ(xev.xkey.y-PAD-O.yplot-0.5); if (best > dist) { best = dist; found = b; } } boxx = object[found].xplot - BOX*0.5; boxy = object[found].yplot - BOX*0.5; memcpy (fbuf[0], fbuf[1], sizeof fbuf[0]); rectangle (boxx, boxy, boxx+BOX, boxy+BOX); identify (found); continue; } key = xev.xkey.keycode; if (key == 50 || key == 62) shift = xev.type == KeyPress; if (key == 37 || key == 105) control = xev.type == KeyPress; if (xev.type == KeyRelease) continue; if (key == 22) { animate = -animate; anim = -anim; continue; } if (key == 65) { if (animate) animate = 0; else goto restart_anim; continue; } switch (key) { case 111: key = 80; break; case 113: key = 83; break; case 114: key = 85; break; case 116: key = 88; break; } if (control) switch (key) { case 83: idir += 358 - 18*!shift; case 85: idir += 361 + 9*!shift; idir %= 360; goto redraw; case 80: if (anim++ < 5) anim++; case 88: if (anim > -5) anim--; restart_anim: animate = anim_speed[abs(anim)]; if (anim < 0) animate = -animate; default: continue; } if ((unsigned)(key-79) < 11) { boxx += dx[key-79] * (shift ? 1:BOX); boxy += dy[key-79] * (shift ? 1:BOX); memcpy (fbuf[0], fbuf[1], sizeof fbuf[0]); rectangle (boxx, boxy, boxx+BOX, boxy+BOX); } if (key == 36 || key == 104) { for (found = -1, b=0; b < nvis; b++) if (O.alt > 0 && (unsigned)O.xplot - boxx <= BOX && (unsigned)O.yplot - boxy <= BOX ) { found = b; break; } identify (found); } if (key == 9) break; if (animate) goto recalc; } XFreeGC (dis, gc); XDestroyWindow (dis, win); XCloseDisplay (dis); } #undef O int main() { char line[64]; int year, month, day, hour, minute; time_t now; struct tm tm; puts("\t\tSKYPLOT: A Computerized Planetarium"); puts("\t\t\tby David Coffin\n"); puts("SKYPLOT can display the sky as it appears from any location on"); puts("Earth for any date or time. The catalog includes all planets"); puts("and 913 stars. Enter the data according to the examples in ().\n"); printf("Do you want to add your own objects? "); fgets (line, sizeof line, stdin); if (toupper(line[0]) == 'Y') userobject(); lat = 42.728763; lon = -72.266711; printf ("\nEnter your latitude and longitude (42.196,-71.304): "); fgets (line, sizeof line, stdin); sscanf (line, "%lf,%lf", &lat, &lon); lat *= M_PI/180; lon *= M_PI/180; idir = 180 * (lat >= 0); now = time(0); localtime_r (&now, &tm); seconds = 62162035200 + now; tzone = tm.tm_gmtoff; get_date (&year, &month, &day, &hour, &minute); printf("\nEnter the date and time (2/25/1988 0:0) (6/20/18 21:35): "); fgets (line, sizeof line, stdin); sscanf (line, "%d/%d/%d %d:%d", &month, &day, &year, &hour, &minute); if (year < 100) year += 2000; printf ("Enter your time zone (GMT=0, EDT=4, EST=5): "); fgets (line, sizeof line, stdin); if (sscanf (line, "%d", &tzone) == 1) tzone *= -3600; set_date (year, month, day, hour, minute); puts ("\nKEYBOARD COMMANDS:"); puts ("<← → ↑ ↓>: Move the box (Shift for more precision)"); puts (": Identify the boxed object"); puts (": Change the view direction"); puts (": Adjust animation speed"); puts (": Reverse animation"); puts (": Pause animation"); puts (": Quit"); loadstars(); event_loop(); }