ecd@cs.umn.edu@ncs-med.UUCP (Elwood C. Downey) (03/11/90)
Posting-number: Volume 11, Issue 6
Submitted-by: ecd@cs.umn.edu@ncs-med.UUCP (Elwood C. Downey)
Archive-name: ephem4.12/part05
# This is the first line of a "shell archive" file.
# This means it contains several files that can be extracted into
# the current directory when run with the sh shell, as follows:
# sh < this_file_name
# This is file 5.
echo x srch.c
sed -e 's/^X//' << 'EOFxEOF' > srch.c
X/* this file contains functions to support iterative ephem searches.
X * we support several kinds of searching and solving algorithms.
X * values used in the evaluations come from the field logging flog.c system.
X * the expressions being evaluated are compiled and executed from compiler.c.
X */
X
X#include <stdio.h>
X#include <math.h>
X#include "screen.h"
X
Xextern char *strcpy();
X
Xstatic int (*srch_f)();
Xstatic int srch_tmscalled;
Xstatic char expbuf[NC]; /* [0] == '\0' when expression is invalid */
Xstatic double tmlimit = 1./60.; /* search accuracy, in hrs; def is one minute */
X
X
Xsrch_setup()
X{
X int srch_minmax(), srch_solve0(), srch_binary();
X static char *chcs[] = {
X "Find extreme", "Find 0", "Binary", "New function", "Accuracy",
X "Stop"
X };
X int fn;
X
X /* let op select algorithm, edit, set accuracy
X * or stop if currently searching
X * algorithms require a function.
X */
X fn = 0;
X ask:
X switch (popup(chcs, fn, srch_f ? 6 : 5)) {
X case 0: if (expbuf[0] == '\0')
X set_function();
X srch_f = expbuf[0] ? srch_minmax : 0;
X break;
X case 1: if (expbuf[0] == '\0')
X set_function();
X srch_f = expbuf[0] ? srch_solve0 : 0;
X break;
X case 2: if (expbuf[0] == '\0')
X set_function();
X srch_f = expbuf[0] ? srch_binary : 0;
X break;
X case 3: srch_f = 0; set_function(); fn = 3; goto ask;
X case 4: srch_f = 0; set_accuracy(); fn = 4; goto ask;
X case 5: srch_f = 0; srch_prstate(0); return;
X default: return;
X }
X
X /* new search */
X srch_tmscalled = 0;
X srch_prstate (0);
X}
X
X/* if searching is in effect call the search type function.
X * it might modify *tmincp according to where it next wants to eval.
X * (remember tminc is in hours, not days).
X * if searching ends for any reason it is also turned off.
X * also, flog the new value.
X * return 0 if caller can continue or -1 if it is time to stop.
X */
Xsrch_eval(mjd, tmincp)
Xdouble mjd;
Xdouble *tmincp;
X{
X char errbuf[128];
X int s;
X double v;
X
X if (!srch_f)
X return (0);
X
X if (execute_expr (&v, errbuf) < 0) {
X srch_f = 0;
X f_msg (errbuf);
X } else {
X s = (*srch_f)(mjd, v, tmincp);
X if (s < 0)
X srch_f = 0;
X (void) flog_log (R_SRCH, C_SRCH, v);
X srch_tmscalled++;
X }
X
X srch_prstate (0);
X return (s);
X}
X
X/* print state of searching. */
Xsrch_prstate (force)
Xint force;
X{
X int srch_minmax(), srch_solve0(), srch_binary();
X static (*last)();
X
X if (force || srch_f != last) {
X f_string (R_SRCH, C_SRCHV,
X srch_f == srch_minmax ? "Extrema" :
X srch_f == srch_solve0 ? " Find 0" :
X srch_f == srch_binary ? " Binary" :
X " off");
X last = srch_f;
X }
X}
X
Xsrch_ison()
X{
X return (srch_f != 0);
X}
X
X/* display current expression. then if type in at least one char make it the
X * current expression IF it compiles ok.
X * TODO: editing?
X */
Xstatic
Xset_function()
X{
X static char prompt[] = "Function: ";
X char newexp[NC];
X int s;
X
X f_prompt (prompt);
X fputs (expbuf, stdout);
X c_pos (R_PROMPT, sizeof(prompt));
X
X s = read_line (newexp, PW-sizeof(prompt));
X if (s >= 0) {
X char errbuf[NC];
X if (s > 0 && compile_expr (newexp, errbuf) < 0)
X f_msg (errbuf);
X else
X strcpy (expbuf, newexp);
X }
X}
X
Xstatic
Xset_accuracy()
X{
X static char p[] = "Desired accuracy ( hrs): ";
X int hrs, mins, secs;
X char buf[NC];
X
X f_prompt (p);
X f_time (R_PROMPT, C_PROMPT+18, tmlimit); /* place in blank spot */
X c_pos (R_PROMPT, sizeof(p));
X if (read_line (buf, PW-sizeof(p)) > 0) {
X f_dec_sexsign (tmlimit, &hrs, &mins, &secs);
X f_sscansex (buf, &hrs, &mins, &secs);
X sex_dec (hrs, mins, secs, &tmlimit);
X }
X}
X
X/* use successive paraboloidal fits to find when expression is at a
X * local minimum or maximum.
X */
Xstatic
Xsrch_minmax(mjd, v, tmincp)
Xdouble mjd;
Xdouble v;
Xdouble *tmincp;
X{
X static double base; /* for better stability */
X static double x1, x2, x3; /* keep in increasing order */
X static double y1, y2, y3;
X double xm, a, b;
X
X if (srch_tmscalled == 0) {
X base = mjd;
X x1 = 0.0;
X y1 = v;
X return (0);
X }
X mjd -= base;
X if (srch_tmscalled == 1) {
X /* put in one of first two slots */
X if (mjd < x1) {
X x2 = x1; y2 = y1;
X x1 = mjd; y1 = v;
X } else {
X x2 = mjd; y2 = v;
X }
X return (0);
X }
X if (srch_tmscalled == 2 || fabs(mjd - x1) < fabs(mjd - x3)) {
X /* closer to x1 so discard x3.
X * or if it's our third value we know to "discard" x3.
X */
X if (mjd > x2) {
X x3 = mjd; y3 = v;
X } else {
X x3 = x2; y3 = y2;
X if (mjd > x1) {
X x2 = mjd; y2 = v;
X } else {
X x2 = x1; y2 = y1;
X x1 = mjd; y1 = v;
X }
X }
X if (srch_tmscalled == 2)
X return (0);
X } else {
X /* closer to x3 so discard x1 */
X if (mjd < x2) {
X x1 = mjd; y1 = v;
X } else {
X x1 = x2; y1 = y2;
X if (mjd < x3) {
X x2 = mjd; y2 = v;
X } else {
X x2 = x3; y2 = y3;
X x3 = mjd; y3 = v;
X }
X }
X }
X
X#ifdef TRACEMM
X { char buf[NC];
X sprintf (buf, "x1=%g y1=%g x2=%g y2=%g x3=%g y3=%g",
X x1, y1, x2, y2, x3, y3);
X f_msg (buf);
X }
X#endif
X a = y1*(x2-x3) - y2*(x1-x3) + y3*(x1-x2);
X if (fabs(a) < 1e-10) {
X /* near-0 zero denominator, ie, curve is pretty flat here,
X * so assume we are done enough.
X * signal this by forcing a 0 tminc.
X */
X *tmincp = 0.0;
X return (-1);
X }
X b = (x1*x1)*(y2-y3) - (x2*x2)*(y1-y3) + (x3*x3)*(y1-y2);
X xm = -b/(2.0*a);
X *tmincp = (xm - mjd)*24.0;
X return (fabs (*tmincp) < tmlimit ? -1 : 0);
X}
X
X/* use secant method to solve for time when expression passes through 0.
X */
Xstatic
Xsrch_solve0(mjd, v, tmincp)
Xdouble mjd;
Xdouble v;
Xdouble *tmincp;
X{
X static double x0, x1; /* x(n-1) and x(n) */
X static double y0, y1; /* y(n-1) and y(n) */
X double x2; /* x(n+1) */
X double df; /* y(n) - y(n-1) */
X
X switch (srch_tmscalled) {
X case 0: x0 = mjd; y0 = v; return(0);
X case 1: x1 = mjd; y1 = v; break;
X default: x0 = x1; y0 = y1; x1 = mjd; y1 = v; break;
X }
X
X df = y1 - y0;
X if (fabs(df) < 1e-10) {
X /* near-0 zero denominator, ie, curve is pretty flat here,
X * so assume we are done enough.
X * signal this by forcing a 0 tminc.
X */
X *tmincp = 0.0;
X return (-1);
X }
X x2 = x1 - y1*(x1-x0)/df;
X *tmincp = (x2 - mjd)*24.0;
X return (fabs (*tmincp) < tmlimit ? -1 : 0);
X}
X
X/* binary search for time when expression changes from its initial state.
X * if the change is outside the initial tminc range, then keep searching in that
X * direction by tminc first before starting to divide down.
X */
Xstatic
Xsrch_binary(mjd, v, tmincp)
Xdouble mjd;
Xdouble v;
Xdouble *tmincp;
X{
X static double lb, ub; /* lower and upper bound */
X static int initial_state;
X int this_state = v >= 0.5;
X
X#define FLUNDEF -9e10
X
X if (srch_tmscalled == 0) {
X if (*tmincp >= 0.0) {
X /* going forwards in time so first mjd is lb and no ub yet */
X lb = mjd;
X ub = FLUNDEF;
X } else {
X /* going backwards in time so first mjd is ub and no lb yet */
X ub = mjd;
X lb = FLUNDEF;
X }
X initial_state = this_state;
X return (0);
X }
X
X if (ub != FLUNDEF && lb != FLUNDEF) {
X if (this_state == initial_state)
X lb = mjd;
X else
X ub = mjd;
X *tmincp = ((lb + ub)/2.0 - mjd)*24.0;
X#ifdef TRACEBIN
X { char buf[NC];
X sprintf (buf, "lb=%g ub=%g tminc=%g mjd=%g is=%d ts=%d",
X lb, ub, *tmincp, mjd, initial_state, this_state);
X f_msg (buf);
X }
X#endif
X /* signal to stop if asking for time change less than TMLIMIT */
X return (fabs (*tmincp) < tmlimit ? -1 : 0);
X } else if (this_state != initial_state) {
X /* gone past; turn around half way */
X if (*tmincp >= 0.0)
X ub = mjd;
X else
X lb = mjd;
X *tmincp /= -2.0;
X return (0);
X } else {
X /* just keep going, looking for first state change but we keep
X * learning the lower (or upper, if going backwards) bound.
X */
X if (*tmincp >= 0.0)
X lb = mjd;
X else
X ub = mjd;
X return (0);
X }
X}
EOFxEOF
len=`wc -c < srch.c`
if expr $len != 7803 > /dev/null
then echo Length of srch.c is $len but it should be 7803.
fi
echo x sun.c
sed -e 's/^X//' << 'EOFxEOF' > sun.c
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given the modified JD, mjd, return the true geocentric ecliptic longitude
X * of the sun for the mean equinox of the date, *lsn, in radians, and the
X * sun-earth distance, *rsn, in AU. (the true ecliptic latitude is never more
X * than 1.2 arc seconds and so may be taken to be a constant 0.)
X * if the APPARENT ecliptic longitude is required, correct the longitude for
X * nutation to the true equinox of date and for aberration (light travel time,
X * approximately -9.27e7/186000/(3600*24*365)*2*pi = -9.93e-5 radians).
X */
Xsunpos (mjd, lsn, rsn)
Xdouble mjd;
Xdouble *lsn, *rsn;
X{
X double t, t2;
X double ls, ms; /* mean longitude and mean anomoay */
X double s, nu, ea; /* eccentricity, true anomaly, eccentric anomaly */
X double a, b, a1, b1, c1, d1, e1, h1, dl, dr;
X
X t = mjd/36525.;
X t2 = t*t;
X a = 100.0021359*t;
X b = 360.*(a-(long)a);
X ls = 279.69668+.0003025*t2+b;
X a = 99.99736042000039*t;
X b = 360*(a-(long)a);
X ms = 358.47583-(.00015+.0000033*t)*t2+b;
X s = .016751-.0000418*t-1.26e-07*t2;
X anomaly (degrad(ms), s, &nu, &ea);
X a = 62.55209472000015*t;
X b = 360*(a-(long)a);
X a1 = degrad(153.23+b);
X a = 125.1041894*t;
X b = 360*(a-(long)a);
X b1 = degrad(216.57+b);
X a = 91.56766028*t;
X b = 360*(a-(long)a);
X c1 = degrad(312.69+b);
X a = 1236.853095*t;
X b = 360*(a-(long)a);
X d1 = degrad(350.74-.00144*t2+b);
X e1 = degrad(231.19+20.2*t);
X a = 183.1353208*t;
X b = 360*(a-(long)a);
X h1 = degrad(353.4+b);
X dl = .00134*cos(a1)+.00154*cos(b1)+.002*cos(c1)+.00179*sin(d1)+
X .00178*sin(e1);
X dr = 5.43e-06*sin(a1)+1.575e-05*sin(b1)+1.627e-05*sin(c1)+
X 3.076e-05*cos(d1)+9.27e-06*sin(h1);
X *lsn = nu+degrad(ls-ms+dl);
X *rsn = 1.0000002*(1-s*cos(ea))+dr;
X range (lsn, 2*PI);
X}
EOFxEOF
len=`wc -c < sun.c`
if expr $len != 1760 > /dev/null
then echo Length of sun.c is $len but it should be 1760.
fi
echo x time.c
sed -e 's/^X//' << 'EOFxEOF' > time.c
X/* get the time from the os.
X *
X * here are two methods I was able to verify; pick one for your system and
X * define exactly one of TZA or TZB:
X * TZA works on our ibm-pc/turbo-c and at&t systems,
X * TZB works on our 4.2 BSD vax.
X *
X * I'm told that on Sun OS 4.0.3 (BSD 4.3?) and Apollo SR 10.1 TZB works if
X * you use <sys/time.h> in place of <time.h>.
X */
X#define TZA
X
X#include <stdio.h>
X#include <time.h>
X
X#include "astro.h"
X#include "circum.h"
X
Xextern char *strncpy();
Xextern long time();
X
Xstatic long c0;
Xstatic double mjd0;
X
X/* save current mjd and corresponding system clock for use by inc_mjd().
X * this establishes the base correspondence between the mjd and system clock.
X */
Xset_t0 (np)
XNow *np;
X{
X mjd0 = mjd;
X time (&c0);
X}
X
X/* fill in n_mjd/tz/tznm from system clock.
X */
Xtime_fromsys (np)
XNow *np;
X{
X extern struct tm *gmtime(), *localtime();
X struct tm *tp;
X long c;
X double day, hr;
X
X time (&c);
X
X tp = localtime (&c);
X settzstuff (tp->tm_isdst ? 1 : 0, np);
X
X /* N.B. gmtime() can return pointer to same area as localtime(), ie,
X * reuse the same static area, so be sure to call it after we are
X * through with local time info.
X */
X tp = gmtime (&c);
X cal_mjd (tp->tm_mon+1, (double)tp->tm_mday, tp->tm_year+1900, &day);
X sex_dec (tp->tm_hour, tp->tm_min, tp->tm_sec, &hr);
X mjd = day + hr/24.0;
X}
X
X/* given whether dst is now in effect (must be strictly 0 or 1), fill in
X * tzname and tz within np.
X */
Xstatic
Xsettzstuff (dst, np)
Xint dst;
XNow *np;
X{
X#ifdef TZA
X extern long timezone;
X extern char *tzname[2];
X
X tzset();
X tz = timezone/3600;
X if (dst)
X tz -= 1.0;
X strncpy (tznm, tzname[dst], sizeof(tznm)-1);
X#endif
X#ifdef TZB
X extern char *timezone();
X struct timeval timev;
X struct timezone timez;
X
X gettimeofday (&timev, &timez);
X tz = timez.tz_minuteswest/60;
X if (dst)
X tz -= 1.0;
X strncpy (tznm, timezone(timez.tz_minuteswest, dst), sizeof(tznm)-1);
X#endif
X tznm[sizeof(tznm)-1] = '\0'; /* insure string is terminated */
X}
X
Xinc_mjd (np, inc)
XNow *np;
Xdouble inc;
X{
X if (inc == RTC) {
X long c;
X time (&c);
X mjd = mjd0 + (c - c0)/SPD;
X } else
X mjd += inc/24.0;
X
X /* round to nearest whole second.
X * without this, you can get fractional days so close to .5 but
X * not quite there that mjd_hr() can return 24.0
X */
X rnd_second (&mjd);
X}
EOFxEOF
len=`wc -c < time.c`
if expr $len != 2295 > /dev/null
then echo Length of time.c is $len but it should be 2295.
fi
echo x utc_gst.c
sed -e 's/^X//' << 'EOFxEOF' > utc_gst.c
X#include "astro.h"
X
X/* given a modified julian date, mjd, and a universally coordinated time, utc,
X * return greenwich mean siderial time, *gst.
X */
Xutc_gst (mjd, utc, gst)
Xdouble mjd;
Xdouble utc;
Xdouble *gst;
X{
X double tnaught();
X static double lastmjd;
X static double t0;
X
X if (mjd != lastmjd) {
X t0 = tnaught (mjd);
X lastmjd = mjd;
X }
X *gst = (1.0/SIDRATE)*utc + t0;
X range (gst, 24.0);
X}
X
X/* given a modified julian date, mjd, and a greenwich mean siderial time, gst,
X * return universally coordinated time, *utc.
X */
Xgst_utc (mjd, gst, utc)
Xdouble mjd;
Xdouble gst;
Xdouble *utc;
X{
X double tnaught();
X static double lastmjd;
X static double t0;
X
X if (mjd != lastmjd) {
X t0 = tnaught (mjd);
X range (&t0, 24.0);
X lastmjd = mjd;
X }
X *utc = gst - t0;
X range (utc, 24.0);
X *utc *= SIDRATE;
X}
X
Xstatic double
Xtnaught (mjd)
Xdouble mjd; /* julian days since 1900 jan 0.5 */
X{
X double dmjd;
X int m, y;
X double d;
X double t, t0;
X
X mjd_cal (mjd, &m, &d, &y);
X cal_mjd (1, 0., y, &dmjd);
X t = dmjd/36525;
X t0 = 6.57098e-2 * (mjd - dmjd) -
X (24 - (6.6460656 + (5.1262e-2 + (t * 2.581e-5))*t) -
X (2400 * (t - (((double)y - 1900)/100))));
X return (t0);
X}
EOFxEOF
len=`wc -c < utc_gst.c`
if expr $len != 1171 > /dev/null
then echo Length of utc_gst.c is $len but it should be 1171.
fi
echo x version.c
sed -e 's/^X//' << 'EOFxEOF' > version.c
X/* N.B. please increment version and date and note each change. */
X
X#include "screen.h"
X
Xstatic char vmsg[] = "Version 4.12 February 20, 1990";
X
X/*
X * 4.12 1/12/90 lay framework for orbital element support for object x.
X * 1/15 fix slight bug related to nutation.
X * plot fields in the same order they were selected.
X * 1/18 add copywrite notice in main.c.
X * note <sys/time.h> in time.c for BSD 4.3.
X * 1/20 work on fixed and parabolic orbital elements.
X * 1/25 work on elliptical orbital elements.
X * 2/1 work on objx's magnitude models.
X * add confirmation question before quitting.
X * 2/6 add d,o,z special speed move chars.
X * 2/8 watch: add LST to night sky and maintain RTC time back in main.
X * 2/12 fix bug in timezone related to daytime flag.
X * add w (week) watch advance key code.
X * add cautionary note about no string[s].h to Readme
X * 2/15 allow for precession moving dec outside range -90..90.
X * 2/19 fix bug that wiggled cursor during plotting in rise/set menu.
X * 2/20 fix bug preventing DAWN/DUSK/LON from being used in search func.
X * 4.11 12/29 add PC_GRAPHICS option in mainmenu.c (no effect on unix version)
X * 1/3/90 fix twilight error when sun never gets as low as -18 degs.
X * 1/4/90 always find alt/az from eod ra/dec, not from precessed values.
X * 1/9/90 lastmjd in plans.c was an int: prevented needless recalcs.
X * 4.10 12/6/89 fix transit times of circumpolar objects that don't rise.
X * fix plotting search function when searching is not on.
X * 12/12 fix Objx rise/set bug.
X * 12/21 don't erase last watch positions until computed all new ones.
X * 12/23 added USE_BIOSCALLS to io.c: Doug McDonald's BIOS calls
X * 12/27 allow dates to be entered as decimal years (for help with plots)
X * 12/27 remove literal ESC chars in strings in io.c.
X * 4.9 11/28/89 provide two forms of non-blocking reads for unix in io.c
X * 11/30/89 take out superfluous ESC testing in chk_arrow().
X * guard better against bogus chars in sel_fld().
X * use %lf in scanf's.
X * command-line arg PROPTS+ adds to settings from config file.
X * change (int) casts in moduloes to (long) for 16bit int systems.
X * 4.8 10/28/89 use doubles everywhere
X * 10/31/89 add direct planet row selection codes.
X * 11/2/89 improve compiler's fieldname parser.
X * 11/3/89 switch from ESC to q for "go on" (CBREAK ESC not very portable)
X * 11/6/89 allow plotting the search function too.
X * 11/8/89 suppress screen updates while plotting and nstep > 1.
X * 11/9/89 fix bug prohibiting plotting venus' sdist and objx's transit.
X * 11/9/89 add option to plot in polar coords.
X * 11/12/89 fix bug related to updating timezone name when it is changed.
X * 11/21/89 fix bug in when to print info about object-x
X * 11/21/89 increase MAXPLTLINES to 10 (to ease plotting all planet seps)
X * 11/22/89 allow setting fields from command line too.
X * 4.7 10/13/89 start adding general searching feature. start with flogging.
X * 10/17/89 add compiler, first menu ideas, get binary srch working.
X * 10/18/89 add parabolic-extrema and secant-0 solvers.
X * 10/23/89 finish up new idea of one-line control and set-up "popup" menus.
X * 4.6 10/29/89 improve transit circumstances by iterating as with rise/set.
X * allow changing lst.
X * show Updating message at better times.
X * avoid overstrikes while watching and add trails option.
X * allow for Turbo-C 2.0 printf bug using %?.0f".
X * 4.5 9/24/89 add third table of all mutual planet angular distances.
X * 4.4 9/21/89 add second planet table with rise/set times.
X * all rise/set times may now use standard or adaptive horizons.
X * 4.3 9/6/89 NM/FM calendar overstikes now use local time (was ut).
X * display elongation of object x.
X * better handling of typo when asking for refraction model.
X * 4.2 7/24/89 specify 7 digits to plot file (not just default to 6)
X * 4.1 7/18/89 use buffered output and fflush in read_char().
X * 4.0 7/8/89 add simple sky and solarsystem plotting (and rearrange fields)
X * change mars' .cfg mnemonic from a to m.
X * clean up nstep/NEW CIR handling
X * quit adding our own .cfg suffixes, but...
X * add looking for $HOME/.ephemrc (Ronald Florence)
X * drop -b
X * no longer support SITE
X * 3.17 6/15/89 misspelt temperature prompt; sun -/= bug. (Mike McCants)
X * change sun() to sunpos() for sake of Sun Microsystems.
X * 3.16 6/9/89 allow JD to be set and plotted.
X * c_clear (LATTIC_C) should use J not j (Alex Pruss)
X * support SIGINT (Ronald Florence)
X * 3.15 6/8/89 forget SIMPLETZ: now TZA and TZB.
X * 3.14 6/6/89 add back borders but allow turning off via -b
X * 3.13 5/26/89 fix Day/Nite picking loc bug.
X * 3.12 5/25/89 add SIMPLETZ option to time.c for systems w/o tzset()
X * files; couldn't plot dayln or niteln.
X * 3.11 5/16/89 local time prompt said utc; add NiteLn; check for bad plot
X * 3.10 4/27/89 allow caps for moving cursor around too
X * 3.9 4/5/89 discard leading termcap delay digits, for now
X * 3.8 3/2/89 shorten displayed precision, add heliocentric lat/long
X * 3.7 2/13/89 change to ^d to quit program.
X * 3.6 2/7/89 try adding .cfg suffix if can't find config file
X * 3.5 2/2/89 sunrise/set times based on actual sun size and pressure/temp
X * 3.4 1/22/89 calendar and all rise/set times based on local date, not utc
X * 3.3 1/6/89 add z to plot files (still don't plot it however)
X * 3.2 1/3/89 if set date/time then time does not inc first step
X * 3.1 1/1/89 user def. graph labels; nstep/stpsz control; genuine c_eol
X * 3.0 12/31/88 add graphing; add version to credits.
X * 2.7 12/30/88 add version to credits.
X * 2.6 12/28/88 twilight defined as 18 rather than 15 degrees below horizon
X * 2.5 12/26/88 remove trace capability; add screen shadowing: ^l.
X * 2.4 10/31/88 add credits banner, -s turns it off; savings time fix.
X * 2.3 9/23/88 exchange Altitude/Elevation titles (no code changes)
X * 2.2 9/20/88 more caution in aaha_aux() guarding acos() arg range
X * 2.1 9/14/88 moon phase always >= 0 to match planets convention
X * 2.0 9/13/88 add version ^v option
X */
X
Xversion()
X{
X f_msg (vmsg);
X}
X
Xstatic char *cre[] = {
X"Ephem - an interactive astronomical ephemeris program",
Xvmsg,
X"",
X"Copyright (c) 1990 by Elwood Charles Downey",
X"",
X"Permission is granted to make and distribute copies of this program",
X"free of charge, provided the copyright notice and this permission",
X"notice are preserved on all copies. All other rights reserved.",
X"",
X"Many formulas and tables are based, with permission, on material found in",
X"\"Astronomy with your Personal Computer\"",
X"by Dr. Peter Duffett-Smith, Cambridge University Press, (c) 1985",
X"",
X"type any key to continue..."
X};
Xcredits()
X{
X int r = 6; /* first row of credits message */
X int l;
X
X c_erase();
X for (l = 0; l < sizeof(cre)/sizeof(cre[0]); l++)
X f_string (r++, (NC - strlen(cre[l]))/2, cre[l]);
X (void) read_char(); /* wait for any char to continue */
X}
EOFxEOF
len=`wc -c < version.c`
if expr $len != 6979 > /dev/null
then echo Length of version.c is $len but it should be 6979.
fi
echo x watch.c
sed -e 's/^X//' << 'EOFxEOF' > watch.c
X/* these functions allow you to watch the sky or the solar system via a
X * simple character-graphics representation on the screen.
X * the interaction starts by using the current time. then control with
X * END returns to table form; or
X * RETURN advances time by one StpSz; or
X * h advances once by 1 hour; or
X * d advances once by 24 hours (1 day); or
X * w advances once by 7 days (1 week); or
X * any other key free runs by StpSz until any key is hit.
X */
X
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X#include "circum.h"
X#include "screen.h"
X
X#define TROW (R_PROMPT+1) /* time/date row */
X#define TCOL C_PROMPT /* " col */
X#define SSZCOL 1 /* column to show solar system z coords */
X
X#define SKYACC 3600. /* desired sky plot accuracy, in arc seconds */
X#define SSACC 3600. /* desired solar system plot accuracy, in arc secs */
X
X/* macros to convert row(col) in range 1..NR(1..NC) to fraction in range 0..1 */
X#define r2fr(r) (((r)-1)/(NR-1))
X#define c2fc(c) (((c)-1)/(NC-1))
X#define fr2r(fr) ((fr)*(NR-1)+1)
X#define fc2c(fc) ((fc)*(NC-1)+1)
X
X/* single-character tag for each body.
X * order must match the #defines in astro.h and screen.h additions.
X */
Xstatic char body_tags[] = "evmjsunpSMx";
X
X/* multiple and single loop prompts */
Xstatic char frprompt[] = "Running... press any key to stop.";
Xstatic char qprompt[] =
X"q to quit, RETURN/h/d/w to step by StpSz/hr/day/wk, or any other to freerun";
X
X/* used to locate, record and then erase last plotted chars */
Xtypedef struct {
X double l_fr, l_fc; /* 2d coords as 0..1 */
X int l_r, l_c; /* screen 2d coords */
X char l_tag; /* char to use to print on screen */
X} LastDraw;
X
Xstatic int trails; /* !0 if want to leave trails */
X
Xwatch (np, tminc, wbodies)
XNow *np; /* time now and on each step */
Xdouble tminc; /* hrs to increment time by each step */
Xint wbodies; /* each bit is !=0 if want that body */
X{
X static char *flds[3] = {
X "Night sky", "Solar system"
X };
X int fn = 0;
X
X ask:
X flds[2] = trails ? "Leave trails" : "No trails";
X switch (popup (flds, fn, 3)) {
X case 0: watch_sky (np, tminc, wbodies); break;
X case 1: watch_solarsystem (np, tminc, wbodies); break;
X case 2: trails ^= 1; fn = 2; goto ask; /* toggle trails and repeat */
X default: break;
X }
X}
X
X/* full night sky view.
X * north is at left and right of screen south at center.
X * 0 elevation is at bottom of screen, zenith at the top.
X * redraw screen when done.
X */
Xstatic
Xwatch_sky (np, tminc, wbodies)
XNow *np; /* time now and on each step */
Xdouble tminc; /* hrs to increment time by each step */
Xint wbodies; /* each bit is !=0 if want */
X{
X static char title[] = "Sky at";
X int tcol = sizeof(title)+1;
X double tminc0 = tminc; /* remember the original */
X /* two draw buffers so we can leave old up while calc new then
X * erase and draw in one quick operation. always calc new in newp
X * buffer and erase previous from lastp. buffers alternate roles.
X */
X LastDraw ld0[NOBJ], ld1[NOBJ], *lp, *lastp = ld0, *newp = ld1;
X int nlast = 0, nnew;
X int once = 1;
X double lmjd, tmp;
X Sky s;
X int p;
X
X c_erase();
X f_string (TROW, TCOL, title);
X
X while (1) {
X if (once)
X print_updating();
X
X /* calculate desired stuff into newp[] */
X nnew = 0;
X for (p = nxtbody(-1); p != -1; p = nxtbody(p))
X if (wbodies & (1<<p)) {
X (void) body_cir (p, SKYACC, np, &s);
X if (s.s_alt >= 0) {
X LastDraw *lnp = newp + nnew;
X lnp->l_fr = 1.0 - s.s_alt/(PI/2);
X lnp->l_fc = s.s_az/(2*PI);
X lnp->l_tag = body_tags[p];
X nnew++;
X }
X }
X set_screencoords (newp, nnew);
X
X /* unless we want trails,
X * erase any previous tags (in same order as written) from lastp[].
X */
X if (!trails)
X for (lp = lastp; --nlast >= 0; lp++)
X f_char (lp->l_r, lp->l_c, ' ');
X
X /* print LOCAL time and date we will be using */
X lmjd = mjd - tz/24.0;
X f_time (TROW, tcol, mjd_hr(lmjd));
X f_date (TROW, tcol+9, mjd_day(lmjd));
X now_lst (np, &tmp);
X f_time (TROW, tcol+20, tmp);
X printf ("LST");
X
X /* now draw new stuff from newp[] */
X for (lp = newp; lp < newp + nnew; lp++)
X f_char (lp->l_r, lp->l_c, lp->l_tag);
X
X /* swap new and last roles and save new count */
X if (newp == ld0)
X newp = ld1, lastp = ld0;
X else
X newp = ld0, lastp = ld1;
X nlast = nnew;
X
X if (once || (chk_char()==0 && read_char()!=0)) {
X if (readwcmd (tminc0, &tminc, &once) < 0)
X break;
X }
X
X /* advance time */
X inc_mjd (np, tminc);
X }
X
X redraw_screen(2);
X}
X
X/* solar system view, "down from the top", first point of aries to the right.
X * always include earth.
X * redraw screen when done.
X */
Xstatic
Xwatch_solarsystem (np, tminc, wbodies)
XNow *np; /* time now and on each step */
Xdouble tminc; /* hrs to increment time by each step */
Xint wbodies;
X{
X /* max au of each planet from sun; in astro.h #defines order */
X static double auscale[] = {.38, .75, 1.7, 5.2, 11., 20., 31., 39.};
X static char title[] = "Solar System at";
X int tcol = sizeof(title)+1;
X double tminc0 = tminc; /* remember the original */
X /* two draw buffers so we can leave old up while calc new then
X * erase and draw in one quick operation. always calc new in newp
X * buffer and erase previous from lastp. buffers alternate roles.
X */
X LastDraw ld0[2*NOBJ], ld1[2*NOBJ], *lp, *lastp = ld0, *newp = ld1;
X int nlast = 0, nnew;
X int once = 1;
X double lmjd;
X double scale;
X Sky s;
X int p;
X
X /* set screen scale: largest au we will have to plot.
X * never make it less than 1 au since we always show earth.
X */
X scale = 1.0;
X for (p = MARS; p <= PLUTO; p++)
X if ((wbodies & (1<<p)) && auscale[p] > scale)
X scale = auscale[p];
X
X c_erase();
X f_string (TROW, TCOL, title);
X
X while (1) {
X if (once)
X print_updating();
X
X /* calculate desired stuff into newp[].
X * fake a sun at center and add earth first.
X * (we get earth's loc when ask for sun)
X */
X nnew = 0;
X set_ss (newp+nnew, 0.0, 0.0, 0.0, 'S');
X nnew += 2;
X (void) body_cir (SUN, SSACC, np, &s);
X set_ss (newp+nnew, s.s_edist/scale, s.s_hlong, 0.0, 'E');
X nnew += 2;
X for (p = MERCURY; p <= PLUTO; p++)
X if (wbodies & (1<<p)) {
X (void) body_cir (p, SSACC, np, &s);
X set_ss (newp+nnew, s.s_sdist/scale, s.s_hlong, s.s_hlat,
X body_tags[p]);
X nnew += 2;
X }
X if (wbodies & (1<<OBJX)) {
X (void) body_cir (OBJX, SSACC, np, &s);
X if (s.s_hlong != NOHELIO && s.s_sdist <= scale) {
X set_ss (newp+nnew, s.s_sdist/scale, s.s_hlong, s.s_hlat,
X body_tags[OBJX]);
X nnew += 2;
X }
X }
X set_screencoords (newp, nnew);
X
X /* unless we want trails,
X * erase any previous tags (in same order as written) from lastp[].
X */
X if (!trails)
X for (lp = lastp; --nlast >= 0; lp++)
X f_char (lp->l_r, lp->l_c, ' ');
X
X /* print LOCAL time and date we will be using */
X lmjd = mjd - tz/24.0;
X f_time (TROW, tcol, mjd_hr(lmjd));
X f_date (TROW, tcol+9, mjd_day(lmjd));
X
X /* now draw new stuff from newp[] */
X for (lp = newp; lp < newp + nnew; lp++)
X f_char (lp->l_r, lp->l_c, lp->l_tag);
X
X /* swap new and last roles and save new count */
X if (newp == ld0)
X newp = ld1, lastp = ld0;
X else
X newp = ld0, lastp = ld1;
X nlast = nnew;
X
X if (once || (chk_char()==0 && read_char()!=0)) {
X if (readwcmd (tminc0, &tminc, &once) < 0)
X break;
X }
X
X /* advance time */
X inc_mjd (np, tminc);
X }
X
X redraw_screen(2);
X}
X
X/* fill in two LastDraw solar system entries,
X * one for the x/y display, one for the z.
X */
Xstatic
Xset_ss (lp, dist, lg, lt, tag)
XLastDraw *lp;
Xdouble dist, lg, lt; /* scaled heliocentric distance, longitude and lat */
Xchar tag;
X{
X lp->l_fr = 0.5 - dist*sin(lg)*0.5;
X lp->l_fc = 0.5 + dist*cos(lg)*0.5/ASPECT;
X lp->l_tag = tag;
X lp++;
X /* row is to show course helio altitude but since we resolve collisions
X * by adjusting columns we can get more detail by smaller variations
X * within one column.
X */
X lp->l_fr = 0.5 - dist*sin(lt)*0.5;
X lp->l_fc = c2fc(SSZCOL) + (1 - lp->l_fr)/NC;
X lp->l_tag = tag;
X}
X
X/* given a list of LastDraw structs with their l_{fr,fc} filled in,
X * fill in their l_{r,c}.
X * TODO: better collision avoidance.
X */
Xstatic
Xset_screencoords (lp, np)
XLastDraw lp[];
Xint np;
X{
X LastDraw *lpi; /* the current basis for comparison */
X LastDraw *lpj; /* the sweep over other existing cells */
X int i; /* index of the current basis cell, lpi */
X int j; /* index of sweep cell, lpj */
X int n; /* total cells placed so far (ie, # to check) */
X
X /* idea is to place each new item onto the screen.
X * after each placement, look for collisions.
X * if find a colliding pair, move the one with the great l_fc to
X * the right one cell, then rescan for more collisions.
X * this will yield a result that is sorted by columns by l_fc.
X * TODO: don't just move to the right, try up too for true 2d adjusts.
X */
X for (n = 0; n < np; n++) {
X lpi = lp + n;
X i = n;
X lpi->l_r = fr2r(lpi->l_fr);
X lpi->l_c = fc2c(lpi->l_fc);
X chk:
X for (j = 0; j < n; j++) {
X lpj = lp + j;
X if (i!=j && lpi->l_r == lpj->l_r && lpi->l_c == lpj->l_c) {
X if (lpj->l_fc > lpi->l_fc) {
X /* move lpj and use it as basis for checks now */
X lpi = lpj;
X i = j;
X }
X if (++lpi->l_c > NC)
X lpi->l_c = 1;
X goto chk;
X }
X }
X }
X}
X
X/* see what the op wants to do now and update prompt/times accordingly.
X * return -1 if we are finished, else 0.
X */
Xstatic int
Xreadwcmd (tminc0, tminc, once)
Xdouble tminc0;
Xdouble *tminc;
Xint *once;
X{
X f_prompt (qprompt);
X
X switch (read_char()) {
X case END: /* back to table */
X return (-1);
X case '\r': /* one StpSz step */
X *tminc = tminc0;
X *once = 1;
X break;
X case 'h': /* one 1-hour step */
X *tminc = 1.0;
X *once = 1;
X break;
X case 'd': /* one 24-hr step */
X *tminc = 24.0;
X *once = 1;
X break;
X case 'w': /* 7 day step */
X *tminc = 7*24.0;
X *once = 1;
X break;
X default: /* free-run */
X *once = 0;
X f_prompt (frprompt);
X }
X return (0);
X}
EOFxEOF
len=`wc -c < watch.c`
if expr $len != 10024 > /dev/null
then echo Length of watch.c is $len but it should be 10024.
fi