downey@cs.umn.edu@dimed1.UUCP (08/31/90)
Posting-number: Volume 14, Issue 77
Submitted-by: downey@cs.umn.edu@dimed1.UUCP
Archive-name: ephem-4.21/part02
#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 2 (of 6)."
# Contents: Readme formats.c listing.c mainmenu.c moon.c riset_c.c
# screen.h srch.c
# Wrapped by allbery@uunet on Thu Aug 30 20:46:31 1990
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'Readme' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'Readme'\"
else
echo shar: Extracting \"'Readme'\" \(6592 characters\)
sed "s/^X//" >'Readme' <<'END_OF_FILE'
XGeneral Notes for version 4.20, August 17, 1990:
X
X1) Start by reading the generic-printer-ready manual in Man.txt. If you have
X source, read on in this file for building hints.
X
X2) Beware that ephem does not gracefully handle the fact that there is no year
X 0 in the Christian calendar. It is best to avoid entering a year of 0 for
X now.
X
XRecent Change History:
X
XChanges from 4.19 to 4.20:
X add g/k and H/G magnitude model options for elliptical objects.
X put moon's geocentric long/lat under Hlong/Hlat columns.
X allow entering negative decimal years.
X init all static mjd storage to unlike times.
X add USE_TERMIO option to io.c.
X
XChanges from 4.18 to 4.19:
X add listing feature, with 'L' hot-key.
X add title for plot file (as well as listing file).
X add some (void) casts for lint sake.
X
XChanges from 4.17 to 4.18:
X fix parabolic comet bug in objx.c (bad lam computation).
X
XChanges from 4.16 to 4.17:
X add 'c' short cut to Menu field.
X display full Dec precision for fixed objx setup.
X increase pluto auscale in watch.c, and guard screen boundries.
X add Pause field.
X fruther improve rise/set and dawn/dusk times.
X add MENU={DATA,RISET,SEP} config/arg option.
X
XChanges from 4.15 to 4.16:
X watch popup now allows changing formats without returning.
X add 'w' short cut to watch field.
X improve labeling a bit in Dome display.
X
XChanges from 4.14 to 4.15:
X move setjmp() in main so it catches error in config file too.
X maintain name of objx/y and generally clean up objx.c.
X fix bug in circum.c related to phase of fixed objects.
X add "Sky dome" watch format and improve labels.
X remember last watch, plot and search popup selections.
X
XChanges from 4.13 to 4.14:
X fix bug watching object y by adding y to body_tags[] in watch.c.
X add ! support (#ifdef BANG in sel_fld()).
X add support for VMS via #ifdef VMS.
X switch to EPHEMCFG environ variable (HOME no longer used).
X fix phase so it works for objects out of the ecliptic.
X
XIf you have source, here is how you build ephem:
X
X1) io.c:
X define UNIX, VMS or TURBO_C in io.c depending on your system. Note that
X TURBO_C also seems to work ok with Microsoft and Lattice C too but these
X have not been tested recently. Note also that the VMS C compiler defines VMS
X automatically so you don't really have to #define it.
X
X Also in io.c and if you use UNIX, you have two choices of methods for doing
X non-blocking reads depending on whether you define USE_NDELAY; use whichever
X works for you. You can also choose between sgtty.h and termio.h; again,
X use whichever works on your system; for the latter, define USE_TERMIO.
X
X MSDOS users can do cursor control with direct BIOS calls (the default),
X or with ANSI.SYS by defining USE_ANSISYS in io.c.
X
X2) time.c:
X Select from two methods of dealing with time from the operating system
X with the TZA/TZB defines in time.c. If you get link undefines related to
X time functions try changing to the other form.
X
X On Sun systems running OS 4.0.3 (or BSD 4.3) or Apollo SR 10.1 use TZB and
X change <time.h> to <sys/time.h>.
X
X For VMS, since it does not support time zone info, do NOT #define EITHER
X of TZA or TZB. This will have the effect of leaving the time zone unchanged
X whenever you set the time via the "Now" option. (This is taken care of
X automatically in time.c by #undef'ining TZA and TZB if VMS is defined.)
X
X3) mainmenu.c:
X if you are compiling for an IBM-PC then try #define PC_GRAPHICS for a
X nicer looking way to draw the screen boundry lines.
X
X4) sel_fld.c:
X if your runtime library supports the system() function (to run a shell
X command) then #define BANG in sel_fld.c, else leave it undefine. When
X defined, this will allow you to jump out of ephem and run any command,
X then resume where you left off.
X
X5) beware that I have not used string.h or strings.h. if your library's
X strlen() and str.*cmp() functions don't return int (such as long), then you
X will have to hand add string.h or your own extern declarations. I have
X included all the necessary declarations for the functions that return
X (char *) such as strcpy(), etc, though.
X
X6) main.c calls sleep() which is not in some IBM-PC runtime C libraries. You
X might kludge up your own call that does a cpu countdown loop. The accuracy
X will only effect the Pause feature, not ephem's actual time mechanisms.
X
X7) Ephem can now be built by simply compiling all the .c files and linking them
X all together. On Unix systems, you must also link with the termcap library
X (-ltermcap) and possibly the auxiliary math library (-lm) if your default C
X library does not include all the required transcendental functions. At the
X end of this file I have included a VMS build script for those building
X ephem on a that system.
X
XThe following files are pretty much just pure transliterations from BASIC
Xinto C from machine-readable copies of the programs in Duffett-Smith's book.
XThey have nothing to do with the rest of ephem so they may be used for
Xcompletely different applications if so desired.
X
X aa_hadec.c anomaly.c astro.h cal_mjd.c comet.c eq_ecl.c moon.c moonnf.c
X nutation.c obliq.c parallax.c pelement.c plans.c precess.c reduce.c
X refract.c sex_dec.c sun.c utc_gst.c
X
X$!========================================================================
X$!
X$! Name : BUILD.COM
X$!
X$! Purpose : compile and link ephem under VMS
X$!
X$! Arguments : P1/P2 = DEBUG: compile with DEBUG info
X$! P1/P2 = LINK : link only
X$!
X$! Created 23-MAR-1990 Karsten Spang
X$!
X$!========================================================================
X$ ON ERROR THEN $ GOTO EXIT
X$ ON CONTROL_Y THEN $ GOTO EXIT
X$ IF P1.EQS."DEBUG" .OR. P2.EQS."DEBUG"
X$ THEN
X$ CC:=CC/DEBUG/NOOPTIMIZE/NOLIST
X$ LINK:=LINK/DEBUG/NOMAP
X$ ELSE
X$ CC:=CC/NOLIST
X$ LINK:=LINK/NOMAP
X$ ENDIF
X$ FILES="MAIN,AA_HADEC,ALTMENUS,ANOMALY,CAL_MJD,CIRCUM,COMET,COMPILER,"+ -
X "EQ_ECL,FLOG,FORMATS,IO,MAINMENU,MOON,MOONNF,NUTATION,"+ -
X "OBJX,OBLIQ,PARALLAX,PELEMENT,PLANS,PLOT,POPUP,PRECESS,REDUCE,"+ -
X "REFRACT,RISET,RISET_C,SEL_FLD,SEX_DEC,SRCH,SUN,TIME,UTC_GST,"+ -
X "VERSION,WATCH"
X$ IF P1.EQS."LINK" .OR. P2.EQS."LINK" THEN GOTO LINK
X$ FILE_NUM=0;
X$COMPILE_LOOP:
X$ FILE=F$ELEMENT(FILE_NUM,",",FILES)
X$ IF FILE.EQS."," THEN GOTO COMPILE_END
X$ CC 'FILE'
X$ FILE_NUM=FILE_NUM+1
X$ GOTO COMPILE_LOOP
X$COMPILE_END:
X$LINK:
X$ LINK/EXE=EPHEM 'FILES',SYS$INPUT/OPT
XSYS$LIBRARY:VAXCRTL/SHARE
X$EXIT:
X$ EXIT
END_OF_FILE
if test 6592 -ne `wc -c <'Readme'`; then
echo shar: \"'Readme'\" unpacked with wrong size!
fi
# end of 'Readme'
fi
if test -f 'formats.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'formats.c'\"
else
echo shar: Extracting \"'formats.c'\" \(7453 characters\)
sed "s/^X//" >'formats.c' <<'END_OF_FILE'
X/* basic formating routines.
X * all the screen oriented printing should go through here.
X */
X
X#include <stdio.h>
X#include <math.h>
X#ifdef VMS
X#include <stdlib.h>
X#endif
X#include "astro.h"
X#include "screen.h"
X
Xextern char *strcpy();
X
X/* suppress screen io if this is true, but always flog stuff.
X */
Xstatic int f_scrnoff;
Xf_on ()
X{
X f_scrnoff = 0;
X}
Xf_off ()
X{
X f_scrnoff = 1;
X}
X
X/* draw n blanks at the given cursor position. */
Xf_blanks (r, c, n)
Xint r, c, n;
X{
X if (f_scrnoff)
X return;
X c_pos (r, c);
X while (--n >= 0)
X putchar (' ');
X}
X
X/* print the given value, v, in "sexadecimal" format at [r,c]
X * ie, in the form A:m.P, where A is a digits wide, P is p digits.
X * if p == 0, then no decimal point either.
X */
Xf_sexad (r, c, a, p, mod, v)
Xint r, c;
Xint a, p; /* left space, min precision */
Xint mod; /* don't let whole portion get this big */
Xdouble v;
X{
X char astr[32], str[32];
X long dec;
X double frac;
X int visneg;
X double vsav = v;
X
X if (v >= 0.0)
X visneg = 0;
X else {
X if (v <= -0.5/60.0*pow(10.0,-1.0*p)) {
X v = -v;
X visneg = 1;
X } else {
X /* don't show as negative if less than the precision showing */
X v = 0.0;
X visneg = 0;
X }
X }
X
X dec = v;
X frac = (v - dec)*60.0;
X (void) sprintf (str, "59.%.*s5", p, "999999999");
X if (frac >= atof (str)) {
X dec += 1;
X frac = 0.0;
X }
X dec %= mod;
X if (dec == 0 && visneg)
X (void) strcpy (str, "-0");
X else
X (void) sprintf (str, "%ld", visneg ? -dec : dec);
X
X /* would just do this if Turbo-C 2.0 %?.0f" worked:
X * sprintf (astr, "%*s:%0*.*f", a, str, p == 0 ? 2 : p+3, p, frac);
X */
X if (p == 0)
X (void) sprintf (astr, "%*s:%02d", a, str, (int)(frac+0.5));
X else
X (void) sprintf (astr, "%*s:%0*.*f", a, str, p+3, p, frac);
X
X (void) flog_log (r, c, vsav, astr);
X
X f_string (r, c, astr);
X}
X
X/* print the given value, t, in sexagesimal format at [r,c]
X * ie, in the form T:mm:ss, where T is nd digits wide.
X * N.B. we assume nd >= 2.
X */
Xf_sexag (r, c, nd, t)
Xint r, c, nd;
Xdouble t;
X{
X char tstr[32];
X int h, m, s;
X int tisneg;
X
X dec_sex (t, &h, &m, &s, &tisneg);
X if (h == 0 && tisneg)
X (void) sprintf (tstr, "%*s-0:%02d:%02d", nd-2, "", m, s);
X else
X (void) sprintf (tstr, "%*d:%02d:%02d", nd, tisneg ? -h : h, m, s);
X
X (void) flog_log (r, c, t, tstr);
X f_string (r, c, tstr);
X}
X
X/* print angle ra, in radians, in ra hours as hh:mm.m at [r,c]
X * N.B. we assume ra is >= 0.
X */
Xf_ra (r, c, ra)
Xint r, c;
Xdouble ra;
X{
X f_sexad (r, c, 2, 1, 24, radhr(ra));
X}
X
X/* print time, t, as hh:mm:ss */
Xf_time (r, c, t)
Xint r, c;
Xdouble t;
X{
X f_sexag (r, c, 2, t);
X}
X
X/* print time, t, as +/-hh:mm:ss (don't show leading +) */
Xf_signtime (r, c, t)
Xint r, c;
Xdouble t;
X{
X f_sexag (r, c, 3, t);
X}
X
X/* print time, t, as hh:mm */
Xf_mtime (r, c, t)
Xint r, c;
Xdouble t;
X{
X f_sexad (r, c, 2, 0, 24, t);
X}
X
X/* print angle, a, in rads, as degress at [r,c] in form ddd:mm */
Xf_angle(r, c, a)
Xint r, c;
Xdouble a;
X{
X f_sexad (r, c, 3, 0, 360, raddeg(a));
X}
X
X/* print angle, a, in rads, as degress at [r,c] in form dddd:mm:ss */
Xf_gangle(r, c, a)
Xint r, c;
Xdouble a;
X{
X f_sexag (r, c, 4, raddeg(a));
X}
X
X/* print the given modified Julian date, jd, as the starting date at [r,c]
X * in the form mm/dd/yyyy.
X */
Xf_date (r, c, jd)
Xint r, c;
Xdouble jd;
X{
X char dstr[32];
X int m, y;
X double d, tmp;
X
X mjd_cal (jd, &m, &d, &y);
X (void) sprintf (dstr, "%2d/%02d/%04d", m, (int)(d), y);
X
X /* shadow to the plot subsystem as years. */
X mjd_year (jd, &tmp);
X (void) flog_log (r, c, tmp, dstr);
X f_string (r, c, dstr);
X}
X
X/* print the given double as a rounded int, with the given format.
X * this is used to plot full precision, but display far less.
X * N.B. caller beware that we really do expect fmt to refer to an int, not
X * a long for example. also beware of range that implies.
X */
Xf_int (row, col, fmt, f)
Xint row, col;
Xchar fmt[];
Xdouble f;
X{
X char str[80];
X int i;
X
X i = (f < 0) ? (int)(f-0.5) : (int)(f+0.5);
X (void) sprintf (str, fmt, i);
X
X (void) flog_log (row, col, f, str);
X f_string (row, col, str);
X}
X
Xf_char (row, col, c)
Xint row, col;
Xchar c;
X{
X if (f_scrnoff)
X return;
X c_pos (row, col);
X putchar (c);
X}
X
Xf_string (r, c, s)
Xint r, c;
Xchar *s;
X{
X if (f_scrnoff)
X return;
X c_pos (r, c);
X (void) fputs (s, stdout);
X}
X
Xf_double (r, c, fmt, f)
Xint r, c;
Xchar *fmt;
Xdouble f;
X{
X char str[80];
X (void) sprintf (str, fmt, f);
X (void) flog_log (r, c, f, str);
X f_string (r, c, str);
X}
X
X/* print prompt line */
Xf_prompt (p)
Xchar *p;
X{
X c_pos (R_PROMPT, C_PROMPT);
X c_eol ();
X c_pos (R_PROMPT, C_PROMPT);
X (void) fputs (p, stdout);
X}
X
X/* clear from [r,c] to end of line, if we are drawing now. */
Xf_eol (r, c)
Xint r, c;
X{
X if (!f_scrnoff) {
X c_pos (r, c);
X c_eol();
X }
X}
X
X/* print a message and wait for op to hit any key */
Xf_msg (m)
Xchar *m;
X{
X f_prompt (m);
X (void) read_char();
X}
X
X/* crack a line of the form X?X?X into its components,
X * where X is an integer and ? can be any character except '0-9' or '-',
X * such as ':' or '/'.
X * only change those fields that are specified:
X * eg: ::10 only changes *s
X * 10 only changes *d
X * 10:0 changes *d and *m
X * if see '-' anywhere, first non-zero component will be made negative.
X */
Xf_sscansex (bp, d, m, s)
Xchar *bp;
Xint *d, *m, *s;
X{
X char c;
X int *p = d;
X int *nonzp = 0;
X int sawneg = 0;
X int innum = 0;
X
X while (c = *bp++)
X if (c >= '0' && c <= '9') {
X if (!innum) {
X *p = 0;
X innum = 1;
X }
X *p = *p*10 + (c - '0');
X if (*p && !nonzp)
X nonzp = p;
X } else if (c == '-') {
X sawneg = 1;
X } else if (c != ' ') {
X /* advance to next component */
X p = (p == d) ? m : s;
X innum = 0;
X }
X
X if (sawneg && nonzp)
X *nonzp = -*nonzp;
X}
X
X/* crack a floating date string, bp, of the form m/d/y, where d may be a
X * floating point number, into its components.
X * leave any component unspecified unchanged.
X * actually, the slashes may be anything but digits or a decimal point.
X * this is functionally the same as f_sscansex() exept we allow for
X * the day portion to be real, and we don't handle negative numbers.
X * maybe someday we could make a combined one and use it everywhere.
X */
Xf_sscandate (bp, m, d, y)
Xchar *bp;
Xint *m, *y;
Xdouble *d;
X{
X char *bp0, c;
X
X bp0 = bp;
X while ((c = *bp++) && (c >= '0' && c <= '9'))
X continue;
X if (bp > bp0+1)
X *m = atoi (bp0);
X if (c == '\0')
X return;
X bp0 = bp;
X while ((c = *bp++) && (c >= '0' && c <= '9' || c == '.'))
X continue;
X if (bp > bp0+1)
X *d = atof (bp0);
X if (c == '\0')
X return;
X bp0 = bp;
X while (c = *bp++)
X continue;
X if (bp > bp0+1)
X *y = atoi (bp0);
X}
X
X/* just like dec_sex() but makes the first non-zero element negative if
X * x is negative (instead of returning a sign flag).
X */
Xf_dec_sexsign (x, h, m, s)
Xdouble x;
Xint *h, *m, *s;
X{
X int n;
X dec_sex (x, h, m, s, &n);
X if (n) {
X if (*h)
X *h = -*h;
X else if (*m)
X *m = -*m;
X else
X *s = -*s;
X }
X}
X
X/* return 1 if bp looks like a decimal year; else 0.
X * any number greater than 12 or less than 0 is assumed to be a year, or any
X * string with exactly one decimal point, an optional minus sign, and nothing
X * else but digits.
X */
Xdecimal_year (bp)
Xchar *bp;
X{
X char c;
X int ndig = 0, ndp = 0, nneg = 0, nchar = 0;
X double y = atof(bp);
X
X while (c = *bp++) {
X nchar++;
X if (c >= '0' && c <= '9')
X ndig++;
X else if (c == '.')
X ndp++;
X else if (c == '-')
X nneg++;
X }
X
X return (y > 12 || y < 0
X || (ndp == 1 && nneg <= 1 && nchar == ndig+ndp+nneg));
X}
END_OF_FILE
if test 7453 -ne `wc -c <'formats.c'`; then
echo shar: \"'formats.c'\" unpacked with wrong size!
fi
# end of 'formats.c'
fi
if test -f 'listing.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'listing.c'\"
else
echo shar: Extracting \"'listing.c'\" \(7232 characters\)
sed "s/^X//" >'listing.c' <<'END_OF_FILE'
X/* code to support the listing capabilities.
X * idea is to let the operator name a listing file and mark some fields for
X * logging. then after each screen update, the logged fields are written to
X * the listing file in the same manner as they appeared on the screen.
X *
X * format of the listing file is one line per screen update.
X */
X
X#include <stdio.h>
X#include <math.h>
X#include "screen.h"
X
Xextern char *strcpy();
Xextern errno;
Xextern char *sys_errlist[];
X#define errsys (sys_errlist[errno])
X
X#define TRACE(x) {FILE *fp = fopen("trace","a"); fprintf x; fclose(fp);}
X
X#define MAXLSTFLDS 10 /* max number of fields we can track.
X * note we can't store more than NFLOGS fields
X * anyway (see flog.c).
X */
X#define FNLEN (14+1) /* longest filename; plus 1 for \0 */
X
Xstatic char lst_filename[FNLEN] = "ephem.lst"; /* default plot file name */
Xstatic FILE *lst_fp; /* the plot file; == 0 means don't plot */
X
X/* store rcfpack()s for each field to track, in l-to-r order */
Xstatic int lstflds[MAXLSTFLDS];
Xstatic int nlstflds; /* number of lstflds[] in actual use */
X
Xstatic int lstsrchfld; /* set when the Search field is to be listed */
X
X/* picked the Listing label:
X * if on, just turn it off.
X * if off, turn on, define fields or select name of file to list to and do it.
X * TODO: more flexibility, more relevance.
X */
Xlisting_setup()
X{
X if (lst_fp)
X lst_turn_off();
X else {
X static char *chcs[] = {
X "Select fields", "Display a listing file", "Begin listing"
X };
X static int fn; /* start with 0, then remember for next time */
X ask:
X switch (popup(chcs, fn, nlstflds > 0 ? 3 : 2)) {
X case 0: fn = 0; lst_select_fields(); goto ask;
X case 1: fn = 1; lst_file(); goto ask;
X case 2: fn = 2; lst_turn_on(); break;
X default: break;
X }
X }
X}
X
X/* write the active listing to the current listing file, if one is open. */
Xlisting()
X{
X if (lst_fp) {
X int n;
X double flx;
X char flstr[32];
X if (!srch_ison() && lstsrchfld) {
X /* if searching is not on but we are listing the search
X * funtion we must evaluate and log it ourselves here and now.
X * lst_turn_on() insured there is a good function to eval.
X * N.B. if searching IS on, we rely on main() having called
X * srch_eval() BEFORE plot() so it is already evaluated.
X */
X double e;
X char errmsg[128];
X if (execute_expr (&e, errmsg) < 0) {
X f_msg (errmsg);
X lst_turn_off();
X return;
X } else {
X (void) sprintf (flstr, "%g", e);
X (void) flog_log (R_SRCH, C_SRCH, e, flstr);
X }
X }
X
X /* list in order of original selection */
X for (n = 0; n < nlstflds; n++)
X if (flog_get (lstflds[n], &flx, flstr) == 0)
X (void) fprintf (lst_fp, "%s ", flstr);
X (void) fprintf (lst_fp, "\n");
X }
X}
X
Xlisting_prstate (force)
Xint force;
X{
X static last;
X int this = lst_fp != 0;
X
X if (force || this != last) {
X f_string (R_LISTING, C_LISTINGV, this ? " on" : "off");
X last = this;
X }
X}
X
Xlisting_ison()
X{
X return (lst_fp != 0);
X}
X
Xstatic
Xlst_reset()
X{
X int *lp;
X
X for (lp = lstflds; lp < &lstflds[nlstflds]; lp++) {
X (void) flog_delete (*lp);
X *lp = 0;
X }
X nlstflds = 0;
X lstsrchfld = 0;
X}
X
X/* let operator select the fields he wants to have in his listing.
X * register them with flog and keep rcfpack() in lstflds[] array.
X * as a special case, set lstsrchfld if Search field is selected.
X */
Xstatic
Xlst_select_fields()
X{
X static char hlp[] = "move and RETURN to select a field, or q to quit";
X static char sry[] = "Sorry; can not list any more fields.";
X int r = R_UT, c = C_UTV; /* TODO: start where main was? */
X int sf = rcfpack (R_SRCH, C_SRCH, 0);
X char buf[64];
X int rcp;
X int i;
X
X lst_reset();
X for (i = 0; i < MAXLSTFLDS; i++) {
X (void) sprintf(buf,"select field for column %d or q to quit", i+1);
X rcp = sel_fld (r, c, alt_menumask()|F_PLT, buf, hlp);
X if (!rcp)
X break;
X if (flog_add (rcp) < 0) {
X f_msg (sry);
X break;
X }
X lstflds[i] = rcp;
X if (rcp == sf)
X lstsrchfld = 1;
X r = unpackr(rcp);
X c = unpackc(rcp);
X }
X if (i == MAXLSTFLDS)
X f_msg (sry);
X nlstflds = i;
X}
X
Xstatic
Xlst_turn_off ()
X{
X (void) fclose (lst_fp);
X lst_fp = 0;
X listing_prstate(0);
X}
X
X/* turn on listing facility.
X * establish a file to use (and thereby set lst_fp, the "listing-is-on" flag).
X * also check that there is a srch function if it is being used.
X */
Xstatic
Xlst_turn_on ()
X{
X int sf = rcfpack(R_SRCH, C_SRCH, 0);
X char fn[FNLEN], fnq[NC];
X char *optype;
X int n;
X
X /* insure there is a valid srch function if we are to list it */
X for (n = 0; n < nlstflds; n++)
X if (lstflds[n] == sf && !prog_isgood()) {
X f_msg ("Listing search function but it is not defined.");
X return;
X }
X
X /* prompt for file name, giving current as default */
X (void) sprintf (fnq, "file to write <%s>: ", lst_filename);
X f_prompt (fnq);
X n = read_line (fn, sizeof(fn)-1);
X
X /* leave plotting off if type END.
X * reuse same fn if just type \n
X */
X if (n < 0)
X return;
X if (n > 0)
X (void) strcpy (lst_filename, fn);
X
X /* give option to append if file already exists */
X optype = "w";
X if (access (lst_filename, 2) == 0) {
X while (1) {
X f_prompt ("files exists; append or overwrite (a/o)?: ");
X n = read_char();
X if (n == 'a') {
X optype = "a";
X break;
X }
X if (n == 'o')
X break;
X }
X }
X
X /* listing is on if file opens ok */
X lst_fp = fopen (lst_filename, optype);
X if (!lst_fp) {
X (void) sprintf (fnq, "can not open %s: %s", lst_filename, errsys);
X f_msg (fnq);
X } else {
X /* add a title if desired */
X static char tp[] = "Title (q to skip): ";
X f_prompt (tp);
X if (read_line (fnq, PW - sizeof(tp)) > 0)
X (void) fprintf (lst_fp, "%s\n", fnq);
X }
X
X listing_prstate (0);
X}
X
X/* ask operator for a listing file to show. if it's ok, do it.
X */
Xstatic
Xlst_file ()
X{
X char fn[FNLEN], fnq[64];
X FILE *lfp;
X int n;
X
X /* prompt for file name, giving current as default */
X (void) sprintf (fnq, "file to read <%s>: ", lst_filename);
X f_prompt (fnq);
X n = read_line (fn, sizeof(fn)-1);
X
X /* forget it if type END.
X * reuse same fn if just type \n
X */
X if (n < 0)
X return;
X if (n > 0)
X (void) strcpy (lst_filename, fn);
X
X /* show it if file opens ok */
X lfp = fopen (lst_filename, "r");
X if (lfp) {
X display_listing_file (lfp);
X (void) fclose (lfp);
X } else {
X char buf[NC];
X (void) sprintf (buf, "can not open %s: %s", lst_filename, errsys);
X f_prompt (buf);
X (void)read_char();
X }
X}
X
X/* display the given listing file on the screen.
X * allow for files longer than the screen.
X * N.B. do whatever you like but redraw the screen when done.
X */
Xstatic
Xdisplay_listing_file (lfp)
XFILE *lfp;
X{
X char buf[128];
X int n;
X
X c_erase();
X n = 0;
X while (1) {
X (void) fgets (buf, sizeof(buf), lfp);
X if (feof(lfp) || (printf ("%.*s\r", NC, buf), ++n == NR-1)) {
X static char p[] = "[Hit any key to continue or q to quit...]";
X static char ep[] = "[End-of-file. Hit any key to resume...]";
X if (feof(lfp)) {
X f_string (NR, (NC-sizeof(ep))/2, ep);
X (void) read_char();
X break;
X } else {
X f_string (NR, (NC-sizeof(p))/2, p);
X if (read_char() == END)
X break;
X }
X c_erase();
X n = 0;
X }
X }
X
X redraw_screen (2); /* full redraw */
X}
END_OF_FILE
if test 7232 -ne `wc -c <'listing.c'`; then
echo shar: \"'listing.c'\" unpacked with wrong size!
fi
# end of 'listing.c'
fi
if test -f 'mainmenu.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'mainmenu.c'\"
else
echo shar: Extracting \"'mainmenu.c'\" \(6861 characters\)
sed "s/^X//" >'mainmenu.c' <<'END_OF_FILE'
X/* printing routines for the main (upper) screen.
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 PC_GRAPHICS */
X#ifdef PC_GRAPHICS
X#define JOINT 207
X#define VERT 179
X#define HORIZ 205
X#else
X#define JOINT '-'
X#define VERT '|'
X#define HORIZ '-'
X#endif
X
Xmm_borders()
X{
X char line[NC+1], *lp;
X register i;
X
X lp = line;
X for (i = 0; i < NC; i++)
X *lp++ = HORIZ;
X *lp = '\0';
X f_string (R_PLANTAB-1, 1, line);
X for (i = R_TOP; i < R_PLANTAB-1; i++)
X f_char (i, COL2-2, VERT);
X f_char (R_PLANTAB-1, COL2-2, JOINT);
X for (i = R_TOP; i < R_PLANTAB-1; i++)
X f_char (i, COL3-2, VERT);
X f_char (R_PLANTAB-1, COL3-2, JOINT);
X for (i = R_LST; i < R_PLANTAB-1; i++)
X f_char (i, COL4-2, VERT);
X f_char (R_PLANTAB-1, COL4-2, JOINT);
X}
X
X/* print the permanent labels, on the top menu and the planet names
X * in the bottom section.
X */
Xmm_labels()
X{
X f_string (R_TZN, C_TZN, "LT");
X f_string (R_UT, C_UT, "UTC");
X f_string (R_JD, C_JD, "JulianDate");
X f_string (R_LISTING, C_LISTING, "Listing");
X f_string (R_WATCH, C_WATCH, "Watch");
X f_string (R_SRCH, C_SRCH, "Search");
X f_string (R_PLOT, C_PLOT, "Plot");
X f_string (R_ALTM, C_ALTM, "Menu");
X
X f_string (R_LST, C_LST, "LST");
X f_string (R_DAWN, C_DAWN, "Dawn");
X f_string (R_DUSK, C_DUSK, "Dusk");
X f_string (R_LON, C_LON, "NiteLn");
X f_string (R_PAUSE, C_PAUSE, "Pause");
X f_string (R_NSTEP, C_NSTEP, "NStep");
X f_string (R_STPSZ, C_STPSZ, "StpSz");
X
X f_string (R_LAT, C_LAT, "Lat");
X f_string (R_LONG, C_LONG, "Long");
X f_string (R_HEIGHT, C_HEIGHT, "Elev");
X f_string (R_TEMP, C_TEMP, "Temp");
X f_string (R_PRES, C_PRES, "AtmPr");
X f_string (R_TZONE, C_TZONE, "TZ");
X f_string (R_EPOCH, C_EPOCH, "Epoch");
X
X f_string (R_PLANTAB, C_OBJ, "Ob");
X f_string (R_SUN, C_OBJ, "Su");
X f_string (R_MOON, C_OBJ, "Mo");
X f_string (R_MERCURY, C_OBJ, "Me");
X f_string (R_VENUS, C_OBJ, "Ve");
X f_string (R_MARS, C_OBJ, "Ma");
X f_string (R_JUPITER, C_OBJ, "Ju");
X f_string (R_SATURN, C_OBJ, "Sa");
X f_string (R_URANUS, C_OBJ, "Ur");
X f_string (R_NEPTUNE, C_OBJ, "Ne");
X f_string (R_PLUTO, C_OBJ, "Pl");
X f_string (R_OBJX, C_OBJ, "X");
X f_string (R_OBJY, C_OBJ, "Y");
X}
X
X/* print all the time/date/where related stuff: the Now structure.
X * print in a nice order, based on the field locations, as much as possible.
X */
Xmm_now (np, all)
XNow *np;
Xint all;
X{
X char buf[32];
X double lmjd = mjd - tz/24.0;
X double jd = mjd + 2415020L;
X double tmp;
X
X (void) sprintf (buf, "%-3.3s", tznm);
X f_string (R_TZN, C_TZN, buf);
X f_time (R_LT, C_LT, mjd_hr(lmjd));
X f_date (R_LD, C_LD, mjd_day(lmjd));
X
X f_time (R_UT, C_UTV, mjd_hr(mjd));
X f_date (R_UD, C_UD, mjd_day(mjd));
X
X (void) sprintf (buf, "%14.5f", jd);
X (void) flog_log (R_JD, C_JDV, jd, buf);
X f_string (R_JD, C_JDV, buf);
X
X now_lst (np, &tmp);
X f_time (R_LST, C_LSTV, tmp);
X
X if (all) {
X f_gangle (R_LAT, C_LATV, lat);
X f_gangle (R_LONG, C_LONGV, -lng); /* + west */
X
X tmp = height * 2.093e7; /* want to see ft, not earth radii */
X (void) sprintf (buf, "%5g ft", tmp);
X (void) flog_log (R_HEIGHT, C_HEIGHTV, tmp, buf);
X f_string (R_HEIGHT, C_HEIGHTV, buf);
X
X tmp = 9./5.*temp + 32.0; /* want to see degrees F, not C */
X (void) sprintf (buf, "%6g F", tmp);
X (void) flog_log (R_TEMP, C_TEMPV, tmp, buf);
X f_string (R_TEMP, C_TEMPV, buf);
X
X tmp = pressure / 33.86; /* want to see in. Hg, not mBar */
X (void) sprintf (buf, "%5.2f in", tmp);
X (void) flog_log (R_PRES, C_PRESV, tmp, buf);
X f_string (R_PRES, C_PRESV, buf);
X
X f_signtime (R_TZONE, C_TZONEV, tz);
X
X if (epoch == EOD)
X f_string (R_EPOCH, C_EPOCHV, "(OfDate)");
X else {
X mjd_year (epoch, &tmp);
X f_double (R_EPOCH, C_EPOCHV, "%8.1f", tmp);
X }
X }
X
X /* print the calendar for local day, if new month/year. */
X mm_calendar (np, all > 1);
X}
X
X/* display dawn/dusk/length-of-night times.
X */
Xmm_twilight (np, force)
XNow *np;
Xint force;
X{
X double dusk, dawn;
X double tmp;
X int status;
X
X if (!twilight_cir (np, &dawn, &dusk, &status) && !force)
X return;
X
X if (status != 0) {
X f_blanks (R_DAWN, C_DAWNV, 5);
X f_blanks (R_DUSK, C_DUSKV, 5);
X f_string (R_LON, C_LONV, "-----");
X return;
X }
X
X f_mtime (R_DAWN, C_DAWNV, dawn);
X f_mtime (R_DUSK, C_DUSKV, dusk);
X tmp = dawn - dusk; range (&tmp, 24.0);
X f_mtime (R_LON, C_LONV, tmp);
X}
X
Xmm_newcir (y)
Xint y;
X{
X static char ncmsg[] = "NEW CIRCUMSTANCES";
X static char nomsg[] = " ";
X static int last_y = -1;
X
X if (y != last_y) {
X f_string (R_NEWCIR, C_NEWCIR, y ? ncmsg : nomsg);
X last_y = y;
X }
X}
X
Xstatic
Xmm_calendar (np, force)
XNow *np;
Xint force;
X{
X static char *mnames[] = {
X "January", "February", "March", "April", "May", "June",
X "July", "August", "September", "October", "November", "December"
X };
X static int last_m, last_y;
X static double last_tz = -100;
X char str[64];
X int m, y;
X double d;
X int f, nd;
X int r;
X double jd0;
X
X /* get local m/d/y. do nothing if still same month and not forced. */
X mjd_cal (mjd_day(mjd-tz/24.0), &m, &d, &y);
X if (m == last_m && y == last_y && tz == last_tz && !force)
X return;
X last_m = m;
X last_y = y;
X last_tz = tz;
X
X /* find day of week of first day of month */
X cal_mjd (m, 1.0, y, &jd0);
X mjd_dow (jd0, &f);
X if (f < 0) {
X /* can't figure it out - too hard before Gregorian */
X int i;
X for (i = 8; --i >= 0; )
X f_string (R_CAL+i, C_CAL, " ");
X return;
X }
X
X /* print header */
X f_blanks (R_CAL, C_CAL, 20);
X (void) sprintf (str, "%s %4d", mnames[m-1], y);
X f_string (R_CAL, C_CAL + (20 - (strlen(mnames[m-1]) + 5))/2, str);
X f_string (R_CAL+1, C_CAL, "Su Mo Tu We Th Fr Sa");
X
X /* find number of days in this month */
X mjd_dpm (jd0, &nd);
X
X /* print the calendar */
X for (r = 0; r < 6; r++) {
X char row[7*3+1], *rp = row;
X int c;
X for (c = 0; c < 7; c++) {
X int i = r*7+c;
X if (i < f || i >= f + nd)
X (void) sprintf (rp, " ");
X else
X (void) sprintf (rp, "%2d ", i-f+1);
X rp += 3;
X }
X row[sizeof(row)-2] = '\0'; /* don't print last blank; causes wrap*/
X f_string (R_CAL+2+r, C_CAL, row);
X }
X
X /* over print the new and full moons for this month.
X * TODO: don't really know which dates to use here (see moonnf())
X * so try several to be fairly safe. have to go back to 4/29/1988
X * to find the full moon on 5/1 for example.
X */
X mm_nfmoon (jd0-3, tz, m, f);
X mm_nfmoon (jd0+15, tz, m, f);
X}
X
Xstatic
Xmm_nfmoon (jd, tzone, m, f)
Xdouble jd, tzone;
Xint m, f;
X{
X static char nm[] = "NM", fm[] = "FM";
X double dm;
X int mm, ym;
X double jdn, jdf;
X int di;
X
X moonnf (jd, &jdn, &jdf);
X mjd_cal (jdn-tzone/24.0, &mm, &dm, &ym);
X if (m == mm) {
X di = dm + f - 1;
X f_string (R_CAL+2+di/7, C_CAL+3*(di%7), nm);
X }
X mjd_cal (jdf-tzone/24.0, &mm, &dm, &ym);
X if (m == mm) {
X di = dm + f - 1;
X f_string (R_CAL+2+di/7, C_CAL+3*(di%7), fm);
X }
X}
END_OF_FILE
if test 6861 -ne `wc -c <'mainmenu.c'`; then
echo shar: \"'mainmenu.c'\" unpacked with wrong size!
fi
# end of 'mainmenu.c'
fi
if test -f 'moon.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'moon.c'\"
else
echo shar: Extracting \"'moon.c'\" \(5143 characters\)
sed "s/^X//" >'moon.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X/* given the mjd, find the geocentric ecliptic longitude, lam, and latitude,
X * bet, and horizontal parallax, hp for the moon.
X * N.B. series for long and lat are good to about 10 and 3 arcseconds. however,
X * math errors cause up to 100 and 30 arcseconds error, even if use double.
X * why?? suspect highly sensitive nature of difference used to get m1..6.
X * N.B. still need to correct for nutation. then for topocentric location
X * further correct for parallax and refraction.
X */
Xmoon (mjd, lam, bet, hp)
Xdouble mjd;
Xdouble *lam, *bet, *hp;
X{
X double t, t2;
X double ld;
X double ms;
X double md;
X double de;
X double f;
X double n;
X double a, sa, sn, b, sb, c, sc, e, e2, l, g, w1, w2;
X double m1, m2, m3, m4, m5, m6;
X
X t = mjd/36525.;
X t2 = t*t;
X
X m1 = mjd/27.32158213;
X m1 = 360.0*(m1-(long)m1);
X m2 = mjd/365.2596407;
X m2 = 360.0*(m2-(long)m2);
X m3 = mjd/27.55455094;
X m3 = 360.0*(m3-(long)m3);
X m4 = mjd/29.53058868;
X m4 = 360.0*(m4-(long)m4);
X m5 = mjd/27.21222039;
X m5 = 360.0*(m5-(long)m5);
X m6 = mjd/6798.363307;
X m6 = 360.0*(m6-(long)m6);
X
X ld = 270.434164+m1-(.001133-.0000019*t)*t2;
X ms = 358.475833+m2-(.00015+.0000033*t)*t2;
X md = 296.104608+m3+(.009192+.0000144*t)*t2;
X de = 350.737486+m4-(.001436-.0000019*t)*t2;
X f = 11.250889+m5-(.003211+.0000003*t)*t2;
X n = 259.183275-m6+(.002078+.000022*t)*t2;
X
X a = degrad(51.2+20.2*t);
X sa = sin(a);
X sn = sin(degrad(n));
X b = 346.56+(132.87-.0091731*t)*t;
X sb = .003964*sin(degrad(b));
X c = degrad(n+275.05-2.3*t);
X sc = sin(c);
X ld = ld+.000233*sa+sb+.001964*sn;
X ms = ms-.001778*sa;
X md = md+.000817*sa+sb+.002541*sn;
X f = f+sb-.024691*sn-.004328*sc;
X de = de+.002011*sa+sb+.001964*sn;
X e = 1-(.002495+7.52e-06*t)*t;
X e2 = e*e;
X
X ld = degrad(ld);
X ms = degrad(ms);
X n = degrad(n);
X de = degrad(de);
X f = degrad(f);
X md = degrad(md);
X
X l = 6.28875*sin(md)+1.27402*sin(2*de-md)+.658309*sin(2*de)+
X .213616*sin(2*md)-e*.185596*sin(ms)-.114336*sin(2*f)+
X .058793*sin(2*(de-md))+.057212*e*sin(2*de-ms-md)+
X .05332*sin(2*de+md)+.045874*e*sin(2*de-ms)+.041024*e*sin(md-ms);
X l = l-.034718*sin(de)-e*.030465*sin(ms+md)+.015326*sin(2*(de-f))-
X .012528*sin(2*f+md)-.01098*sin(2*f-md)+.010674*sin(4*de-md)+
X .010034*sin(3*md)+.008548*sin(4*de-2*md)-e*.00791*sin(ms-md+2*de)-
X e*.006783*sin(2*de+ms);
X l = l+.005162*sin(md-de)+e*.005*sin(ms+de)+.003862*sin(4*de)+
X e*.004049*sin(md-ms+2*de)+.003996*sin(2*(md+de))+
X .003665*sin(2*de-3*md)+e*.002695*sin(2*md-ms)+
X .002602*sin(md-2*(f+de))+e*.002396*sin(2*(de-md)-ms)-
X .002349*sin(md+de);
X l = l+e2*.002249*sin(2*(de-ms))-e*.002125*sin(2*md+ms)-
X e2*.002079*sin(2*ms)+e2*.002059*sin(2*(de-ms)-md)-
X .001773*sin(md+2*(de-f))-.001595*sin(2*(f+de))+
X e*.00122*sin(4*de-ms-md)-.00111*sin(2*(md+f))+.000892*sin(md-3*de);
X l = l-e*.000811*sin(ms+md+2*de)+e*.000761*sin(4*de-ms-2*md)+
X e2*.000704*sin(md-2*(ms+de))+e*.000693*sin(ms-2*(md-de))+
X e*.000598*sin(2*(de-f)-ms)+.00055*sin(md+4*de)+.000538*sin(4*md)+
X e*.000521*sin(4*de-ms)+.000486*sin(2*md-de);
X l = l+e2*.000717*sin(md-2*ms);
X *lam = ld+degrad(l);
X range (lam, 2*PI);
X
X g = 5.12819*sin(f)+.280606*sin(md+f)+.277693*sin(md-f)+
X .173238*sin(2*de-f)+.055413*sin(2*de+f-md)+.046272*sin(2*de-f-md)+
X .032573*sin(2*de+f)+.017198*sin(2*md+f)+.009267*sin(2*de+md-f)+
X .008823*sin(2*md-f)+e*.008247*sin(2*de-ms-f);
X g = g+.004323*sin(2*(de-md)-f)+.0042*sin(2*de+f+md)+
X e*.003372*sin(f-ms-2*de)+e*.002472*sin(2*de+f-ms-md)+
X e*.002222*sin(2*de+f-ms)+e*.002072*sin(2*de-f-ms-md)+
X e*.001877*sin(f-ms+md)+.001828*sin(4*de-f-md)-e*.001803*sin(f+ms)-
X .00175*sin(3*f);
X g = g+e*.00157*sin(md-ms-f)-.001487*sin(f+de)-e*.001481*sin(f+ms+md)+
X e*.001417*sin(f-ms-md)+e*.00135*sin(f-ms)+.00133*sin(f-de)+
X .001106*sin(f+3*md)+.00102*sin(4*de-f)+.000833*sin(f+4*de-md)+
X .000781*sin(md-3*f)+.00067*sin(f+4*de-2*md);
X g = g+.000606*sin(2*de-3*f)+.000597*sin(2*(de+md)-f)+
X e*.000492*sin(2*de+md-ms-f)+.00045*sin(2*(md-de)-f)+
X .000439*sin(3*md-f)+.000423*sin(f+2*(de+md))+
X .000422*sin(2*de-f-3*md)-e*.000367*sin(ms+f+2*de-md)-
X e*.000353*sin(ms+f+2*de)+.000331*sin(f+4*de);
X g = g+e*.000317*sin(2*de+f-ms+md)+e2*.000306*sin(2*(de-ms)-f)-
X .000283*sin(md+3*f);
X w1 = .0004664*cos(n);
X w2 = .0000754*cos(c);
X *bet = degrad(g)*(1-w1-w2);
X
X *hp = .950724+.051818*cos(md)+.009531*cos(2*de-md)+.007843*cos(2*de)+
X .002824*cos(2*md)+.000857*cos(2*de+md)+e*.000533*cos(2*de-ms)+
X e*.000401*cos(2*de-md-ms)+e*.00032*cos(md-ms)-.000271*cos(de)-
X e*.000264*cos(ms+md)-.000198*cos(2*f-md);
X *hp = *hp+.000173*cos(3*md)+.000167*cos(4*de-md)-e*.000111*cos(ms)+
X .000103*cos(4*de-2*md)-.000084*cos(2*md-2*de)-
X e*.000083*cos(2*de+ms)+.000079*cos(2*de+2*md)+.000072*cos(4*de)+
X e*.000064*cos(2*de-ms+md)-e*.000063*cos(2*de+ms-md)+
X e*.000041*cos(ms+de);
X *hp = *hp+e*.000035*cos(2*md-ms)-.000033*cos(3*md-2*de)-
X .00003*cos(md+de)-.000029*cos(2*(f-de))-e*.000029*cos(2*md+ms)+
X e2*.000026*cos(2*(de-ms))-.000023*cos(2*(f-de)+md)+
X e*.000019*cos(4*de-ms-md);
X *hp = degrad(*hp);
X}
END_OF_FILE
if test 5143 -ne `wc -c <'moon.c'`; then
echo shar: \"'moon.c'\" unpacked with wrong size!
fi
# end of 'moon.c'
fi
if test -f 'riset_c.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'riset_c.c'\"
else
echo shar: Extracting \"'riset_c.c'\" \(8287 characters\)
sed "s/^X//" >'riset_c.c' <<'END_OF_FILE'
X/* find rise and set circumstances, ie, riset_cir() and related functions. */
X
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X#include "circum.h"
X#include "screen.h" /* just for SUN and MOON */
X
X#define TRACE(x) {FILE *fp = fopen("trace","a"); fprintf x; fclose(fp);}
X
X#define STDREF degrad(34./60.) /* nominal horizon refraction amount */
X#define TWIREF degrad(18.) /* twilight horizon displacement */
X#define TMACC (15./3600.) /* convergence accuracy, hours */
X
X/* find where and when a body, p, will rise and set and
X * its transit circumstances. all times are local, angles rads e of n.
X * return 0 if just returned same stuff as previous call, else 1 if new.
X * status is set from the RS_* #defines in circum.h.
X * also used to find astro twilight by calling with dis of 18 degrees.
X */
Xriset_cir (p, np, force, hzn, ltr, lts, ltt, azr, azs, altt, status)
Xint p; /* one of the body defines in astro.h or screen.h */
XNow *np;
Xint force; /* set !=0 to force computations */
Xint hzn; /* STDHZN or ADPHZN or TWILIGHT */
Xdouble *ltr, *lts; /* local rise and set times */
Xdouble *ltt; /* local transit time */
Xdouble *azr, *azs; /* local rise and set azimuths, rads e of n */
Xdouble *altt; /* local altitude at transit */
Xint *status; /* one or more of the RS_* defines */
X{
X typedef struct {
X Now l_now;
X double l_ltr, l_lts, l_ltt, l_azr, l_azs, l_altt;
X int l_hzn;
X int l_status;
X } Last;
X /* must be in same order as the astro.h/screen.h #define's */
X static Last last[NOBJ] = {
X {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD},
X {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD}, {NOMJD}
X };
X Last *lp;
X int new;
X
X lp = last + p;
X if (!force && same_cir (np, &lp->l_now) && same_lday (np, &lp->l_now)
X && lp->l_hzn == hzn) {
X *ltr = lp->l_ltr;
X *lts = lp->l_lts;
X *ltt = lp->l_ltt;
X *azr = lp->l_azr;
X *azs = lp->l_azs;
X *altt = lp->l_altt;
X *status = lp->l_status;
X new = 0;
X } else {
X *status = 0;
X iterative_riset (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status);
X lp->l_ltr = *ltr;
X lp->l_lts = *lts;
X lp->l_ltt = *ltt;
X lp->l_azr = *azr;
X lp->l_azs = *azs;
X lp->l_altt = *altt;
X lp->l_status = *status;
X lp->l_hzn = hzn;
X lp->l_now = *np;
X new = 1;
X }
X return (new);
X}
X
Xstatic
Xiterative_riset (p, np, hzn, ltr, lts, ltt, azr, azs, altt, status)
Xint p;
XNow *np;
Xint hzn;
Xdouble *ltr, *lts, *ltt; /* local times of rise, set and transit */
Xdouble *azr, *azs, *altt;/* local azimuths of rise, set and transit altitude */
Xint *status;
X{
X#define MAXPASSES 6
X double lstr, lsts, lstt; /* local sidereal times of rising/setting */
X double mjd0; /* mjd estimates of rise/set event */
X double lnoon; /* mjd of local noon */
X double x; /* discarded tmp value */
X Now n; /* just used to call now_lst() */
X double lst; /* lst at local noon */
X double diff, lastdiff; /* iterative improvement to mjd0 */
X int pass;
X int rss;
X
X /* first approximation is to find rise/set times of a fixed object
X * in its position at local noon.
X */
X lnoon = mjd_day(mjd - tz/24.0) + (12.0 + tz)/24.0; /*mjd of local noon*/
X n.n_mjd = lnoon;
X n.n_lng = lng;
X now_lst (&n, &lst); /* lst at local noon */
X mjd0 = lnoon;
X stationary_riset (p,mjd0,np,hzn,&lstr,&lsts,&lstt,&x,&x,&x,&rss);
X chkrss:
X switch (rss) {
X case 0: break;
X case 1: *status = RS_NEVERUP; return;
X case -1: *status = RS_CIRCUMPOLAR; goto transit;
X default: *status = RS_ERROR; return;
X }
X
X /* find a better approximation to the rising circumstances based on
X * more passes, each using a "fixed" object at the location at
X * previous approximation of the rise time.
X */
X lastdiff = 1000.0;
X for (pass = 1; pass < MAXPASSES; pass++) {
X diff = (lstr - lst)*SIDRATE; /* next guess at rise time wrt noon */
X if (diff > 12.0)
X diff -= 24.0*SIDRATE; /* not tomorrow, today */
X else if (diff < -12.0)
X diff += 24.0*SIDRATE; /* not yesterday, today */
X mjd0 = lnoon + diff/24.0; /* next guess at mjd of rise */
X stationary_riset (p,mjd0,np,hzn,&lstr,&x,&x,azr,&x,&x,&rss);
X if (rss != 0) goto chkrss;
X if (fabs (diff - lastdiff) < TMACC)
X break;
X lastdiff = diff;
X }
X if (pass == MAXPASSES)
X *status |= RS_NORISE; /* didn't converge - no rise today */
X else {
X *ltr = 12.0 + diff;
X if (p != MOON &&
X (*ltr <= 24.0*(1.0-SIDRATE) || *ltr >= 24.0*SIDRATE))
X *status |= RS_2RISES;
X }
X
X /* find a better approximation to the setting circumstances based on
X * more passes, each using a "fixed" object at the location at
X * previous approximation of the set time.
X */
X lastdiff = 1000.0;
X for (pass = 1; pass < MAXPASSES; pass++) {
X diff = (lsts - lst)*SIDRATE; /* next guess at set time wrt noon */
X if (diff > 12.0)
X diff -= 24.0*SIDRATE; /* not tomorrow, today */
X else if (diff < -12.0)
X diff += 24.0*SIDRATE; /* not yesterday, today */
X mjd0 = lnoon + diff/24.0; /* next guess at mjd of set */
X stationary_riset (p,mjd0,np,hzn,&x,&lsts,&x,&x,azs,&x,&rss);
X if (rss != 0) goto chkrss;
X if (fabs (diff - lastdiff) < TMACC)
X break;
X lastdiff = diff;
X }
X if (pass == MAXPASSES)
X *status |= RS_NOSET; /* didn't converge - no set today */
X else {
X *lts = 12.0 + diff;
X if (p != MOON &&
X (*lts <= 24.0*(1.0-SIDRATE) || *lts >= 24.0*SIDRATE))
X *status |= RS_2SETS;
X }
X
X transit:
X /* find a better approximation to the transit circumstances based on
X * more passes, each using a "fixed" object at the location at
X * previous approximation of the transit time.
X */
X lastdiff = 1000.0;
X for (pass = 1; pass < MAXPASSES; pass++) {
X diff = (lstt - lst)*SIDRATE; /*next guess at transit time wrt noon*/
X if (diff > 12.0)
X diff -= 24.0*SIDRATE; /* not tomorrow, today */
X else if (diff < -12.0)
X diff += 24.0*SIDRATE; /* not yesterday, today */
X mjd0 = lnoon + diff/24.0; /* next guess at mjd of transit */
X stationary_riset (p,mjd0,np,hzn,&x,&x,&lstt,&x,&x,altt,&rss);
X if (fabs (diff - lastdiff) < TMACC)
X break;
X lastdiff = diff;
X }
X if (pass == MAXPASSES)
X *status |= RS_NOTRANS; /* didn't converge - no transit today */
X else {
X *ltt = 12.0 + diff;
X if (p != MOON &&
X (*ltt <= 24.0*(1.0-SIDRATE) || *ltt >= 24.0*SIDRATE))
X *status |= RS_2TRANS;
X }
X}
X
Xstatic
Xstationary_riset (p, mjd0, np, hzn, lstr, lsts, lstt, azr, azs, altt, status)
Xint p;
Xdouble mjd0;
XNow *np;
Xint hzn;
Xdouble *lstr, *lsts, *lstt;
Xdouble *azr, *azs, *altt;
Xint *status;
X{
X extern void bye();
X double dis;
X Now n;
X Sky s;
X
X /* find object p's topocentric ra/dec at mjd0
X * (this must include parallax)
X */
X n = *np;
X n.n_mjd = mjd0;
X (void) body_cir (p, 0.0, &n, &s);
X if (epoch != EOD)
X precess (epoch, mjd0, &s.s_ra, &s.s_dec);
X if (s.s_edist > 0) {
X /* parallax, if we can */
X double ehp, lst, ha;
X if (p == MOON)
X ehp = asin (6378.14/s.s_edist);
X else
X ehp = (2.*6378./146e6)/s.s_edist;
X now_lst (&n, &lst);
X ha = hrrad(lst) - s.s_ra;
X ta_par (ha, s.s_dec, lat, height, ehp, &ha, &s.s_dec);
X s.s_ra = hrrad(lst) - ha;
X range (&s.s_ra, 2*PI);
X }
X
X switch (hzn) {
X case STDHZN:
X /* nominal atmospheric refraction.
X * then add nominal moon or sun semi-diameter, as appropriate.
X * other objects assumes to be negligibly small.
X */
X dis = STDREF;
X if (p == MOON || p == SUN)
X dis += degrad (32./60./2.);
X break;
X case TWILIGHT:
X if (p != SUN) {
X f_msg ("Non-sun twilight bug!");
X bye();
X }
X dis = TWIREF;
X break;
X case ADPHZN:
X /* adaptive includes actual refraction conditions and also
X * includes object's semi-diameter.
X */
X unrefract (pressure, temp, 0.0, &dis);
X dis = -dis;
X dis += degrad(s.s_size/3600./2.0);
X break;
X }
X
X riset (s.s_ra, s.s_dec, lat, dis, lstr, lsts, azr, azs, status);
X transit (s.s_ra, s.s_dec, np, lstt, altt);
X}
X
X
X/* find when and how hi object at (r,d) is when it transits. */
Xstatic
Xtransit (r, d, np, lstt, altt)
Xdouble r, d; /* ra and dec, rads */
XNow *np; /* for refraction info */
Xdouble *lstt; /* local sidereal time of transit */
Xdouble *altt; /* local, refracted, altitude at time of transit */
X{
X *lstt = radhr(r);
X *altt = PI/2 - lat + d;
X if (*altt > PI/2)
X *altt = PI - *altt;
X refract (pressure, temp, *altt, altt);
X}
END_OF_FILE
if test 8287 -ne `wc -c <'riset_c.c'`; then
echo shar: \"'riset_c.c'\" unpacked with wrong size!
fi
# end of 'riset_c.c'
fi
if test -f 'screen.h' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'screen.h'\"
else
echo shar: Extracting \"'screen.h'\" \(5258 characters\)
sed "s/^X//" >'screen.h' <<'END_OF_FILE'
X/* screen layout details
X *
X * it looks better if the fields are drawn in some nice order so it you
X * rearrange the fields, check the menu printing functions.
X */
X
X/* size of screen */
X#define NR 24
X#define NC 80
X
X#define ASPECT (4./3.) /* screen width to height dimensions ratio */
X
X#define GAP 6 /* gap between field name and value */
X
X#define COL1 1
X#define COL2 27
X#define COL3 44
X#define COL4 61 /* calendar */
X
X#define R_PROMPT 1 /* prompt row */
X#define C_PROMPT COL1
X
X#define R_NEWCIR 2
X#define C_NEWCIR ((NC-17)/2) /* 17 is length of the message */
X
X#define R_TOP 3 /* first row of top menu items */
X
X#define R_TZN (R_TOP+0)
X#define C_TZN COL1
X#define R_LT R_TZN
X#define C_LT (C_TZN+GAP-2)
X#define R_LD R_TZN
X#define C_LD (C_TZN+13)
X
X#define R_UT (R_TOP+1)
X#define C_UT COL1
X#define C_UTV (C_UT+GAP-2)
X#define R_UD R_UT
X#define C_UD (C_UT+13)
X
X#define R_JD (R_TOP+2)
X#define C_JD COL1
X#define C_JDV (C_JD+GAP+3)
X
X#define R_LST (R_TOP)
X#define C_LST COL2
X#define C_LSTV (C_LST+GAP)
X
X#define R_LAT (R_TOP+0)
X#define C_LAT COL3
X#define C_LATV (C_LAT+4)
X
X#define R_DAWN (R_TOP+2)
X#define C_DAWN COL2
X#define C_DAWNV (C_DAWN+GAP+3)
X
X#define R_STPSZ (R_TOP+7)
X#define C_STPSZ COL2
X#define C_STPSZV (C_STPSZ+GAP-1)
X
X#define R_HEIGHT (R_TOP+2)
X#define C_HEIGHT COL3
X#define C_HEIGHTV (C_HEIGHT+GAP)
X
X#define R_PRES (R_TOP+4)
X#define C_PRES COL3
X#define C_PRESV (C_PRES+GAP)
X
X#define R_WATCH (R_TOP+3)
X#define C_WATCH COL1
X
X#define R_LISTING (R_TOP+4)
X#define C_LISTING COL1
X#define C_LISTINGV (C_LISTING+20)
X
X#define R_SRCH (R_TOP+5)
X#define C_SRCH COL1
X#define C_SRCHV (C_SRCH+16)
X
X#define R_PLOT (R_TOP+6)
X#define C_PLOT COL1
X#define C_PLOTV (C_PLOT+20)
X
X#define R_ALTM (R_TOP+7)
X#define C_ALTM COL1
X#define C_ALTMV (C_ALTM+10)
X
X#define R_TZONE (R_TOP+5)
X#define C_TZONE COL3
X#define C_TZONEV (C_TZONE+GAP-1)
X
X#define R_LONG (R_TOP+1)
X#define C_LONG COL3
X#define C_LONGV (C_LONG+4)
X
X#define R_DUSK (R_TOP+3)
X#define C_DUSK COL2
X#define C_DUSKV (C_DUSK+GAP+3)
X
X#define R_NSTEP (R_TOP+6)
X#define C_NSTEP COL2
X#define C_NSTEPV (C_NSTEP+GAP)
X
X#define R_TEMP (R_TOP+3)
X#define C_TEMP COL3
X#define C_TEMPV (C_TEMP+GAP)
X
X#define R_EPOCH (R_TOP+6)
X#define C_EPOCH COL3
X#define C_EPOCHV (C_EPOCH+GAP)
X
X#define R_PAUSE (R_TOP+7)
X#define C_PAUSE COL3
X#define C_PAUSEV (C_PAUSE+GAP)
X
X#define R_MNUDEP (R_TOP+6)
X#define C_MNUDEP COL3
X#define C_MNUDEPV (C_EPOCH+GAP)
X
X#define R_LON (R_TOP+4)
X#define C_LON COL2
X#define C_LONV (C_LON+GAP+3)
X
X#define R_CAL R_TOP
X#define C_CAL COL4
X
X/* planet rows, for all menus */
X#define R_PLANTAB (R_TOP+9)
X#define R_SUN (R_PLANTAB+1)
X#define R_MOON (R_PLANTAB+2)
X#define R_MERCURY (R_PLANTAB+3)
X#define R_VENUS (R_PLANTAB+4)
X#define R_MARS (R_PLANTAB+5)
X#define R_JUPITER (R_PLANTAB+6)
X#define R_SATURN (R_PLANTAB+7)
X#define R_URANUS (R_PLANTAB+8)
X#define R_NEPTUNE (R_PLANTAB+9)
X#define R_PLUTO (R_PLANTAB+10)
X#define R_OBJX (R_PLANTAB+11)
X#define R_OBJY (R_PLANTAB+12)
X
X/* menu 1 info table */
X#define C_OBJ 1
X#define C_RA 4
X#define C_DEC 12
X#define C_AZ 19
X#define C_ALT 26
X#define C_HLONG 33
X#define C_HLAT 40
X#define C_EDIST 47
X#define C_SDIST 54
X#define C_ELONG 61
X#define C_SIZE 68
X#define C_MAG 73
X#define C_PHASE 78
X
X/* menu 2 screen items */
X#define C_RISETM 7
X#define C_RISEAZ 18
X#define C_TRANSTM 29
X#define C_TRANSALT 40
X#define C_SETTM 51
X#define C_SETAZ 62
X#define C_TUP 73
X
X/* menu 3 items */
X#define C_SUN 4
X#define C_MOON 10
X#define C_MERCURY 17
X#define C_VENUS 23
X#define C_MARS 30
X#define C_JUPITER 36
X#define C_SATURN 43
X#define C_URANUS 49
X#define C_NEPTUNE 56
X#define C_PLUTO 62
X#define C_OBJX 69
X#define C_OBJY 75
X
X#define PW (NC-C_PROMPT+1) /* total prompt line width */
X
X/* macros to pack a row/col and menu selection flags all into an int.
X * (use this rather than a structure because we can compare them so easily.
X * could use bit fields and a union, but then can't init them or use switch.)
X * bit field defs: [15..14]=menu [13..12]=flags [11..7]=row [6..0]=column.
X * see sel_fld.c.
X * F_MNUX also used in main to manage which bottom menu is up.
X */
X#define F_CHG (1<<12) /* field may be picked for changing */
X#define F_PLT (1<<13) /* field may be picked for plotting or listng */
X#define F_MMNU (0<<14) /* field is on main menu */
X#define F_MNU1 (1<<14) /* field is on menu 1 */
X#define F_MNU2 (2<<14) /* field is on menu 2 */
X#define F_MNU3 (3<<14) /* field is on menu 3 */
X#define rcfpack(r,c,f) ((f) | ((r) << 7) | (c))
X#define unpackr(p) (((p) >> 7) & 0x1f)
X#define unpackc(p) ((p) & 0x7f)
X#define unpackrc(p) ((p) & 0xfff)
X#define tstpackf(p,f) (((p) & ((f)&0x3000)) && \
X (((p)&0xc000) == ((f)&0xc000) || ((p)&0xc000) == 0))
X
X/* additions to the planet defines from astro.h.
X * must not conflict, and must fit in range 0..15.
X */
X#define SUN (PLUTO+1)
X#define MOON (PLUTO+2)
X#define OBJX (PLUTO+3) /* the user-defined object */
X#define OBJY (PLUTO+4) /* the user-defined object */
X#define NOBJ (OBJY+1) /* total number of objects */
X
X#define cntrl(x) ((x) & 037)
X#define QUIT cntrl('d') /* char to exit program */
X#define HELP '?' /* char to give help message */
X#define REDRAW cntrl('l') /* char to redraw (like vi) */
X#define VERSION cntrl('v') /* char to display version number */
X#define END 'q' /* char to quit current mode */
END_OF_FILE
if test 5258 -ne `wc -c <'screen.h'`; then
echo shar: \"'screen.h'\" unpacked with wrong size!
fi
# end of 'screen.h'
fi
if test -f 'srch.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'srch.c'\"
else
echo shar: Extracting \"'srch.c'\" \(8127 characters\)
sed "s/^X//" >'srch.c' <<'END_OF_FILE'
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 static int fn; /* start with 0, then remember for next time */
X
X /* let op select algorithm, edit, set accuracy
X * or stop if currently searching
X * algorithms require a function.
X */
X ask:
X switch (popup(chcs, fn, srch_f ? 6 : 5)) {
X case 0: fn = 0;
X if (expbuf[0] == '\0')
X set_function();
X srch_f = expbuf[0] ? srch_minmax : 0;
X if (srch_f)
X break;
X else
X goto ask;
X case 1: fn = 1;
X if (expbuf[0] == '\0')
X set_function();
X srch_f = expbuf[0] ? srch_solve0 : 0;
X if (srch_f)
X break;
X else
X goto ask;
X case 2: fn = 2;
X if (expbuf[0] == '\0')
X set_function();
X srch_f = expbuf[0] ? srch_binary : 0;
X if (srch_f)
X break;
X else
X goto ask;
X case 3: fn = 3; srch_f = 0; set_function(); goto ask;
X case 4: fn = 4; srch_f = 0; set_accuracy(); 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 (void) 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 (void) 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 x_1, x_2, x_3; /* keep in increasing order */
X static double y_1, y_2, y_3;
X double xm, a, b;
X
X if (srch_tmscalled == 0) {
X base = mjd;
X x_1 = 0.0;
X y_1 = v;
X return (0);
X }
X mjd -= base;
X if (srch_tmscalled == 1) {
X /* put in one of first two slots */
X if (mjd < x_1) {
X x_2 = x_1; y_2 = y_1;
X x_1 = mjd; y_1 = v;
X } else {
X x_2 = mjd; y_2 = v;
X }
X return (0);
X }
X if (srch_tmscalled == 2 || fabs(mjd - x_1) < fabs(mjd - x_3)) {
X /* closer to x_1 so discard x_3.
X * or if it's our third value we know to "discard" x_3.
X */
X if (mjd > x_2) {
X x_3 = mjd; y_3 = v;
X } else {
X x_3 = x_2; y_3 = y_2;
X if (mjd > x_1) {
X x_2 = mjd; y_2 = v;
X } else {
X x_2 = x_1; y_2 = y_1;
X x_1 = mjd; y_1 = v;
X }
X }
X if (srch_tmscalled == 2)
X return (0);
X } else {
X /* closer to x_3 so discard x_1 */
X if (mjd < x_2) {
X x_1 = mjd; y_1 = v;
X } else {
X x_1 = x_2; y_1 = y_2;
X if (mjd < x_3) {
X x_2 = mjd; y_2 = v;
X } else {
X x_2 = x_3; y_2 = y_3;
X x_3 = mjd; y_3 = v;
X }
X }
X }
X
X#ifdef TRACEMM
X { char buf[NC];
X sprintf (buf, "x_1=%g y_1=%g x_2=%g y_2=%g x_3=%g y_3=%g",
X x_1, y_1, x_2, y_2, x_3, y_3);
X f_msg (buf);
X }
X#endif
X a = y_1*(x_2-x_3) - y_2*(x_1-x_3) + y_3*(x_1-x_2);
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 = (x_1*x_1)*(y_2-y_3) - (x_2*x_2)*(y_1-y_3) + (x_3*x_3)*(y_1-y_2);
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, x_1; /* x(n-1) and x(n) */
X static double y_0, y_1; /* y(n-1) and y(n) */
X double x_2; /* x(n+1) */
X double df; /* y(n) - y(n-1) */
X
X switch (srch_tmscalled) {
X case 0: x0 = mjd; y_0 = v; return(0);
X case 1: x_1 = mjd; y_1 = v; break;
X default: x0 = x_1; y_0 = y_1; x_1 = mjd; y_1 = v; break;
X }
X
X df = y_1 - y_0;
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 x_2 = x_1 - y_1*(x_1-x0)/df;
X *tmincp = (x_2 - 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}
END_OF_FILE
if test 8127 -ne `wc -c <'srch.c'`; then
echo shar: \"'srch.c'\" unpacked with wrong size!
fi
# end of 'srch.c'
fi
echo shar: End of archive 2 \(of 6\).
cp /dev/null ark2isdone
MISSING=""
for I in 1 2 3 4 5 6 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 6 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0