rck@ihuxx.UUCP (Kukuk) (01/24/86)
This program computes all sorts of data concerning the Sun, Moon, and planets. Refer to the README for more details. #!/bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #!/bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Makefile # README # conv.c # glob.c # io.c # main.c # misc.c # p.h # sm.c # time.c # wand.c # This archive created: Fri Jan 24 13:15:16 1986 export PATH; PATH=/bin:$PATH if test -f 'Makefile' then echo shar: over-writing existing file "'Makefile'" fi cat << \SHAR_EOF > 'Makefile' SRC=main.c time.c conv.c wand.c sm.c misc.c glob.c io.c p: main.o time.o conv.o wand.o sm.o misc.o glob.o io.o Makefile cc main.o time.o conv.o wand.o sm.o misc.o glob.o io.o \ -o p -lm -lcurses main.o time.o conv.o wand.o sm.o misc.o glob.o io.o : p.h print: @prdef -args -c82 p.h $(SRC) cf: @cf -n -B $(SRC) | pr -l88 -h p lint: lint $(SRC) -lm SHAR_EOF if test -f 'README' then echo shar: over-writing existing file "'README'" fi cat << \SHAR_EOF > 'README' "P" is a program which displays all sorts of astronomical data about the Sun, Moon, and planets. The algorithms were taken from Duffet-Smith's book, "Practical Astronomy with your Calculator". All of the accuracy limitations described in the book apply, of course, to this program. (Positions +- ~1 deg.; rise/set times +- ~5 min.) The program displays the following data in the upper part of the screen: - RA (in hours and minutes) of Sun, Moon, and planets - DEC (in degrees) of the same - Azimuth and elevation (in degrees) of same - Rise and set times of same - Distance (in astronomical units) of same - Light travel time (in hours and minutes) of same - Apparent diameter (in arc seconds) of same - Phase of same - Apparent magnitude of same A matrix of angular separation (in degrees) of the planets is displayed in the lower half of the screen. A line near the bottom of the screen displays the date and the time. The time is presented in Local, Greenwich Mean, Greenwich Siderial, and Local Siderial formats. The bottom line is used for input when the program is run in interactive mode. The Lunar and Solar eclipse calculations are included, the comet calculations are not. Here are some examples of running the program: $ p In this mode, "p" will prompt you for the time (HH:MM:SS) and date (MM/DD/YYYY). The program will display the values and exit. Use a negative YYYY for years B. C. $ p -n In this mode, "p" uses the current time and date for the calculations. $ p -n -l The program will use the current time and date for the first calculation and then repeat the calculations using a default delta of one day. The program is running in interactive mode, so you can vary the delta as you watch. The interactive commands are: >>> r<cr> for reverse (or negate) the delta value >>> <number><X><cr> where <number> is a floating point number, <X> is one of NULL for day (default) s for second m for minute h for hour d for day w for week y for year >>> n<cr> reset the internal date to "now" >>> e<cr> enter elapsed time mode; this only works for reasonable values of delta. >>> ne<cr> leave elapsed time mode >>> nl<cr> leave looping mode >>> q<cr> quit the program If you want printed output instead of displayed output, add "-p" to the command line. The default delta of one day in looping mode can be overridden by adding the desired delta after the "-l" (with no white space) on the command line. The format of the delta is that of the interactive delta described above. Necessary local mods: There are four #define's, HERELONG, HERELAT, HERELEV, and HEREZONE, that must be defined to be your longitude, latitude, elevation, and timezone, respectively. Look in p.h for more details. This program has been run only on AT&T versions of UN*X. Enjoy. Ron Kukuk SHAR_EOF if test -f 'conv.c' then echo shar: over-writing existing file "'conv.c'" fi cat << \SHAR_EOF > 'conv.c' #include "p.h" /* * (Sect. 24; Pg. 39) */ REAL ratoha(alpha) REAL alpha; { REAL lsthour; lsthour = hmstohr(&now.lst); lsthour -= alpha; lsthour = tcanon(lsthour); return(lsthour); } /* * (Sect. 24; Pg. 39) */ REAL hatora(H) REAL H; { REAL lsthour; lsthour = hmstohr(&now.lst); lsthour -= H; lsthour = tcanon(lsthour); return(lsthour); } /* * (Sect. 25; Pg. 40) */ struct ae * eqtohzn(eqp) register struct equat *eqp; { static struct ae ae; REAL a, A, delta, H; H = ratoha(eqp->ra); /* its hour angle */ H *= 15.0; /* hour angle to degrees */ delta = eqp->dec; a = (REAL) sin(delta/RAD) * (REAL) sin(HERELAT/RAD) + (REAL)cos(delta/RAD) * (REAL)cos(HERELAT/RAD) * (REAL)cos(H/RAD); a = (REAL) asin(a) * RAD; A = ((REAL) sin(delta/RAD) - (REAL) sin(HERELAT/RAD) * (REAL) sin(a/RAD)) / ((REAL) cos(HERELAT/RAD) * (REAL) cos(a/RAD)); A = (REAL) acos(A) * RAD; if (sin(H/RAD) > 0.0) A = 360.0 - A; ae.az = A; ae.el = a; return(&ae); } /* * (Sect. 26; Pg. 42) */ struct equat * hzntoeq(aep) register struct ae *aep; { static struct equat equat; REAL a, A, delta, H; A = aep->az; a = aep->el; delta = (REAL) sin(a/RAD) * (REAL) sin(HERELAT/RAD) + (REAL) cos(a/RAD) * (REAL) cos(HERELAT/RAD) * (REAL) cos(A/RAD); delta = (REAL) asin(delta) * RAD; H = ((REAL) sin(a/RAD) - (REAL) sin(HERELAT/RAD) * (REAL) sin(delta/RAD)) / ((REAL) cos(HERELAT/RAD) * (REAL) cos(delta/RAD)); H = (REAL) acos(H) * RAD; if (sin(A/RAD) > 0.0) H = 360.0 - H; equat.ra = H / 15.0; /* convert deg to hour angle */ /* * Then convert hour angle to RA. */ equat.ra = hatora(equat.ra); equat.dec = delta; return(&equat); } /* * (Sect. 27; Pg. 44) */ struct equat * ectoeq(ecp) register struct eclipt *ecp; { REAL alpha, beta, delta, lambda, epsilon, x, y; static struct equat equat; beta = ecp->lat / RAD; lambda = ecp->lon / RAD; epsilon = obliq(now.jd) / RAD; delta = (REAL) sin(beta) * (REAL) cos(epsilon) + (REAL) cos(beta) * (REAL) sin(epsilon) * (REAL) sin(lambda); delta = asin(delta) * RAD; y = (REAL) sin(lambda) * (REAL) cos(epsilon) - (REAL) tan(beta) * (REAL) sin(epsilon); x = (REAL) cos(lambda); alpha = atan2a(y, x) * RAD; alpha /= 15.0; /* degrees to RA */ equat.ra = alpha; equat.dec = delta; return(&equat); } /* * (Sect. 28; Pg. 46) */ struct eclipt * eqtoec(eqp) register struct equat *eqp; { REAL alpha, beta, delta, lambda, epsilon, x, y; static struct eclipt eclipt; delta = eqp->dec / RAD; alpha = (eqp->ra * 15.0) / RAD; epsilon = obliq(now.jd) / RAD; beta = (REAL) sin(delta) * (REAL) cos(epsilon) - (REAL) cos(delta) * (REAL) sin(epsilon) * (REAL) sin(alpha); beta = asin(beta) * RAD; y = (REAL) sin(alpha) * (REAL) cos(epsilon) + (REAL) tan(delta) * (REAL) sin(epsilon); x = (REAL) cos(alpha); lambda = atan2a(y, x) * RAD; eclipt.lat = beta; eclipt.lon = lambda; return(&eclipt); } REAL obliq(jd) REAL jd; { REAL T, de, epsilon; jd -= 2415020.0; /* JD for 1900 */ T = jd / 36525.0; de = 46.845 * T + 0.0059 * T * T - 0.00181 * T * T * T; de /= 3600.0; epsilon = 23.452294 - de; return(epsilon); } SHAR_EOF if test -f 'glob.c' then echo shar: over-writing existing file "'glob.c'" fi cat << \SHAR_EOF > 'glob.c' #include "p.h" int dmsize[12] = { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 }; char month[12][3] = { "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec" }; char days[7][3] = { "Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat" }; /* * Everything worth knowing about Sol. The ecliptic field must * be kept up to date by calling sunorbit(); */ struct sun sun; /* * Everything worth knowing about Luna. */ struct moon moon; /* * Planetary equatorial coordinates; kept and used for angular * separation calculations. */ struct equat wandeq[8]; /* * Planetary orbital data. */ struct pod pod[] = { /* Mercury */ { 0.24085, 231.2973, 77.1442128, 0.2056306, 0.3870986, 7.0043579, 48.0941733, 6.74, 1.918e-6 }, /* Venus */ { 0.61521, 355.73352, 131.2895792, 0.0067826, 0.7233316, 3.394435, 76.4997524, 16.92, 1.721e-5 }, /* Mars */ { 1.88089, 126.30783, 335.6908166, 0.0933865, 1.5236883, 1.8498011, 49.4032001, 9.36, 4.539e-6 }, /* Jupiter */ { 11.86224, 146.966365, 14.0095493, 0.0484658, 5.202561, 1.3041819, 100.2520175, 196.74, 1.994e-4 }, /* Saturn */ { 29.45771, 165.322242, 92.6653974, 0.0556155, 9.554747, 2.4893741, 113.4888341, 165.60, 1.740e-4 }, /* Uranus */ { 84.01247, 228.0708551, 172.7363288, 0.0463232, 19.21814, 0.7729895, 73.8768642, 65.80, 7.768e-5 }, /* Neptune */ { 164.79558, 260.3578998, 47.8672148, 0.0090021, 30.10957, 1.7716017, 131.5606494, 62.20, 7.597e-5 }, /* Pluto (osculating elements for 1/2/1980) */ { 250.9, 209.439, 222.972, 0.25387, 39.78459, 17.137, 109.941, 8.20, 4.073e-6 } }; /* * Keep Earth orbital data separate for ease in subscripting the others. */ struct pod epod = { { 1.00004, 98.833540, 102.596403, 0.016718, 1.0, 0.0, 0.0, 0.0, 0.0 } }; /* * CRT screen coordinates for displaying the results. */ struct pt pt[] = { { 1, 8 }, /* Mercury */ { 2, 8 }, /* Venus */ { 3, 8 }, /* Mars */ { 4, 8 }, /* Jupiter */ { 5, 8 }, /* Saturn */ { 6, 8 }, /* Uranus */ { 7, 8 }, /* Neptune */ { 8, 8 }, /* Pluto */ { 10, 8 }, /* Sun */ { 11, 8 }, /* Moon */ { 13, 0 }, /* Separation matrix */ { 22, 0 }, /* Time */ }; char *name[] = { "Mercury", "Venus", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto", "Sun", "Moon", "Angsep", "" }; /* * This is the now structure. Its contents should be kept up to date * by calling tempus(). */ struct now now; /* * Orbital position of Earth (heliocentric coodinates) computed * from the time in the now structure. Kept up to date by calling * earthorbit(). */ struct planpos earth; int iter; /* iteration count for Keplerian calculations */ int flags; /* file control flags */ int print; /* hard copy output flag */ int loop; /* keep looping */ int ritenow; /* compute results for today's date and time */ int request; /* user has request pending */ extern int daylight; /* enable DST conversions */ int eltime; /* sleep away the actual interval */ REAL deltat; /* loop increment in REAL days */ int donin; /* input line is complete; parse it */ SHAR_EOF if test -f 'io.c' then echo shar: over-writing existing file "'io.c'" fi cat << \SHAR_EOF > 'io.c' #include "p.h" /* * Get the requested local time into the local time part of the * now structure. */ gettime() { char buf[40]; int nitems; do { fprintf(stderr, "Time (hh:mm:ss): "); gets(buf); nitems = sscanf(buf, "%2d:%2d:%2d", &now.lt.tm_hour, &now.lt.tm_min, &now.lt.tm_sec); } while (nitems != 3); do { fprintf(stderr, "Date (mm/dd/yyyy): "); gets(buf); nitems = sscanf(buf, "%2d/%2d/%5d", &now.lt.tm_mon, &now.lt.tm_mday, &now.lt.tm_year); } while (nitems != 3); /* * Add day of year. */ dayofyr(&now.lt); /* * Maybe adjust for DST. */ dstime(&now.lt); } struct equat * getequat() { static struct equat eq; struct deg deg; char buf[40]; int nitems; do { printf("Hour angle (hh:mm:ss): "); gets(buf); nitems = sscanf(buf, "%2d:%2d:%2d", °.deg, °.min, °.sec); } while (nitems != 3); eq.ra = dmstodeg(°); do { printf("Decl (deg:min:ss): "); gets(buf); nitems = sscanf(buf, "%2d:%2d:%2d", °.deg, °.min, °.sec); } while (nitems != 3); eq.dec = dmstodeg(°); return(&eq); } /* * Compute and display angular separation data. */ prtsep() { register p1, p2, p3; REAL seps[8][8]; REAL smsep; /* * Compute the angular separations. */ for (p1 = MERCURY; p1 <= PLUTO; p1++) { for (p2 = p1 + 1; p2 <= PLUTO; p2++) { seps[p1][p2] = angsep(&wandeq[p1], &wandeq[p2]); } } smsep = angsep(&sun.equat, &moon.equat); if ( ! print) { /* * Display the angular separations. */ for (p1 = MERCURY; p1 <= PLUTO; p1++) { for (p2 = p1 + 1; p2 <= PLUTO; p2++) { move(pt[SEPMAT].row+p1+1, p2*8); if (seps[p1][p2] <= 5.0) { standout(); printw("%7.3f", seps[p1][p2]); standend(); } else printw("%7.3f", seps[p1][p2]); } } move(pt[SEPMAT].row+6, 8); printw("%7.3f", smsep); } else { /* hard copy */ printf("\n "); /* * Print the labels. */ for (p1 = VENUS; p1 <= PLUTO; p1++) printf("%7.7s ", name[p1]); putchar('\n'); for (p1 = MERCURY; p1 < PLUTO; p1++) { printf(" "); for (p3 = 0; p3 < p1; p3++) printf(" "); for (p2 = p1 + 1; p2 <= PLUTO; p2++) { printf("%7.3f ", seps[p1][p2]); } printf("%s\n", name[p1]); } printf("\tSun-Moon: %7.3f\n", smsep); dashes(); } } prtshdr() { register p1; if (print) return; /* * Toss the label. */ prtname(SEPMAT); /* * Print the labels. (Note the limits.) */ for (p1 = VENUS; p1 <= PLUTO; p1++) { move(pt[SEPMAT].row, (p1-1)*8+8); printw("%7.7s", name[p1]); } for (p1 = MERCURY; p1 < PLUTO; p1++) { /* move(pt[SEPMAT].row+p1+1, 1); printw("."); */ move(pt[SEPMAT].row+p1+1, 64); printw("%s", name[p1]); } move(pt[SEPMAT].row+5, 8); printw(" Sun"); move(pt[SEPMAT].row+6, 16); printw("Moon"); } static char *hdrhdr = " rh rm dec azd eld rise set disAU h mn dias phse magnit"; prthdr() { if (print) { dashes(); printf(hdrhdr); putchar('\n'); putchar('\n'); } else { move(0, 0); printw(hdrhdr); } } dashes() { register i; for (i = 0; i < 76; i++) putchar('-'); putchar('\n'); } prtname(plano) { if (print) { if (plano == MERCURY) prthdr(); if (plano >= SUN) putchar('\n'); printf(name[plano]); if (plano != TIME) putchar('\t'); } else { move(pt[plano].row, 0); printw(name[plano]); } } prtmp(tmp) register struct tm *tmp; { switch (tmp->tm_isdst) { case DST: case LT: case LST: case GMT: case GST: printf("%2dh %2dm %2ds", tmp->tm_hour, tmp->tm_min, tmp->tm_sec); switch (tmp->tm_isdst) { case DST: printf( " DST"); break; case LT: printf( " LT "); break; case LST: printf( " LST"); break; case GMT: printf( " GMT"); break; case GST: printf( " GST"); break; } putchar('\n'); break; default: printf("%2dh %2dm %2ds isdst: %d\n", tmp->tm_hour, tmp->tm_min, tmp->tm_sec, tmp->tm_isdst); } } prtdate(tmp) register struct tm *tmp; { register int yr; register char *day; day = days[tmp->tm_wday]; yr = tmp->tm_year; printf("%3s %2d/%2d/%4d", day, tmp->tm_mon, tmp->tm_mday, yr); printf(" wday: %d yday: %3d\n", tmp->tm_wday, tmp->tm_yday); } prtae(aep) register struct ae *aep; { struct deg az, el; degtodms(&az, aep->az); degtodms(&el, aep->el); printf("Az: %3dd %2dm %2ds ", az.deg, az.min, az.sec); printf("El: %c%2dd %2dm %2ds \n", el.sg, el.deg, el.min, el.sec); } prteq(eqp) register struct equat *eqp; { struct deg ra, dec; degtodms(&ra, eqp->ra); degtodms(&dec, eqp->dec); printf("Ra: %2dh %2dm %2ds ", ra.deg, ra.min, ra.sec); printf("Dec: %c%2dd %2dm %2ds\n", dec.sg, dec.deg, dec.min, dec.sec); } static char *planfmt = "%2d %2d %3d %3d %3d %02d:%02d %02d:%02d %3s %5.2f %d %2d %4.1f %4.2f %6.2f"; prtplan(planp, plano) register struct planet *planp; { struct tm tm, *tmp; struct deg ra, dec, az, el; struct deg rt, st; char *timek(); degtodms(&ra, planp->equat.ra); degtodms(&dec, planp->equat.dec); if (dec.sg == '-') dec.deg = -dec.deg; degtodms(&az, planp->hzn.az); degtodms(&el, planp->hzn.el); if (el.sg == '-') el.deg = -el.deg; hrtohms(&tm, planp->rs.risetime); tmp = lsttolt(&tm); rt.deg = tmp->tm_hour; rt.min = tmp->tm_min; hrtohms(&tm, planp->rs.settime); tmp = lsttolt(&tm); st.deg = tmp->tm_hour; st.min = tmp->tm_min; if (print) { printf(planfmt, ra.deg, ra.min, dec.deg, az.deg, el.deg, rt.deg, rt.min, st.deg, st.min, timek(tmp), planp->rho, planp->ltt.deg, planp->ltt.min, planp->th, planp->F, planp->m); putchar('\n'); } else { move(pt[plano].row, pt[plano].col); printw(planfmt, ra.deg, ra.min, dec.deg, az.deg, el.deg, rt.deg, rt.min, st.deg, st.min, timek(tmp), planp->rho, planp->ltt.deg, planp->ltt.min, planp->th, planp->F, planp->m); } } char * timek(tmp) register struct tm *tmp; { char *p; switch (tmp->tm_isdst) { case LT: p = "LT "; break; case DST: p = "DST"; break; case GMT: p = "GMT"; break; case GST: p = "GST"; break; case LST: p = "LST"; break; default: p = "???"; break; } return(p); } static char *sunfmt = "%2d %2d %3d %3d %3d %02d:%02d %02d:%02d %3s %5.4fAU %2dm %2ds"; prtsun() { struct tm tm, *tmp; struct deg ra, dec, az, el; struct deg rt, st; struct deg adia; char *timek(); degtodms(&ra, sun.equat.ra); degtodms(&dec, sun.equat.dec); if (dec.sg == '-') dec.deg = -dec.deg; degtodms(&az, sun.hzn.az); degtodms(&el, sun.hzn.el); if (el.sg == '-') el.deg = -el.deg; hrtohms(&tm, sun.rs.risetime); tmp = lsttolt(&tm); rt.deg = tmp->tm_hour; rt.min = tmp->tm_min; hrtohms(&tm, sun.rs.settime); tmp = lsttolt(&tm); st.deg = tmp->tm_hour; st.min = tmp->tm_min; degtodms(&adia, sun.th); if (print) { printf(sunfmt, ra.deg, ra.min, dec.deg, az.deg, el.deg, rt.deg, rt.min, st.deg, st.min, timek(tmp), sun.dist, adia.min, adia.sec); putchar('\n'); } else { move(pt[SUN].row, pt[SUN].col); printw(sunfmt, ra.deg, ra.min, dec.deg, az.deg, el.deg, rt.deg, rt.min, st.deg, st.min, timek(tmp), sun.dist, adia.min, adia.sec); } } static char *moonfmt = "%2d %2d %3d %3d %3d %02d:%02d %02d:%02d \ %3s %6.0fkm %2dm %2ds %4.2f %4.1fda"; prtmoon() { struct tm tm, *tmp; struct deg ra, dec, az, el; struct deg rt, st; struct deg adia; char *timek(); degtodms(&ra, moon.equat.ra); degtodms(&dec, moon.equat.dec); if (dec.sg == '-') dec.deg = -dec.deg; degtodms(&az, moon.hzn.az); degtodms(&el, moon.hzn.el); if (el.sg == '-') el.deg = -el.deg; hrtohms(&tm, moon.rs.risetime); tmp = lsttolt(&tm); rt.deg = tmp->tm_hour; rt.min = tmp->tm_min; hrtohms(&tm, moon.rs.settime); tmp = lsttolt(&tm); st.deg = tmp->tm_hour; st.min = tmp->tm_min; degtodms(&adia, moon.th); if (print) { printf(moonfmt, ra.deg, ra.min, dec.deg, az.deg, el.deg, rt.deg, rt.min, st.deg, st.min, timek(tmp), moon.dist, adia.min, adia.sec, moon.F, moon.age); putchar('\n'); } else { move(pt[MOON].row, pt[MOON].col); printw(moonfmt, ra.deg, ra.min, dec.deg, az.deg, el.deg, rt.deg, rt.min, st.deg, st.min, timek(tmp), moon.dist, adia.min, adia.sec, moon.F, moon.age); } } prtimes() { register int yr; register char *day, *mon, *ab; char *dfmt = "%3.3s, %3.3s %2d, %4d %2.2s"; day = days[now.lt.tm_wday]; yr = now.lt.tm_year; mon = month[now.lt.tm_mon-1]; if (yr > 0) ab = "AD"; else { ab = "BC"; yr = -yr; } if (print) printf(dfmt, day, mon, now.lt.tm_mday, yr, ab); else { move(pt[TIME].row, pt[TIME].col); printw(dfmt, day, mon, now.lt.tm_mday, yr, ab); } if ( ! print) move(pt[TIME].row, 22); else putchar(' '); prtim(&now.lt); if ( ! print) move(pt[TIME].row, 36); else printf(" "); prtim(&now.gmt); if ( ! print) move(pt[TIME].row, 50); else printf(" "); prtim(&now.gst); if ( ! print) move(pt[TIME].row, 64); else printf(" "); prtim(&now.lst); /* * For pr page alignment. */ if (print) putchar('\n'); } prtim(tmp) register struct tm *tmp; { char *tfmt = "%02d:%02d:%02d %s"; if (print) printf(tfmt, tmp->tm_hour, tmp->tm_min, tmp->tm_sec, timek(tmp)); else printw(tfmt, tmp->tm_hour, tmp->tm_min, tmp->tm_sec, timek(tmp)); } static char inb[30]; static int bufsub; static char buf; /* * Scan stdin for a character. */ scan() { if (print) return; while (read(0, &buf, 1) == 1) { move(23, bufsub + 4); if (buf != ' ' && buf != '\n') { inb[bufsub++] = buf; addch(buf); } else { inb[bufsub] = '\0'; printw(" "); donin++; bufsub = 0; } refresh(); } } /* * Process the user's request, if it's there. */ dorequest() { REAL dt; if ( ! donin) return(-1); else donin = 0; switch (inb[0]) { case 'e': /* enable elapsed time mode */ if ((deltat * 1440.0 * 60.0) < 999.0) eltime++; break; case 'n': if (inb[1]) switch (inb[1]) { case 'e': /* disable elapsed time mode */ eltime = 0; move(pt[SEPMAT].row+7, SECCOL); printw(" "); break; case 'l': /* disable looping */ loop = 0; break; } else { settime(); return(0); /* don't ticktock */ } break; case 'q': die(); case 'r': /* reverse the flow of time */ if ( ! eltime) deltat = -deltat; break; default: /* look for a new delta t */ dt = parsint(inb); if (eltime && (dt * 1440.0 * 60.0 > 999.0)) dt = 0.0; if (dt) /* * Running in reverse? */ if (deltat < 0.0) deltat = -dt; else deltat = dt; break; } return(1); } SHAR_EOF if test -f 'main.c' then echo shar: over-writing existing file "'main.c'" fi cat << \SHAR_EOF > 'main.c' #include <fcntl.h> #include <signal.h> #include "p.h" main(argc, argv) char **argv; { register int wanderer; register struct planet *planp; register int bump; int die(); while (argc-- > 1) { if (strcmp(argv[1], "-p") == NULL) argv++, print++; if (strncmp(argv[1], "-l", 2) == NULL) { /* * Look for the loop increment. */ deltat = parsint(&argv[1][2]); if (deltat == 0.0) deltat = 1.0; argv++, loop = 1; } if (strcmp(argv[1], "-n") == NULL) argv++, ritenow++; if (strcmp(argv[1], "-e") == NULL) argv++, eltime++; } /* * Who want's to sleep this long anyway? */ if (eltime) if ((long) (deltat * 1440.0 * 60.0) > 999) eltime = 0; if (ritenow) settime(); else gettime(); if ( ! print) { initscr(); erase(); sleep(1); refresh(); signal(SIGINT, die); signal(SIGHUP, die); signal(SIGBUS, die); raw(); noecho(); flags = fcntl(0, F_GETFL, O_NDELAY); fcntl(0, F_SETFL, O_NDELAY); clear(); sleep(1); move(23, 0); printw(">>> "); prthdr(); prtshdr(); for (wanderer = MERCURY; wanderer <= PLUTO; wanderer++) prtname(wanderer); prtname(SUN); prtname(MOON); prtname(TIME); move(0, 0); refresh(); } do { /* * Process user requests only in interactive mode. */ if ( ! print) bump = dorequest(); gestalt(); /* do everything */ for (wanderer = MERCURY; wanderer <= PLUTO; wanderer++) { if ( ! print) bump = dorequest(); if (print) prtname(wanderer); /* * Compute everything about a planet. */ planp = plandata(wanderer); /* * Save equatorial coordinates for angular * separation calculations. */ wandeq[wanderer] = planp->equat; /* * Print the planetary data. */ prtplan(planp, wanderer); /* * Scan for user input. */ scan(); } if (print) prtname(SUN); sun.equat = *ectoeq(&sun.eclipt); sun.hzn = *eqtohzn(&sun.equat); sun.rs = *sunriseset(); prtsun(); if (print) prtname(MOON); moondata(); prtmoon(); if (print) prtname(TIME); prtimes(); eclipse(); prtsep(); scan(); if ( ! print) { move(pt[SEPMAT].row+7, 0); printw("%#-9.3g", deltat); refresh(); } if ( ! print) { if (bump != 0) ticktock(); } else ticktock(); } while (loop); if (print) putchar('\n'); else { move(23, 0); refresh(); } if ( ! print) die(); else exit(0); } /* * Parse the interval pointed to by arg. */ REAL parsint(sptr) register char *sptr; { register char chr; register indx; REAL dt; dt = atof(sptr); indx = strlen(sptr) - 1; chr = sptr[indx]; /* * Delta interval is in days; modify based * on optional units. */ switch (chr) { case 's': if (dt < 1.0) dt = 1.0; dt /= 1440.0 * 60.0; break; case 'm': dt /= 1440.0; break; case 'h': dt /= 24.0; break; case 'w': dt *= 7.0; break; case 'y': dt *= 365.2422; break; default: dt = 0.0; break; } return(dt); } die() { signal(SIGINT, SIG_IGN); signal(SIGHUP, SIG_IGN); signal(SIGBUS, SIG_IGN); move(COLS - 1, LINES - 1); echo(); noraw(); fcntl(0, F_SETFL, flags); endwin(); exit(0); } SHAR_EOF if test -f 'misc.c' then echo shar: over-writing existing file "'misc.c'" fi cat << \SHAR_EOF > 'misc.c' #include "p.h" /* * Do everything necessary to keep everything * up to date. */ gestalt() { tempus(); /* fugit */ earthorbit(); /* don't loose it */ sunorbit(); /* or this, either */ } /* * Keep the now structure up to date. */ tempus() { now.gmt = now.lt; /* copy over local time */ lttogmt(&now.gmt, HEREZONE); /* generate GMT */ now.gst = now.gmt; gmttogst(&now.gst); /* generate GST */ now.lst = now.gst; gsttolst(&now.lst, HERELONG); /* generate LST */ now.jd = julianday(&now.lt); /* * In case it's not there. */ now.lt.tm_wday = dayofwk(now.jd); } /* * Go through the routines to convert LST to LT using * the now structure. */ struct tm * lsttolt(tmp) register struct tm *tmp; { static struct tm tm; tm = now.lt; /* use the now date */ tm.tm_hour = tmp->tm_hour; tm.tm_min = tmp->tm_min; tm.tm_sec = tmp->tm_sec; lsttogst(&tm, HERELONG); gsttogmt(&tm); gmttolt(&tm, HEREZONE); return(&tm); } /* * Count the number of days from the epoch, Jan. 0, 1980. * Convert date to Julian day number, then subtract the bias for * Jan. 0, 1980. Note that the number of days could be negative. */ REAL epoch80(tmp) register struct tm *tmp; { REAL e80; /* * First convert to Julian day. */ e80 = julianday(tmp); /* * Then subtract the bias for the Jan 0, 1980, epoch. */ e80 -= 2444238.5; return(e80); } /* * Return the number of days in the year of interest. */ dysize(yr) { if ((yr % 4) == 0 && yr % 100 && (yr % 400) == 0) return(366); return(365); } /* * Compute and return a floating remainder. */ REAL flrem(a, b) REAL a, b; { REAL c; long d; c = a / b; d = (long) c; c -= (REAL) d; return(c * b); } /* * Compute and return day of the month as a floating number. */ REAL flday(tmp) register struct tm *tmp; { REAL day; day = hmstohr(tmp) / 24.0; day += (REAL) tmp->tm_mday; return(day); } /* * Fix the floating argument day into appropriate fields * of the time structure. */ fixday(tmp, day) register struct tm *tmp; REAL day; { REAL hours; tmp->tm_mday = (int) day; /* * Extract the fraction of day and convert to REAL hours. */ hours = flrem(day, 1.0) * 24.0; tmp->tm_hour = (int) hours; hrtohms(tmp, hours); } /* * Compute arctan of two arguments resolving the ambiguity * as shown in Fig. 10, pg. 45. */ REAL atan2a(y, x) REAL y, x; { REAL a, lo, hi; a = (REAL) atan(y/x); if (y >= 0.0 && x >= 0.0) { lo = 0.0; hi = 90.0; } else if (y < 0.0 && x >= 0.0) { lo = 270.0; hi = 360.0; } else if (y < 0.0 && x < 0.0) { lo = 180.0; hi = 270.0; } else { lo = 90.0; hi = 180.0; } /* * Go to degrees. */ a *= RAD; if (lo <= a && a < hi) return(a/RAD); if (a < lo) while (a < lo) a += 90.0; else while (a >= hi) a -= 90.0; return(a/RAD); } /* * Compute an object's eccentric anomaly using Kepler's method. * This method is valid for e's < 0.1. */ REAL kepler(M, e) REAL M, e; { REAL E, delta, ad; iter = 0; E = M; while (1) { if (iter++ > 100) /* not converging */ break; delta = E - e * (REAL) sin(E) - M; if (delta < 0.0) ad = -delta; else ad = delta; if (ad < 10e-6) break; E -= delta / (1.0 - e * (REAL) cos(E)); } return(E); } /* * Bring the argument angle in to the range 0 - 359.999 ... * This routine works for negative angles. */ REAL dcanon(angle) REAL angle; { if (angle > 0.0) while (angle >= 360.0) angle -= 360.0; else while (angle < 0.0) angle += 360.0; return(angle); } /* * Make the time argument canonical. */ REAL tcanon(tim) REAL tim; { if (tim > 0.0) while (tim >= 24.0) tim -= 24.0; else while (tim < 0.0) tim += 24.0; return(tim); } SHAR_EOF if test -f 'p.h' then echo shar: over-writing existing file "'p.h'" fi cat << \SHAR_EOF > 'p.h' #include <stdio.h> #include <time.h> #include <math.h> #include <curses.h> /* * Compute and display various data on the Sun, Moon, and * planets. The source book for this program is * * "Practical astronomy with your calculator" * by * Peter Duffet-Smith * * R. C. Kukuk * 11/83 */ typedef double REAL; extern int dmsize[]; extern char month[12][3]; extern char days[7][3]; /* * The time structure, tm, is used as the basic data structure for * holding time and date information. The following differences * are noted: * * tm_year contains the actual year. * * tm_mon ranges between 1 and 12. * * tm_yday ranges between 1 and 365 (366). * * tm_isdst contains more than just a daylight savings time flag. * It also contains an indicator describing what kind of time is * currently being stored in the structure. The kinds are listed * below. * * When conversions are made between the various time formats, not all * of the fields may be valid. For example, only the fields relating * to the time of day are valid for GMT, GST, and LST. */ #define DST 1 /* Local daylight savings time */ #define LT 0 /* Local time */ #define LST -1 /* Local siderial time */ #define GST -2 /* Greenwich siderial time */ #define GMT -3 /* Greenwich mean time */ #define NAPLONG 88.1 /* west long is pos */ #define NAPLAT 41.7 /* west lat is pos */ #define NAPELEV 198.0 /* meters */ #define NAPZONE 6.0 /* Central Standard Timezone */ #define HERELONG NAPLONG #define HERELAT NAPLAT #define HERELEV NAPELEV #define HEREZONE NAPZONE #define NORISE 25.0 #define NOSET -1.0 #define PI 3.141592654 #define RAD (360.0/(2.0*PI)) #define NOW 1 #define NOTNOW 0 /* * Coordinate structures. */ struct equat { REAL ra; /* in hours, min, and sec */ REAL dec; }; struct ae { REAL az; REAL el; }; struct eclipt { REAL lat; REAL lon; }; struct deg { char sg; /* sign for neg angles */ int deg; /* or hour angle or ra */ int min; int sec; }; struct rs { REAL risetime; REAL settime; REAL riseaz; REAL setaz; }; /* * CRT screen coordinates for displaying the results. */ struct pt { short row; short col; }; extern struct pt pt[]; /* * Planetary orbital position data structure. */ struct planpos { REAL M; /* Mean anonaly */ REAL v; /* True anomaly */ REAL L; /* Heliocentric longitude */ REAL R; /* Radius vector (AU) */ }; /* * Everything you wanted to know about a planet. */ struct planet { struct equat equat; /* Location */ struct ae hzn; /* Current horizon location */ struct rs rs; /* Rising & setting times and az's */ REAL rho; /* Distance (AU) */ struct deg ltt; /* Light travel time (hr, min) */ REAL th; /* Angular diameter */ REAL F; /* Phase */ REAL m; /* Apparent brightness */ }; /* * Everything worth knowing about Sol. The ecliptic field must * be kept up to date by calling sunorbit(); */ struct sun { struct equat equat; /* Location */ struct eclipt eclipt; /* "Funny" heliocentric coordinates */ struct ae hzn; /* If you've lost it */ struct rs rs; /* "Sunrise, sunset, swiftly go the days ..." */ REAL th; /* Angular diameter */ REAL dist; /* Distance in AU */ REAL M; /* Needed for Lunar calculations */ }; extern struct sun sun; /* * Everything worth knowing about Luna. */ struct moon { struct equat equat; /* Location */ struct ae hzn; /* If you've lost it */ struct rs rs; /* Moon rise and set */ REAL dist; /* Distance in km */ REAL th; /* Angular diameter */ REAL age; /* Age of phase */ REAL F; /* Phase */ REAL lpp; /* Needed for phase calculations */ REAL Mpm; /* Needed for distance calcs */ REAL Ec; /* " " " " */ REAL heliolon; /* Needed for eclipse calculations */ REAL Np; /* " " " " */ }; extern struct moon moon; /* * Planetary equatorial coordinates; kept and used for angular * separation calculations. */ extern struct equat wandeq[]; /* * Orbital data for Sol. * (Pg. 82) */ #define S_EPSG 278.83354 #define S_OMEGG 282.596403 #define S_E 0.016718 #define S_R0 1.495985e8 #define S_TH0 0.533128 /* * Orbital data for Luna. * (Pg. 140) */ #define L_L0 64.975464 #define L_P0 349.383063 #define L_N0 151.950429 #define L_I 5.145396 #define L_E 0.05490 #define L_TH0 0.5181 #define L_A 384401.0 #define L_PI0 0.9507 /* * Appropriate subscripts. */ #define MERCURY 0 #define VENUS 1 #define MARS 2 #define JUPITER 3 #define SATURN 4 #define URANUS 5 #define NEPTUNE 6 #define PLUTO 7 #define SUN 8 #define MOON 9 #define SEPMAT 10 #define TIME 11 #define SECCOL 9 /* * Planetary orbital data. */ struct pod { REAL Tp; /* Period in tropical years */ REAL epsilon; /* Longitude at epoch (1/0/1980) */ REAL w; /* Longitude of perihelion */ REAL e; /* Eccentricity of orbit */ REAL a; /* Semi-major axis (AU) */ REAL i; /* Inclination of orbit */ REAL omega; /* Longitude of ascending node */ REAL th0; /* Angular size at 1 AU */ REAL A; /* Brightness factor */ }; extern struct pod pod[]; /* * Keep Earth orbital data separate for ease in subscripting the others. */ extern struct pod epod; extern char *name[]; REAL flrem(), flday(), julianday(), hmstohr(), dmstodeg(), epoch80(); REAL s12b(); REAL hatora(), ratoha(); REAL obliq(), kepler(), dcanon(), tcanon(); REAL angsep(), atan2a(), parsint(); double atof(); struct tm *jdtodate(), *lsttolt(); struct ae *eqtohzn(); struct equat *getequat(), *hzntoeq(), *ectoeq(), *moonorbit(); struct eclipt *eqtoec(); struct rs *riseset(), *sunriseset(), *atmref(), *moonriseset(); struct planpos *orbit(); struct planet *plandata(); /* * This is the now structure. Its contents should be kept up to date * by calling tempus(). */ struct now { struct tm lt; struct tm gmt; struct tm gst; struct tm lst; REAL jd; }; extern struct now now; /* * Orbital position of Earth (heliocentric coodinates) computed * from the time in the now structure. Kept up to date by calling * earthorbit(). */ extern struct planpos earth; extern int iter; /* iteration count for Keplerian calculations */ extern int flags; /* file control flags */ extern int print; /* hard copy output flag */ extern int loop; /* keep looping */ extern int ritenow; /* compute results for today's date and time */ extern int request; /* user has request pending */ extern int daylight; /* enable DST conversions */ extern REAL deltat; /* loop increment in REAL days */ extern int eltime; /* sleep away the actual interval */ extern int donin; SHAR_EOF if test -f 'sm.c' then echo shar: over-writing existing file "'sm.c'" fi cat << \SHAR_EOF > 'sm.c' #include "p.h" /* * (Sect. 42; Pg. 80) */ sunorbit() { REAL E, N, D, M, v; D = epoch80(&now.lt); N = (360.0 / 365.2422) * D; N = dcanon(N); M = N + S_EPSG - S_OMEGG; M = dcanon(M); E = kepler(M/RAD, S_E); v = (REAL) sqrt((1.0 + S_E) / (1.0 - S_E)) * (REAL) tan(E/2); v = (REAL) atan(v); v *= 2 * RAD; sun.eclipt.lon = dcanon(v + S_OMEGG); sun.eclipt.lat = 0.0; sun.M = M; /* needed for calculating Luna's orbit */ /* * Do the incidentals. */ sun.dist = (1.0 - S_E * S_E) / (1.0 + S_E * (REAL) cos(v / RAD)); sun.th = S_TH0 * (1.0 + S_E * (REAL) cos(v/RAD)) / (1.0 - S_E * S_E); } /* * (Sect. 45; Pg. 88) */ struct rs * sunriseset() { static struct rs rs; struct rs *rsp, rs1, rs2; struct equat *eqp, eq1, eq2; struct eclipt eclipt; REAL dayfrac; /* * Make a copy of the sun's location. */ eclipt = sun.eclipt; /* * Compute the fraction of the current day so the Sun * can be "moved back" to it's position on the previous * midnight. */ dayfrac = hmstohr(&now.lt) / 24.0; /* * Then move it back. */ eclipt.lon -= dayfrac * 0.985647; eqp = ectoeq(&eclipt); /* * Copy over first set of equatorial coords. */ eq1 = *eqp; /* * Then move the Sun to its position on midnight, tomorrow. */ eclipt.lon += 0.985647; eqp = ectoeq(&eclipt); /* * Copy over second set of equatorial coords. */ eq2 = *eqp; rsp = riseset(&eq1); rs1 = *rsp; rsp = riseset(&eq2); rs2 = *rsp; rs.risetime = (24.07 * rs1.risetime) / (24.07 + rs1.risetime - rs2.risetime); rs.settime = (24.07 * rs1.settime) / (24.07 + rs1.settime - rs2.settime); rs.riseaz = (rs1.riseaz + rs2.riseaz) / 2.0; rs.setaz = (rs1.setaz + rs2.setaz) / 2.0; /* * Make a rough correction to allow for the Sun's * finite diameter and refraction. (+- 5 minutes.) */ rs.risetime -= 5.0/60.0; rs.settime += 5.0/60.0; /* * The rising and setting times are LST. */ return(&rs); } /* * (Sect. 61; Pg. 139) */ struct equat * moonorbit(tmp, nowflag) register struct tm *tmp; { static struct eclipt eclipt; REAL D, l, Mm, N, Ev, Ae, A3, Mpm, Ec; REAL A4, lp, V, lpp, Np, x, y; REAL t; D = epoch80(tmp); l = dcanon(13.176396 * D + L_L0); Mm = dcanon(l - 0.1114041 * D - L_P0); N = dcanon(L_N0 - 0.0529539 * D); Ev = 1.2739 * (REAL) sin(2.0 * (l/RAD - sun.eclipt.lon/RAD) - Mm/RAD); Ae = 0.1858 * (REAL) sin(sun.M / RAD); A3 = 0.37 * (REAL) sin(sun.M / RAD); Mpm = Mm + Ev - Ae - A3; Ec = 6.2886 * (REAL) sin(Mpm / RAD); A4 = 0.214 * (REAL) sin(2.0 * Mpm / RAD); lp = l + Ev + Ec - Ae + A4; V = 0.6583 * (REAL) sin(2.0 * (lp/RAD - sun.eclipt.lat/RAD)); lpp = lp + V; Np = N - 0.16 * (REAL) sin(sun.M / RAD); t = lpp / RAD - Np / RAD; y = (REAL) sin(t) * (REAL) cos(L_I/RAD); x = (REAL) cos(t); eclipt.lon = (REAL) atan2a(y, x) * RAD + Np; x = (REAL) sin(t) * (REAL) sin(L_I/RAD); eclipt.lat = (REAL) asin(x) * RAD; /* * Stash away the other quantities that will be needed * elsewhere. Don't stash the quantities computed as a * result of the moon rise and set calls. */ if (nowflag == NOW) { moon.lpp = lpp; moon.Mpm = Mpm; moon.Np = Np; moon.Ec = Ec; moon.heliolon = eclipt.lon; } /* * Return equatorial coords. */ return(ectoeq(&eclipt)); } /* * (Sect. 63; Pg. 145) */ moondata() { /* * Compute the Moon's orbit. */ moon.equat = *moonorbit(&now.lt, NOW); /* * Apply precession correction. */ precess(&moon.equat); /* * And then convert to horizon coords. */ moon.hzn = *eqtohzn(&moon.equat); /* * Do age and Fase. */ moon.age = dcanon(moon.lpp - sun.eclipt.lon) / 13.0; moon.F = 0.5 * (1.0 - (REAL) cos(moon.age * 13.0 / RAD)); /* * Find the distance in km. */ moon.dist = L_A * (1.0 - L_E * L_E) / (1.0 + L_E * (REAL) cos((moon.Mpm + moon.Ec)/RAD)); /* * Now the moon's apparent diameter, in degrees. */ moon.th = L_TH0 * L_A / moon.dist; /* * Last, but not least, compute moonrise and moonset. */ moon.rs = *moonriseset(); } #define NOECL 0 #define LUNAR 1 #define SOLAR 2 #define MAYBE 3 #define MUST 4 /* * Look at the Solar and Lunar coordinates and determine if a * Lunar or Solar eclipse is possible. If so, display an * appropriate message. If a message has been displayed, and * the conditions for the eclipse no longer hold, take down * the message. */ eclipse() { static int eclmsg; static REAL svdelta; REAL londiff, nodediff; register int ecltype, echance; /* * Make the angles canonical, and make the node angle * within 0 <= nodediff < 180. */ londiff = dcanon(moon.heliolon - sun.eclipt.lon); nodediff = dcanon(moon.lpp - moon.Np); if (nodediff >= 180.0) nodediff -= 180.0; /* * If an eclipse is possible, which one will it be? */ if (londiff >= 353.0 || londiff <= 7.0) ecltype = SOLAR; else if (londiff >= 173.0 && londiff <= 187.0) ecltype = LUNAR; else ecltype = NOECL; /* * Slow down the passage of time in order to take a * closer look. Or speed up because the configuration * is past. * The nodediff test is needed so we don't slow down * every opposition or conjunction, in spite of the * location of the nodes. */ if (ecltype != NOECL && nodediff <= 20.0) { if (svdelta == 0.0) { svdelta = deltat; /* * Back up before slowing down. */ deltat = -(deltat * 0.75); ticktock(); gestalt(); deltat = svdelta; /* * Then decide how much to slow down. */ if (svdelta < 0.1) deltat /= 10.0; if (svdelta >= 0.1) deltat /= 10.0; if (svdelta >= 1.0) deltat /= 4.0; if (svdelta >= 7.0) deltat /= 4.0; } } else { /* * Full speed ahead. */ if (svdelta) { /* * If we've reversed directions during the * slowed interval, carry over the new * sign. */ if ((svdelta * deltat) < 0.0) /* different signs */ deltat = -svdelta; else deltat = svdelta; svdelta = 0.0; } } /* * But is one the Moon's orbital nodes in the right place? */ if (ecltype == SOLAR) { echance = NOECL; if (nodediff <= 18.52) echance = MAYBE; if (nodediff <= 15.52) echance = MUST; } if (ecltype == LUNAR) { echance = NOECL; if (nodediff <= 12.25) echance = MAYBE; if (nodediff <= 9.5) echance = MUST; } /* * Now generate a message based on the flags. */ if ( ! print) { if (ecltype != NOECL && echance != NOECL) { if (eclmsg == 0) { standout(); move(pt[SUN].row, 63); printw("%s", ecltype == LUNAR ? "Lunar " : "Solar "); move(pt[SUN].row, 69); printw("Eclipse"); eclmsg = 1; standend(); } move(pt[SUN].row-1, 67); if (echance == MAYBE) { standout(); printw("Maybe"); standend(); } else printw(" "); } else if (eclmsg) { move(pt[SUN].row-1, 67); printw(" "); move(pt[SUN].row, 63); printw(" "); move(pt[SUN].row, 69); printw(" "); eclmsg = 0; } } } /* * (Sect. 66; Pg. 150) * * (Oh, by the way, there's something wrong with the algorithm, * or maybe the implementation of it ...) */ struct rs * moonriseset() { static struct rs mrs; struct equat eq1, eq2; struct rs rs1, rs2; struct tm *tmp; REAL jday; /* * If the Moon is very close the the horizon, the * rising and setting times will be way off; don't * do the calculations. Just return the old values. * * But make at least one calculation first time through. */ if (moon.hzn.el > -8.0 && moon.hzn.el < 8.0 && mrs.setaz != 0.0) return(&mrs); /* * First get today's julian date. */ jday = now.jd; /* * Convert it to midnight yesterday. */ jday -= 0.5; jday -= flrem(jday, 1.0); jday += 0.5; /* * Make a date out of it. */ tmp = jdtodate(jday); /* * Find the Moon's coords at that time. */ eq1 = *moonorbit(tmp, NOTNOW); /* * Advance the clock 12 hours and redo the calc. */ tmp = jdtodate(jday + 0.5); eq2 = *moonorbit(tmp, NOTNOW); /* * Apply precession correction. */ precess(&eq1); precess(&eq2); /* * Correct the coordinates for geocentric parallax. */ geopalx(&eq1); geopalx(&eq2); /* * Compute the rising and setting times and az's. */ rs1 = *riseset(&eq1); rs2 = *riseset(&eq2); /* * Interpolate to find the moon rise and set times. */ mrs.risetime = (12.03 * rs1.risetime) / (12.03 + rs1.risetime - rs2.risetime); mrs.settime = (12.03 * rs1.settime) / (12.03 + rs1.settime - rs2.settime); /* * Just average to find azimuths. */ mrs.riseaz = (rs1.riseaz + rs2.riseaz) / 2.0; mrs.setaz = (rs1.setaz + rs2.setaz) / 2.0; /* * Make a rough correction to allow for the Moon's * finite diameter. (+- 2 minutes.) */ mrs.risetime -= 2.0/60.0; mrs.settime += 2.0/60.0; /* * The rising and setting times are LST. */ return(&mrs); } SHAR_EOF if test -f 'time.c' then echo shar: over-writing existing file "'time.c'" fi cat << \SHAR_EOF > 'time.c' #include "p.h" /* * Compute the day of the year. Put it in the time structure. * And return it. * (Sect. 3; Pg. 6) */ dayofyr(tmp) register struct tm *tmp; { register int day, mon; day = 1; mon = tmp->tm_mon; /* * Leap year */ if (dysize(tmp->tm_year) == 366 && mon >= 3) day += 1; while(--mon) day += dmsize[mon-1]; tmp->tm_yday = day + tmp->tm_mday; return(tmp->tm_yday); } /* * (Sect. 4; Pg. 9) */ REAL julianday(tmp) register struct tm *tmp; { short y, m; register long A, B, C, D; REAL day; y = tmp->tm_year; m = tmp->tm_mon; if (m <= 2) { y--; m += 12; } /* * Beware of 1582 Oct 15 */ if (tmp->tm_year <= 1582 && tmp->tm_mon <= 10 && tmp->tm_mday <= 15) B = 0; else { A = y / 100; B = 2 - A + (A / 4); } C = (long) (365.25 * (REAL) y); D = (long) (30.6001 * ((REAL) m + 1.0)); /* * Get the floating equivalent of the day of the month. */ day = flday(tmp); day += 1720994.5; day += (REAL) B + (REAL) C + (REAL) D; return(day); } /* * (Sect. 5; Pg. 11) */ struct tm * jdtodate(jday) REAL jday; { static struct tm jtm; REAL F, d; register long A, B, C, D, E, G, H, I; jday += 0.5; I = (long) jday; F = flrem(jday, 1.0); if (I > 2299160) { A = (long) (((REAL) I - 1867216.25) / 36524.25); B = I + 1 + A - (A / 4); } else B = I; /* note the misprint in the book */ C = B + 1524; D = (long) (((REAL) C - 122.1) / 365.25); E = (long) (365.25 * (REAL) D); G = (long) (((REAL) C - (REAL) E) / 30.6001); H = (long) ((REAL) G * 30.6001); d = (REAL) C - (REAL) E + F - H; if (G <= 13) jtm.tm_mon = G - 1; else jtm.tm_mon = G - 13; if (jtm.tm_mon <= 2) jtm.tm_year = D - 4715; else jtm.tm_year = D - 4716; fixday(&jtm, d); /* * Compute the day of the year into the structure. */ dayofyr(&jtm); /* * Also determine the day of the week. */ jtm.tm_wday = dayofwk(jday); /* * Possible adjust for DST. */ dstime(&jtm); return(&jtm); } /* * (Sect. 6; Pg. 12) */ dayofwk(jday) REAL jday; { REAL A; register int wkday; A = (jday + 1.5) / 7.0; wkday = (int) (flrem(A, 1.0) * 7.0); return(wkday); } /* * (Sect. 7; Pg. 13) */ REAL hmstohr(tmp) register struct tm *tmp; { REAL sec, min, hr; sec = (REAL) tmp->tm_sec; min = (REAL) tmp->tm_min; hr = (REAL) tmp->tm_hour; sec /= 60.0; min += sec; min /= 60.0; hr += min; return(hr); } /* * (Sect. 7; Pg. 13) */ REAL dmstodeg(degp) register struct deg *degp; { REAL sec, min, hr; sec = (REAL) degp->sec; min = (REAL) degp->min; hr = (REAL) degp->deg; sec /= 60.0; min += sec; min /= 60.0; hr += min; return(hr); } /* * (Sect. 8; Pg. 14) */ hrtohms(tmp, hours) register struct tm *tmp; REAL hours; { REAL F; tmp->tm_hour = (int) hours; F = flrem(hours, 1.0); F *= 60.0; tmp->tm_min = (int) F; F = flrem(F, 1.0); F *= 60.0; tmp->tm_sec = (int) F; } /* * (Sect. 8; Pg. 14) */ degtodms(degp, hours) register struct deg *degp; REAL hours; { REAL F; if (hours < 0.0) { hours = -hours; degp->sg = '-'; } else degp->sg = ' '; degp->deg = (int) hours; F = flrem(hours, 1.0); F *= 60.0; degp->min = (int) F; F = flrem(F, 1.0); F *= 60.0; degp->sec = (int) F; } /* * (Sect. 9; Pg. 16) */ lttogmt(tmp, zone) register struct tm *tmp; REAL zone; { REAL dectime; dectime = hmstohr(tmp); /* * Correct for DST, if necessary. Use the dst flag in the * now structure. */ if (now.lt.tm_isdst > 0) zone -= 1.0; dectime += zone; dectime = tcanon(dectime); hrtohms(tmp, dectime); tmp->tm_isdst = GMT; } /* * (Sect. 10; Pg. 18) */ gmttolt(tmp, zone) register struct tm *tmp; REAL zone; { REAL dectime; dectime = hmstohr(tmp); /* * Correct for DST, if necessary. Use the dst flag in the * now structure. */ if (now.lt.tm_isdst > 0) { zone -= 1.0; tmp->tm_isdst = DST; } else tmp->tm_isdst = LT; dectime -= zone; dectime = tcanon(dectime); hrtohms(tmp, dectime); } /* * (Sect. 12; Pg. 20) */ #define S12A 0.0657098 #define S12C 1.002738 #define S12D 0.997270 gmttogst(tmp) register struct tm *tmp; { REAL ydays, gmthours; REAL T0, B; ydays = (REAL) dayofyr(tmp); ydays *= S12A; B = s12b(tmp); T0 = ydays - B; gmthours = hmstohr(tmp) * S12C; T0 += gmthours; T0 = tcanon(T0); hrtohms(tmp, T0); tmp->tm_isdst = GST; } /* * Compute B in section 12. */ REAL s12b(tmp) register struct tm *tmp; { REAL B, JD, R, S, T, U; static struct tm jan0; jan0.tm_year = tmp->tm_year; jan0.tm_mon = 1; /* January 0 */ JD = julianday(&jan0); S = JD - 2415020.0; T = S / 36525.0; R = 6.6460656 + (2400.051262 * T) + (0.00002581 * T * T); U = R - (24.0 * (REAL) (tmp->tm_year - 1900)); B = 24.0 - U; return(B); } /* * (Sect. 13; Pg. 22) */ gsttogmt(tmp) register struct tm *tmp; { REAL gmthours, gsthours, ydays, T0, B; ydays = (REAL) dayofyr(tmp); ydays *= S12A; B = s12b(tmp); T0 = ydays - B; T0 = tcanon(T0); gsthours = hmstohr(tmp); gsthours -= T0; gsthours = tcanon(gsthours); gmthours = gsthours * S12D; hrtohms(tmp, gmthours); tmp->tm_isdst = GMT; } /* * (Sect. 14; Pg. 24) */ gsttolst(tmp, longitude) register struct tm *tmp; REAL longitude; { REAL dectime, diff; dectime = hmstohr(tmp); diff = longitude / 15.0; dectime -= diff; /* W long is positive, W is negative */ dectime = tcanon(dectime); hrtohms(tmp, dectime); tmp->tm_isdst = LST; } /* * (Sect. 15; Pg. 25) */ lsttogst(tmp, longitude) register struct tm *tmp; REAL longitude; { REAL dectime, diff; dectime = hmstohr(tmp); diff = longitude / 15.0; dectime += diff; /* W long is positive, W is negative */ dectime = tcanon(dectime); hrtohms(tmp, dectime); tmp->tm_isdst = GST; } /* * Update the local time in the now structure. Or sleep the equivalent * amount if we're running in real time. */ ticktock() { short tpast, tickint, ticktime; long cursec; REAL jday, dt; double fabs(); /* * Sleep the interval away. */ if (eltime) { /* * Convert julian day to seconds. */ jday = julianday(&now.lt) * 1440.0 * 60.0; /* * Find the number of seconds into the day. */ jday = flrem(jday, 1440.0 * 60.0); cursec = (long) (jday + 0.5); tickint = (unsigned) ((REAL) fabs(deltat) * 1440.0 * 60.0); tpast = cursec % tickint; ticktime = tickint - tpast; /* * Beware! We don't have one second resolution. */ if (ticktime == 1) ticktime += tickint; /* * Keep the REAL equivalent of the ticktime in deltat * units. */ dt = ((REAL) ticktime) / (1440.0 * 60.0); /* * Break up the sleep to look for input. */ while (ticktime > 0) { move(pt[SEPMAT].row+7, SECCOL); printw("(%d sec) ", ticktime); move(0, 0); refresh(); if (ticktime > 5) { sleep(5); ticktime -= 5; } else { sleep(ticktime); break; } scan(); /* * If a request was processed, break out of the * sleep loop; this request may affect sleep time. */ if (dorequest() >= 0) break; } } else dt = deltat; jday = julianday(&now.lt) + dt; /* * Copy back the entire structure. */ now.lt = *jdtodate(jday); } /* * Generate a new localtime for the now structure. */ settime() { struct tm *localtime(); long tvec; time(&tvec); now.lt = *localtime(&tvec); now.lt.tm_mon++; now.lt.tm_yday++; now.lt.tm_year += 1900; } /* * The following table is used for 1974 and 1975 and * gives the day number of the first day after the Sunday of the * change. */ static struct { int daylb; int dayle; } daytab[] = { 5, 333, /* 1974: Jan 6 - last Sun. in Nov */ 58, 303, /* 1975: Last Sun. in Feb - last Sun in Oct */ }; /* * Determine if the time pointed to by tmp should be set to DST. * If so, do it. */ dstime(tmp) register struct tm *tmp; { register int dayno, daylbegin, daylend; /* * Don't try to convert to DST for times and dates earlier * than 1/1/1900. */ if (tmp->tm_year < 1900) return; dayno = tmp->tm_yday; daylbegin = 119; /* last Sun in Apr */ daylend = 303; /* Last Sun in Oct */ if(tmp->tm_year == 1974 || tmp->tm_year == 1975) { daylbegin = daytab[tmp->tm_year-1974].daylb; daylend = daytab[tmp->tm_year-1974].dayle; } daylbegin = sunday(tmp, daylbegin); daylend = sunday(tmp, daylend); if(daylight && (dayno > daylbegin || (dayno == daylbegin && tmp->tm_hour >= 2)) && (dayno < daylend || (dayno == daylend && tmp->tm_hour < 1))) { tmp->tm_isdst = DST; } else tmp->tm_isdst = LT; } /* * The argument is a 0-origin day number. * The value is the day number of the first * Sunday on or after the day. */ static int sunday(t, d) register struct tm *t; register int d; { if(d >= 58) d += dysize(t->tm_year) - 365; return(d - (d - t->tm_yday + t->tm_wday + 700) % 7); } SHAR_EOF if test -f 'wand.c' then echo shar: over-writing existing file "'wand.c'" fi cat << \SHAR_EOF > 'wand.c' #include "p.h" /* * (Sect. 31; Pg. 52) */ REAL angsep(eqp1, eqp2) register struct equat *eqp1, *eqp2; { REAL alph1, alph2, del1, del2; REAL d; /* * Hours to degrees to rads. */ alph1 = eqp1->ra * 15.0 / RAD; alph2 = eqp2->ra * 15.0 / RAD; del1 = eqp1->dec / RAD; del2 = eqp2->dec / RAD; d = (REAL) sin(del1) * (REAL) sin(del2) + (REAL) cos(del1) * (REAL) cos(del2) * (REAL) cos(alph1 - alph2); d = (REAL) acos(d) * RAD; return(d); } /* * (Sect. 32; Pg. 54) */ struct rs * riseset(eqp) register struct equat *eqp; { REAL Ar, As, alpha, delta, H; static struct rs rs; struct rs *drsp; delta = eqp->dec; alpha = eqp->ra; /* alpha is in hours */ Ar = (REAL) sin(delta/RAD) / (REAL) cos(HERELAT/RAD); rs.riseaz = rs.setaz = -1.0; if (Ar <= -1.0) { norise: rs.risetime = rs.settime = NORISE; return(&rs); } if (Ar >= 1.0) { noset: rs.risetime = rs.settime = NOSET; return(&rs); } Ar = (REAL) acos(Ar) * RAD; As = 360.0 - Ar; rs.riseaz = Ar; rs.setaz = As; H = -(REAL) tan(HERELAT/RAD) * (REAL) tan(delta/RAD); if (H >= 1.0) goto norise; if (H <= -1.0) goto noset; H = (REAL) acos(H) * RAD; H /= 15.0; /* degrees to hour angle */ /* * Now correct for atmospheric refraction. */ drsp = atmref(eqp); rs.riseaz -= drsp->riseaz; rs.setaz += drsp->setaz; rs.risetime = tcanon(24.0 + alpha - H - drsp->risetime); rs.settime = tcanon(alpha + H + drsp->settime); /* * The rising and setting times are LST. */ return(&rs); } /* * Correct equatorial coordinates for precession. * (Sect. 33, Pg. 58) */ precess(eqp) register struct equat *eqp; { REAL da, dd; REAL N; REAL d, a; a = eqp->ra / RAD; /* hours to rads */ d = eqp->dec / RAD; /* degrees to rads */ /* * Use epoch 1950 precession values. */ N = (REAL) now.lt.tm_year - 1950.0; /* forget fractions */ #define ATERM (3.07327 / 3600.0) /* seconds to hours */ #define BTERM (1.33617 / 3600.0) /* seconds to hours */ #define CTERM (20.0426 / 3600.0) /* arcseconds to degrees */ /* * Compute the deltas. */ da = (ATERM + BTERM * (REAL) sin(a) * (REAL) tan(d)) * N; dd = CTERM * (REAL) cos(a) * N; eqp->ra += da; /* in hours */ eqp->dec += dd; /* in degrees */ } /* * (Sect. 34; Pg. 60) */ struct rs * atmref(eqp) register struct equat *eqp; { static struct rs drs; REAL x, y, psi, dA, dt; x = 0.566667; psi = (REAL) sin(HERELAT/RAD) / (REAL) cos(eqp->dec/RAD); psi = (REAL) acos(psi) * RAD; dA = (REAL) tan(x/RAD) / (REAL) tan(psi/RAD); dA = (REAL) asin(dA) * RAD; y = (REAL) sin(x/RAD) / (REAL) sin(psi/RAD); y = (REAL) asin(y) * RAD; dt = (240.0 * y) / (REAL) cos(eqp->dec/RAD); dt /= 3600.0; drs.riseaz = drs.setaz = dA; drs.risetime = drs.settime = dt; return(&drs); } /* * Correct the equatorial coordinates of the Moon for geocentric * parallax. * (Sect. 35, 36; Pg. 63, 66) */ geopalx(eqp) register struct equat *eqp; { REAL u, hp, rhosinphi, rhocosphi; REAL dalpha, deltap, r; REAL n, d, H, Hp; /* * First find rhosinphi and rhocosphi. */ u = (REAL) tan(HERELAT/RAD); u = (REAL) atan(0.996647 * u); hp = HERELEV / 6378140.0; rhosinphi = 0.996647 * (REAL) sin(u) + hp * (REAL) sin(HERELAT/RAD); rhocosphi = (REAL) cos(u) + hp * (REAL) cos(HERELAT/RAD); /* * Then use them in the parallax calculations. */ H = ratoha(eqp->ra) * 15.0; r = moon.dist / 6378.16; n = rhocosphi * (REAL) sin(H/RAD); d = r * (REAL) cos(eqp->dec/RAD) * rhocosphi * (REAL) cos(H/RAD); /* * Delta alpha to degrees, then hours; then apply it. */ dalpha = (REAL) atan(n/d) * RAD; eqp->ra -= dalpha / 15.0; Hp = H + dalpha; n = r * (REAL) sin(eqp->dec/RAD) - rhosinphi; d = r * (REAL) cos(eqp->dec/RAD) * (REAL) cos(H/RAD) - rhocosphi; deltap = (REAL) atan((REAL) cos(Hp/RAD) * n / d); /* * Replace the old declination with the new. */ eqp->dec = deltap; } /* * Compute the planetary position data for Earth. */ earthorbit() { earth = *orbit(&epod); } /* * Compute the planetary data for the given planet number. * (Sect. 50, 53, 54, 56; Pp. 98, 113, 115, 118) */ struct planet * plandata(plano) { REAL a, x, y, psi, lp, rp, A, rho; static struct planet planet; register struct planpos *ppos; struct eclipt ecl; struct equat *eqp; /* * Compute the heliocentric coordinates for the requested * planet. Earth has already been done. */ ppos = orbit(&pod[plano]); a = (ppos->L - pod[plano].omega) / RAD; psi = (REAL) sin(a) * (REAL) sin(pod[plano].i/RAD); psi = (REAL) asin(psi) * RAD; y = (REAL) sin(a) * (REAL) cos(pod[plano].i/RAD); x = (REAL) cos(a); a = (REAL) atan2a(y, x) * RAD; lp = a + pod[plano].omega; rp = ppos->R * (REAL) cos(psi/RAD); /* * Different geometry for inner and outer planets. */ if (plano == MERCURY || plano == VENUS) { a = (earth.L - lp) / RAD; A = (rp * (REAL) sin(a)) / (earth.R - rp * (REAL) cos(a)); A = (REAL) atan(A) * RAD; ecl.lon = dcanon(180.0 + earth.L + A); ecl.lat = (rp * (REAL) tan(psi/RAD) * (REAL) sin((ecl.lon - lp)/RAD)) / (earth.R * (REAL) sin((lp - earth.L)/RAD)); ecl.lat *= RAD; } else { a = (lp - earth.L) / RAD; x = (earth.R * (REAL) sin(a)) / (rp - earth.R * (REAL) cos(a)); x = (REAL) atan(x) * RAD + lp; ecl.lon = dcanon(x); x = (rp * (REAL) tan(psi/RAD) * (REAL) sin((ecl.lon - lp)/RAD)) / (earth.R * (REAL) sin(a)); ecl.lat = (REAL) atan(x) * RAD; } /* * Convert heliocentric (ecliptic) coordinates to * equatorial ones. */ eqp = ectoeq(&ecl); /* * Apply precession correction. */ precess(eqp); /* * Stash the equatorial coords. */ planet.equat.ra = eqp->ra; planet.equat.dec = eqp->dec; /* * Compute the horizon coordinates. */ planet.hzn = *eqtohzn(eqp); /* * Compute the rise and set times and azimuths. * (The times are LST.) */ planet.rs = *riseset(eqp); /* * Determine the distance. */ rho = earth.R * earth.R + ppos->R * ppos->R - 2.0 * earth.R * ppos->R * (REAL) cos(ppos->L/RAD - earth.L/RAD); rho = (REAL) sqrt(rho); planet.rho = rho; /* * Do light travel time. */ degtodms(&planet.ltt, rho * 0.1386); /* * Do angular diameter. */ planet.th = pod[plano].th0 / rho; /* * Determine the phase. */ planet.F = 0.5 * (1.0 + (REAL) cos((ecl.lon - ppos->L) / RAD)); /* * And last, but not least, determine the apparent * brightness. */ x = (ppos->R * rho) / (pod[plano].A * (REAL) sqrt(planet.F)); planet.m = 5.0 * (REAL) log10(x) - 26.7; return(&planet); } struct planpos * orbit(plp) register struct pod *plp; { static struct planpos scr; REAL e80, N, eps, w, e, a; e80 = epoch80(&now.lt); eps = plp->epsilon; e = plp->e; w = plp->w; a = plp->a; N = (360.0 / 365.2422) * (e80 / plp->Tp); N = dcanon(N); scr.M = N + eps - w; scr.L = N + (360.0/PI) * e * (REAL) sin(scr.M/RAD) + eps; scr.L = dcanon(scr.L); scr.v = scr.L - w; scr.R = a * (1.0 - e * e) / (1.0 + e * (REAL) cos(scr.v/RAD)); return(&scr); } #ifdef KEEPOUT /* * This function is under construction. Invoke at your own risk. */ orbit(plp) register struct pod *plp; { REAL E, N, D, M, v; D = epoch80(&now.lt); N = (360.0 / 365.2422) * (D / plp->Tp); N = dcanon(N); M = N + plp->epsilon - plp->w; M = dcanon(M); E = kepler(M/RAD, plp->e); v = (REAL) sqrt((1.0 + plp->e) / (1.0 - plp->e)) * (REAL) tan(E/2); v = (REAL) atan(v); v *= 2 * RAD; eclipt.lon = dcanon(v + plp->epsilon); eclipt.lat = 0.0; return(&eclipt); } #endif SHAR_EOF # End of shell archive exit 0