allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc) (11/28/89)
Posting-number: Volume 9, Issue 39
Submitted-by: ecd@umn-cs.cs.umn.edu@ncs-med.UUCP (Elwood C. Downey)
Archive-name: ephem2/part02
# This is a "shell archive" file; run it with sh.
# This is file 2.
echo x screen.h
cat > screen.h << 'xXx'
/* screen layout details
*
* it looks better if the fields are drawn in some nice order so it you
* rearrange the fields, check the menu printing functions.
*/
/* size of screen */
#define NR 24
#define NC 80
#define ASPECT (4./3.) /* screen width to height dimensions ratio */
#define GAP 6 /* gap between field name and value */
#define COL1 1
#define COL2 27
#define COL3 44
#define COL4 61 /* calendar */
#define R_PROMPT 1 /* prompt row */
#define C_PROMPT COL1
#define R_NEWCIR 2
#define C_NEWCIR ((NC-17)/2) /* 17 is length of the message */
#define R_TOP 3 /* first row of top menu items */
#define R_TZN (R_TOP+0)
#define C_TZN COL1
#define R_LT R_TZN
#define C_LT (C_TZN+GAP-2)
#define R_LD R_TZN
#define C_LD (C_TZN+13)
#define R_UT (R_TOP+1)
#define C_UT COL1
#define C_UTV (C_UT+GAP-2)
#define R_UD R_UT
#define C_UD (C_UT+13)
#define R_JD (R_TOP+2)
#define C_JD COL1
#define C_JDV (C_JD+GAP+3)
#define R_LST (R_TOP)
#define C_LST COL2
#define C_LSTV (C_LST+GAP)
#define R_LAT (R_TOP+0)
#define C_LAT COL3
#define C_LATV (C_LAT+4)
#define R_DAWN (R_TOP+2)
#define C_DAWN COL2
#define C_DAWNV (C_DAWN+GAP+3)
#define R_STPSZ (R_TOP+7)
#define C_STPSZ COL2
#define C_STPSZV (C_STPSZ+GAP-1)
#define R_HEIGHT (R_TOP+2)
#define C_HEIGHT COL3
#define C_HEIGHTV (C_HEIGHT+GAP)
#define R_PRES (R_TOP+4)
#define C_PRES COL3
#define C_PRESV (C_PRES+GAP)
#define R_WATCH (R_TOP+4)
#define C_WATCH COL1
#define R_SRCH (R_TOP+5)
#define C_SRCH COL1
#define C_SRCHV (C_SRCH+16)
#define R_PLOT (R_TOP+6)
#define C_PLOT (COL1)
#define C_PLOTV (C_PLOT+20)
#define R_ALTM (R_TOP+7)
#define C_ALTM COL1
#define C_ALTMV (C_ALTM+10)
#define R_TZONE (R_TOP+5)
#define C_TZONE COL3
#define C_TZONEV (C_TZONE+GAP-1)
#define R_LONG (R_TOP+1)
#define C_LONG COL3
#define C_LONGV (C_LONG+4)
#define R_DUSK (R_TOP+3)
#define C_DUSK COL2
#define C_DUSKV (C_DUSK+GAP+3)
#define R_NSTEP (R_TOP+6)
#define C_NSTEP COL2
#define C_NSTEPV (C_NSTEP+GAP)
#define R_TEMP (R_TOP+3)
#define C_TEMP COL3
#define C_TEMPV (C_TEMP+GAP)
#define R_EPOCH (R_TOP+6)
#define C_EPOCH COL3
#define C_EPOCHV (C_EPOCH+GAP)
#define R_MNUDEP (R_TOP+6)
#define C_MNUDEP COL3
#define C_MNUDEPV (C_EPOCH+GAP)
#define R_LON (R_TOP+4)
#define C_LON COL2
#define C_LONV (C_LON+GAP+3)
#define R_CAL R_TOP
#define C_CAL COL4
/* menu 1 info table */
#define R_PLANTAB (R_TOP+9)
#define R_SUN (R_PLANTAB+2)
#define R_MOON (R_PLANTAB+3)
#define R_MERCURY (R_PLANTAB+4)
#define R_VENUS (R_PLANTAB+5)
#define R_MARS (R_PLANTAB+6)
#define R_JUPITER (R_PLANTAB+7)
#define R_SATURN (R_PLANTAB+8)
#define R_URANUS (R_PLANTAB+9)
#define R_NEPTUNE (R_PLANTAB+10)
#define R_PLUTO (R_PLANTAB+11)
#define R_OBJX (R_PLANTAB+12)
#define C_OBJ 1
#define C_RA 4
#define C_DEC 12
#define C_AZ 19
#define C_ALT 26
#define C_HLONG 33
#define C_HLAT 40
#define C_EDIST 47
#define C_SDIST 54
#define C_ELONG 61
#define C_SIZE 68
#define C_MAG 73
#define C_PHASE 78
/* menu 2 screen items */
#define C_RISETM 10
#define C_RISEAZ 18
#define C_TRANSTM 31
#define C_TRANSALT 39
#define C_SETTM 52
#define C_SETAZ 60
#define C_TUP 72
/* menu 3 items */
#define C_SUN 4
#define C_MOON 11
#define C_MERCURY 18
#define C_VENUS 25
#define C_MARS 32
#define C_JUPITER 39
#define C_SATURN 46
#define C_URANUS 53
#define C_NEPTUNE 60
#define C_PLUTO 67
#define C_OBJX 74
#define MAXOBJXNM 2 /* object x's name. does not include trail 0 */
#define PW (NC-C_PROMPT+1) /* total prompt line width */
/* macros to pack a row/col and menu selection flags all into an int.
* (use this rather than a structure because we can compare them so easily.
* could use bit fields and a union, but then can't init them or use switch.)
* bit field defs: [15..14]=menu [13..12]=flags [11..7]=row [6..0]=column.
* see sel_fld.c.
* F_MNUX also used in main to manage which bottom menu is up.
*/
#define F_CHG (1<<12) /* field may be picked for changing */
#define F_PLT (1<<13) /* field may be picked for plotting */
#define F_MMNU (0<<14) /* field is on main menu */
#define F_MNU1 (1<<14) /* field is on menu 1 */
#define F_MNU2 (2<<14) /* field is on menu 2 */
#define F_MNU3 (3<<14) /* field is on menu 3 */
#define rcfpack(r,c,f) ((f) | ((r) << 7) | (c))
#define unpackr(p) (((p) >> 7) & 0x1f)
#define unpackc(p) ((p) & 0x7f)
#define unpackrc(p) ((p) & 0xfff)
#define tstpackf(p,f) (((p) & ((f)&0x3000)) && \
(((p)&0xc000) == ((f)&0xc000) || ((p)&0xc000) == 0))
/* additions to the planet defines from astro.h.
* must not conflict, and must fit in range 0..15.
*/
#define SUN (PLUTO+1)
#define MOON (PLUTO+2)
#define OBJX (PLUTO+3) /* the user-defined object */
#define NOBJ (OBJX+1) /* total number of objects */
#define cntrl(x) ((x) & 037)
#define QUIT cntrl('d') /* char to exit program */
#define HELP '?' /* char to give help message */
#define REDRAW cntrl('l') /* char to redraw (like vi) */
#define VERSION cntrl('v') /* char to display version number */
#define END 'q' /* char to quit current mode */
xXx
echo x aa_hadec.c
cat > aa_hadec.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given latitude (n+, radians), lat, altitude (up+, radians), alt, and
* azimuth (angle round to the east from north+, radians),
* return hour angle (radians), ha, and declination (radians), dec.
*/
aa_hadec (lat, alt, az, ha, dec)
double lat;
double alt, az;
double *ha, *dec;
{
aaha_aux (lat, az, alt, ha, dec);
}
/* given latitude (n+, radians), lat, hour angle (radians), ha, and declination
* (radians), dec,
* return altitude (up+, radians), alt, and
* azimuth (angle round to the east from north+, radians),
*/
hadec_aa (lat, ha, dec, alt, az)
double lat;
double ha, dec;
double *alt, *az;
{
aaha_aux (lat, ha, dec, az, alt);
}
/* the actual formula is the same for both transformation directions so
* do it here once for each way.
* N.B. all arguments are in radians.
*/
static
aaha_aux (lat, x, y, p, q)
double lat;
double x, y;
double *p, *q;
{
static double lastlat = -1000.;
static double sinlastlat, coslastlat;
double sy, cy;
double sx, cx;
double sq, cq;
double a;
double cp;
/* latitude doesn't change much, so try to reuse the sin and cos evals.
*/
if (lat != lastlat) {
sinlastlat = sin (lat);
coslastlat = cos (lat);
lastlat = lat;
}
sy = sin (y);
cy = cos (y);
sx = sin (x);
cx = cos (x);
/* define GOODATAN2 if atan2 returns full range -PI through +PI.
*/
#ifdef GOODATAN2
*q = asin ((sy*sinlastlat) + (cy*coslastlat*cx));
*p = atan2 (-cy*sx, -cy*cx*sinlastlat + sy*coslastlat);
#else
#define EPS (1e-20)
sq = (sy*sinlastlat) + (cy*coslastlat*cx);
*q = asin (sq);
cq = cos (*q);
a = coslastlat*cq;
if (a > -EPS && a < EPS)
a = a < 0 ? -EPS : EPS; /* avoid / 0 */
cp = (sy - (sinlastlat*sq))/a;
if (cp >= 1.0) /* the /a can be slightly > 1 */
*p = 0.0;
else if (cp <= -1.0)
*p = PI;
else
*p = acos ((sy - (sinlastlat*sq))/a);
if (sx>0) *p = 2.0*PI - *p;
#endif
}
xXx
echo x altmenus.c
cat > altmenus.c << 'xXx'
/* printing routines for the three alternative bottom half menus,
* "menu1", "menu2" and "menu3".
*/
#include <stdio.h>
#include <math.h>
#include "astro.h"
#include "circum.h"
#include "screen.h"
static int altmenu = F_MNU1; /* which alternate menu is up; one of F_MNUi */
static int alt2_stdhzn; /* whether to use STDHZN (aot ADPHZN) horizon algthm */
static int alt3_geoc; /* whether to use geocentric (aot topocentric) vantage*/
/* table of screen rows given a body #define from astro/h or screen.h */
static short bodyrow[NOBJ] = {
R_MERCURY, R_VENUS, R_MARS, R_JUPITER, R_SATURN,
R_URANUS, R_NEPTUNE, R_PLUTO, R_SUN, R_MOON, R_OBJX
};
/* table of screen cols for third menu format, given body #define ... */
static short bodycol[NOBJ] = {
C_MERCURY, C_VENUS, C_MARS, C_JUPITER, C_SATURN,
C_URANUS, C_NEPTUNE, C_PLUTO, C_SUN, C_MOON, C_OBJX
};
/* let op decide which alternate menu should be up,
* including any menu-specific setup they might require.
* return 0 if things changed to require updating the alt menu; else -1.
*/
altmenu_setup()
{
static char *flds[5] = {
"Data menu", "Rise/Set menu", "Separations menu"
};
int newmenu = altmenu, newhzn = alt2_stdhzn, newgeoc = alt3_geoc;
int new;
int fn = altmenu == F_MNU1 ? 0 : altmenu == F_MNU2 ? 1 : 2;
ask:
flds[3]= newhzn ? "Standard hzn" : "Adaptive hzn";
flds[4]= newgeoc? "Geocentric" : "Topocentric";
switch (popup (flds, fn, 5)) {
case 0: newmenu = F_MNU1; break;
case 1: newmenu = F_MNU2; break;
case 2: newmenu = F_MNU3; break;
case 3: newhzn ^= 1; fn = 3; goto ask;
case 4: newgeoc ^= 1; fn = 4; goto ask;
default: return (-1);
}
new = 0;
if (newmenu != altmenu) {
altmenu = newmenu;
new++;
}
if (newhzn != alt2_stdhzn) {
alt2_stdhzn = newhzn;
if (newmenu == F_MNU2)
new++;
}
if (newgeoc != alt3_geoc) {
alt3_geoc = newgeoc;
if (newmenu == F_MNU3)
new++;
}
return (new ? 0 : -1);
}
/* erase the info for the given planet */
alt_nobody (p)
int p;
{
c_pos (bodyrow[p], C_RA);
c_eol();
}
alt_body (b, force, np)
int b; /* which body, ala astro.h and screen.h defines */
int force; /* if !0 then draw for sure, else just if changed since last */
Now *np;
{
switch (altmenu) {
case F_MNU1: alt1_body (b, force, np); break;
case F_MNU2: alt2_body (b, force, np); break;
case F_MNU3: alt3_body (b, force, np); break;
}
}
/* draw the labels for the current alternate menu format */
alt_labels ()
{
switch (altmenu) {
case F_MNU1: alt1_labels (); break;
case F_MNU2: alt2_labels (); break;
case F_MNU3: alt3_labels (); break;
}
}
/* erase the labels for the current alternate menu format */
alt_nolabels ()
{
int i;
f_string (R_ALTM, C_ALTMV, " ");
for (i = R_PLANTAB; i < R_SUN; i++) {
c_pos (i, C_RA);
c_eol();
}
}
alt_menumask()
{
return (altmenu);
}
/* handy function to return the next planet in the order in which they are
* displayed in the lower half of the screen.
* input is a given planet, return is the next planet.
* if input is not legal, then first planet is returned; when input is the
* last planet, then -1 is returned.
* typical usage is something like:
* for (p = nxtbody(-1); p != -1; p = nxtbody(p))
*/
nxtbody(c)
int c;
{
static short nxtpl[NOBJ] = {
VENUS, MARS, JUPITER, SATURN, URANUS,
NEPTUNE, PLUTO, OBJX, MOON, MERCURY, -1
};
if (c < MERCURY || c >= NOBJ)
return (SUN);
else
return (nxtpl[c]);
}
static
alt1_labels()
{
f_string (R_ALTM, C_ALTMV, " Planet Data");
f_string (R_PLANTAB, C_RA, "R.A.");
f_string (R_PLANTAB+1, C_RA, "Hr:Mn.d");
f_string (R_PLANTAB, C_DEC, "Dec");
f_string (R_PLANTAB+1, C_DEC, "Deg:Mn");
f_string (R_PLANTAB, C_AZ, "Az");
f_string (R_PLANTAB+1, C_AZ, "Deg E");
f_string (R_PLANTAB, C_ALT, "Alt");
f_string (R_PLANTAB+1, C_ALT, "Deg Up");
f_string (R_PLANTAB, C_HLONG,"Helio");
f_string (R_PLANTAB+1, C_HLONG,"Long");
f_string (R_PLANTAB, C_HLAT, "Helio");
f_string (R_PLANTAB+1, C_HLAT, "Lat");
f_string (R_PLANTAB, C_EDIST,"Ea Dst");
f_string (R_PLANTAB+1, C_EDIST,"AU(mi)");
f_string (R_PLANTAB, C_SDIST,"Sn Dst");
f_string (R_PLANTAB+1, C_SDIST,"AU");
f_string (R_PLANTAB, C_ELONG,"Elong");
f_string (R_PLANTAB+1, C_ELONG,"Deg E");
f_string (R_PLANTAB, C_SIZE, "Size");
f_string (R_PLANTAB+1, C_SIZE, "ArcS");
f_string (R_PLANTAB, C_MAG, "VMag");
f_string (R_PLANTAB, C_PHASE,"Phs");
f_char (R_PLANTAB+1, C_PHASE,'%');
}
static
alt2_labels()
{
f_string (R_ALTM, C_ALTMV, "Rise/Set Info");
f_string (R_PLANTAB, C_RISETM+5, "Rise");
f_string (R_PLANTAB+1, C_RISETM+1, "Time");
f_string (R_PLANTAB+1, C_RISEAZ+2, "Az");
f_string (R_PLANTAB, C_TRANSTM+3, "Transit");
f_string (R_PLANTAB+1, C_TRANSTM+1, "Time");
f_string (R_PLANTAB+1, C_TRANSALT+2, "Alt");
f_string (R_PLANTAB, C_SETTM+5, "Set");
f_string (R_PLANTAB+1, C_SETTM+1, "Time");
f_string (R_PLANTAB+1, C_SETAZ+2, "Az");
f_string (R_PLANTAB, C_TUP, "Hrs Up");
}
static
alt3_labels()
{
f_string (R_ALTM, C_ALTMV, " Separations");
f_string (R_PLANTAB, C_SUN+2, "Sun");
f_string (R_PLANTAB, C_MOON+1, "Moon");
f_string (R_PLANTAB, C_MERCURY+1, "Merc");
f_string (R_PLANTAB, C_VENUS+1, "Venus");
f_string (R_PLANTAB, C_MARS+1, "Mars");
f_string (R_PLANTAB, C_JUPITER+2, "Jup");
f_string (R_PLANTAB, C_SATURN, "Saturn");
f_string (R_PLANTAB, C_URANUS, "Uranus");
f_string (R_PLANTAB, C_NEPTUNE+2, "Nep");
f_string (R_PLANTAB, C_PLUTO+1, "Pluto");
if (objx_ison()) {
char xnm[MAXOBJXNM+1];
char buf[10];
objx_get ((double *)0, (double *)0, (double *)0, xnm);
sprintf (buf, "%-2.2s", xnm); /* set all spaces */
f_string (R_PLANTAB, C_OBJX, buf);
}
}
/* print body info in first menu format */
static
alt1_body (p, force, np)
int p; /* which body, as in astro.h/screen.h defines */
int force; /* whether to print for sure or only if things have changed */
Now *np;
{
Sky sky;
double as = plot_ison() || srch_ison() ? 0.0 : 60.0;
int row = bodyrow[p];
if (body_cir (p, as, np, &sky) || force) {
if (p == OBJX)
pr_objxname();
f_ra (row, C_RA, sky.s_ra);
f_angle (row, C_DEC, sky.s_dec);
if (p != MOON && p != OBJX) {
f_angle (row, C_HLONG, sky.s_hlong);
f_angle (row, C_HLAT, sky.s_hlat);
}
if (p == MOON) {
/* distance is on km, show in miles */
f_double (R_MOON, C_EDIST, "%6.0f", sky.s_edist/1.609344);
} else if (p != OBJX) {
/* show distance in au */
f_double (row, C_EDIST,(sky.s_edist>=10.0)?"%6.3f":"%6.4f",
sky.s_edist);
}
if (p != OBJX)
f_double (row, C_SDIST, (sky.s_sdist>=10.0)?"%6.3f":"%6.4f",
sky.s_sdist);
f_double (row, C_ELONG, "%6.1f", sky.s_elong);
if (p != OBJX) {
f_double (row, C_SIZE, (p==MOON||p==SUN)?"%4.0f":"%4.1f",
sky.s_size);
f_double (row, C_MAG, (p==MOON||p==SUN)?"%4.0f":"%4.1f",
sky.s_mag);
if (p != SUN) {
char buf[10];
/* would just do this if Turbo-C 2.0 "%?.0f" worked:
* f_double (row, C_PHASE, "%3.0f", sky.s_phase);
*/
sprintf (buf, "%3d", (int)(sky.s_phase+0.5));
f_string (row, C_PHASE, buf);
}
}
}
f_angle (row, C_AZ, sky.s_az);
f_angle (row, C_ALT, sky.s_alt);
}
/* print body info in the second menu format */
static
alt2_body (p, force, np)
int p; /* which body, as in astro.h/screen.h defines */
int force; /* whether to print for sure or only if things have changed */
Now *np;
{
double ltr, lts, ltt, azr, azs, altt;
int row = bodyrow[p];
int status;
double tmp;
int today_tup = 0;
if (!riset_cir (p, np, alt2_stdhzn?STDHZN:ADPHZN, <r, <s, <t,
&azr, &azs, &altt, &status) && !force)
return;
alt_nobody (p);
if (p == OBJX)
pr_objxname();
if (status & RS_ERROR) {
/* can not find where body is! */
f_string (row, C_RISETM, "?Error?");
return;
}
if (status & RS_CIRCUMPOLAR) {
/* body is up all day */
f_string (row, C_RISETM, "Circumpolar");
f_mtime (row, C_TRANSTM, ltt);
if (status & RS_2TRANS)
f_char (row, C_TRANSTM+5, '+');
f_angle (row, C_TRANSALT, altt);
f_string (row, C_TUP, "24:00"); /*f_mtime() changes to 0:00 */
return;
}
if (status & RS_NEVERUP) {
/* body never up at all today */
f_string (row, C_RISETM, "Never up");
f_mtime (row, C_TUP, 0.0);
return;
}
if (status & RS_NORISE) {
/* object does not rise as such today */
f_string (row, C_RISETM, "Never rises");
ltr = 0.0; /* for TUP */
today_tup = 1;
} else {
f_mtime (row, C_RISETM, ltr);
if (status & RS_2RISES) {
/* object rises more than once today */
f_char (row, C_RISETM+5, '+');
}
f_angle (row, C_RISEAZ, azr);
}
if (status & RS_NOTRANS)
f_string (row, C_TRANSTM, "No transit");
else {
f_mtime (row, C_TRANSTM, ltt);
if (status & RS_2TRANS)
f_char (row, C_TRANSTM+5, '+');
f_angle (row, C_TRANSALT, altt);
}
if (status & RS_NOSET) {
/* object does not set as such today */
f_string (row, C_SETTM, "Never sets");
lts = 24.0; /* for TUP */
today_tup = 1;
} else {
f_mtime (row, C_SETTM, lts);
if (status & RS_2SETS)
f_char (row, C_SETTM+5, '+');
f_angle (row, C_SETAZ, azs);
}
tmp = lts - ltr;
if (tmp < 0)
tmp = 24.0 + tmp;
f_mtime (row, C_TUP, tmp);
if (today_tup)
f_char (row, C_TUP+5, '+');
}
/* print body info in third menu format. this may be either the geocentric
* or topocentric angular separation between object p and each of the others.
* the latter, of course, includes effects of refraction and so can change
* quite rapidly near the time of each planets rise or set.
* for now, we don't save old values so we always redo everything and ignore
* the "force" argument. this isn't that bad since body_cir() has memory and
* will avoid most computations as we hit them again in the lower triangle.
*/
/*ARGSUSED*/
static
alt3_body (p, force, np)
int p; /* which body, as in astro.h/screen.h defines */
int force; /* whether to print for sure or only if things have changed */
Now *np;
{
int row = bodyrow[p];
Sky skyp, skyq;
double spy, cpy, px, *qx, *qy;
int wantx = objx_ison();
double as = plot_ison() || srch_ison() ? 0.0 : 60.0;
int q;
(void) body_cir (p, as, np, &skyp);
if (alt3_geoc) {
/* use ra for "x", dec for "y". */
spy = sin (skyp.s_dec);
cpy = cos (skyp.s_dec);
px = skyp.s_ra;
qx = &skyq.s_ra;
qy = &skyq.s_dec;
} else {
/* use azimuth for "x", altitude for "y". */
spy = sin (skyp.s_alt);
cpy = cos (skyp.s_alt);
px = skyp.s_az;
qx = &skyq.s_az;
qy = &skyq.s_alt;
}
if (p == OBJX)
pr_objxname();
for (q = nxtbody(-1); q != -1; q = nxtbody(q))
if (q != p && (q != OBJX || wantx)) {
double csep;
(void) body_cir (q, as, np, &skyq);
csep = spy*sin(*qy) + cpy*cos(*qy)*cos(px-*qx);
f_angle (row, bodycol[q], acos(csep));
}
}
static
pr_objxname()
{
char n[MAXOBJXNM+1];
char buf[10];
objx_get ((double *)0, (double *)0, (double *)0, n);
sprintf (buf, "%-2.2s", n); /* set all spaces */
f_string (R_OBJX, C_OBJ, buf);
}
xXx
echo x anomaly.c
cat > anomaly.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
#define TWOPI (2*PI)
/* given the mean anomaly, ma, and the eccentricity, s, of elliptical motion,
* find the true anomaly, *nu, and the eccentric anomaly, *ea.
* all angles in radians.
*/
anomaly (ma, s, nu, ea)
double ma, s;
double *nu, *ea;
{
double m, dla, fea;
m = ma-TWOPI*(int)(ma/TWOPI);
fea = m;
while (1) {
dla = fea-(s*sin(fea))-m;
if (fabs(dla)<1e-6)
break;
dla /= 1-(s*cos(fea));
fea -= dla;
}
*nu = 2*atan(sqrt((1+s)/(1-s))*tan(fea/2));
*ea = fea;
}
xXx
echo x cal_mjd.c
cat > cal_mjd.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* given a date in months, mn, days, dy, years, yr,
* return the modified Julian date (number of days elapsed since 1900 jan 0.5),
* *mjd.
*/
cal_mjd (mn, dy, yr, mjd)
int mn, yr;
double dy;
double *mjd;
{
int b, d, m, y;
long c;
m = mn;
y = (yr < 0) ? yr + 1 : yr;
if (mn < 3) {
m += 12;
y -= 1;
}
if (yr < 1582 || yr == 1582 && (mn < 10 || mn == 10 && dy < 15))
b = 0;
else {
int a;
a = y/100;
b = 2 - a + a/4;
}
if (y < 0)
c = (long)((365.25*y) - 0.75) - 694025L;
else
c = (long)(365.25*y) - 694025L;
d = 30.6001*(m+1);
*mjd = b + c + d + dy - 0.5;
}
/* given the modified Julian date (number of days elapsed since 1900 jan 0.5,),
* mjd, return the calendar date in months, *mn, days, *dy, and years, *yr.
*/
mjd_cal (mjd, mn, dy, yr)
double mjd;
int *mn, *yr;
double *dy;
{
double d, f;
double i, a, b, ce, g;
d = mjd + 0.5;
i = floor(d);
f = d-i;
if (f == 1) {
f = 0;
i += 1;
}
if (i > -115860.0) {
a = floor((i/36524.25)+.9983573)+14;
i += 1 + a - floor(a/4.0);
}
b = floor((i/365.25)+.802601);
ce = i - floor((365.25*b)+.750001)+416;
g = floor(ce/30.6001);
*mn = g - 1;
*dy = ce - floor(30.6001*g)+f;
*yr = b + 1899;
if (g > 13.5)
*mn = g - 13;
if (*mn < 2.5)
*yr = b + 1900;
if (*yr < 1)
*yr -= 1;
}
/* given an mjd, set *dow to 0..6 according to which dayof the week it falls
* on (0=sunday) or set it to -1 if can't figure it out.
*/
mjd_dow (mjd, dow)
double mjd;
int *dow;
{
/* cal_mjd() uses Gregorian dates on or after Oct 15, 1582.
* (Pope Gregory XIII dropped 10 days, Oct 5..14, and improved the leap-
* year algorithm). however, Great Britian and the colonies did not
* adopt it until Sept 14, 1752 (they dropped 11 days, Sept 3-13,
* due to additional accumulated error). leap years before 1752 thus
* can not easily be accounted for from the cal_mjd() number...
*/
if (mjd < -53798.5) {
/* pre sept 14, 1752 too hard to correct */
*dow = -1;
return;
}
*dow = ((int)floor(mjd-.5) + 1) % 7; /* 1/1/1900 (mjd 0.5) is a Monday*/
if (*dow < 0)
*dow += 7;
}
/* given a mjd, return the the number of days in the month. */
mjd_dpm (mjd, ndays)
double mjd;
int *ndays;
{
static short dpm[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
int m, y;
double d;
mjd_cal (mjd, &m, &d, &y);
*ndays = (m==2 && ((y%4==0 && y%100!=0)||y%400==0)) ? 29 : dpm[m-1];
}
/* given a mjd, return the year as a double. */
mjd_year (mjd, yr)
double mjd;
double *yr;
{
int m, y;
double d;
double e0, e1; /* mjd of start of this year, start of next year */
mjd_cal (mjd, &m, &d, &y);
cal_mjd (1, 1.0, y, &e0);
cal_mjd (12, 31.0, y, &e1); e1 += 1.0;
*yr = y + (mjd - e0)/(e1 - e0);
}
xXx
echo x circum.c
cat > circum.c << 'xXx'
/* fill in a Sky struct with all we know about each object.
*(object-x is in objx.c)
*/
#include <stdio.h>
#include <math.h>
#include "astro.h"
#include "circum.h"
#include "screen.h" /* just for SUN and MOON */
/* find body p's circumstances now.
* to save some time the caller may specify a desired accuracy, in arc seconds.
* if, based on its mean motion, it would not have moved this much since the
* last time we were called we only recompute altitude and azimuth and avoid
* recomputing the planet's heliocentric position. use 0.0 for best possible.
* return 0 if only alt/az changes, else 1 if all other stuff updated too.
* TODO: beware of fast retrograde motions.
*/
body_cir (p, as, np, sp)
int p;
double as;
Now *np;
Sky *sp;
{
typedef struct {
double l_dpas; /* mean days per arc second */
Now l_now; /* when l_sky was found */
Sky l_sky;
} Last;
/* must be in same order as the astro.h #define's */
static Last last[8] =
{{.000068},{.00017},{.00053},{.0034},{.0092},{.027},{.046},{.069}};
double lst, alt, az;
Last *lp;
int new;
switch (p) {
case SUN:
return (sun_cir (as, np, sp));
case MOON:
return (moon_cir (as, np, sp));
case OBJX:
return (objx_cir (as, np, sp));
}
/* if less than l_every days from last time for this planet
* just redo alt/az.
*/
lp = last + p;
if (same_cir(np,&lp->l_now) && about_now(np,&lp->l_now,as*lp->l_dpas)) {
*sp = lp->l_sky;
new = 0;
} else {
double lpd0, psi0; /* heliocentric ecliptic longitude and latitude */
double rp0; /* dist from sun */
double rho0; /* dist from earth */
double lam, bet; /* geocentric ecliptic long and lat */
double dia, mag; /* angular diameter at 1 AU and magnitude */
double lsn, rsn; /* true geoc lng of sun, dist from sn to earth*/
double deps, dpsi;
double a, ca, sa;
double el; /* elongation */
double f; /* phase from earth */
lp->l_now = *np;
plans (mjd, p, &lpd0, &psi0, &rp0, &rho0, &lam, &bet, &dia, &mag);
nutation (mjd, &deps, &dpsi); /* correct for nutation */
lam += dpsi;
sunpos (mjd, &lsn, &rsn);
/* correct for 20.4" aberration */
a = lsn-lam;
ca = cos(a);
sa = sin(a);
lam -= degrad(20.4/3600)*ca/cos(bet);
bet -= degrad(20.4/3600)*sa*sin(bet);
ecl_eq (mjd, bet, lam, &sp->s_ra, &sp->s_dec);
if (epoch != EOD)
precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
sp->s_edist = rho0;
sp->s_sdist = rp0;
elongation (lam, bet, lsn, &el);
el = raddeg(el);
sp->s_elong = el;
sp->s_size = dia/rho0;
f = 0.5*(1+cos(lam-lpd0));
sp->s_phase = f*100.0; /* percent */
sp->s_mag = 5.0*log(rp0*rho0/sqrt(f))/log(10.0) + mag;
sp->s_hlong = lpd0;
sp->s_hlat = psi0;
new = 1;
}
/* alt, az; correct for refraction, in place */
now_lst (np, &lst);
hadec_aa (lat, hrrad(lst) - sp->s_ra, sp->s_dec, &alt, &az);
refract (pressure, temp, alt, &alt);
sp->s_alt = alt;
sp->s_az = az;
lp->l_sky = *sp;
return (new);
}
/* find local times when sun is 18 degrees below horizon.
* return 0 if just returned same stuff as previous call, else 1 if new.
*/
twilight_cir (np, dawn, dusk, status)
Now *np;
double *dawn, *dusk;
int *status;
{
static Now last_now;
static double last_dawn, last_dusk;
static int last_status;
int new;
if (same_cir (np, &last_now) && same_lday (np, &last_now)) {
*dawn = last_dawn;
*dusk = last_dusk;
*status = last_status;
new = 0;
} else {
double x;
(void) riset_cir (SUN,np,TWILIGHT,dawn,dusk,&x,&x,&x,&x,status);
last_dawn = *dawn;
last_dusk = *dusk;
last_status = *status;
last_now = *np;
new = 1;
}
return (new);
}
/* find sun's circumstances now.
* as is the desired accuracy, in arc seconds; use 0.0 for best possible.
* return 0 if only alt/az changes, else 1 if all other stuff updated too.
*/
sun_cir (as, np, sp)
double as;
Now *np;
Sky *sp;
{
static Sky last_sky;
static Now last_now;
double lst, alt, az;
int new;
if (same_cir (np, &last_now) && about_now (np, &last_now, as*.00028)) {
*sp = last_sky;
new = 0;
} else {
double lsn, rsn;
double deps, dpsi;
last_now = *np;
sunpos (mjd, &lsn, &rsn); /* sun's true ecliptic long
* and dist
*/
nutation (mjd, &deps, &dpsi); /* correct for nutation */
lsn += dpsi-degrad(20.4/3600); /* and 20.4" aberration */
sp->s_edist = rsn;
sp->s_sdist = 0.0;
sp->s_elong = 0.0;
sp->s_size = raddeg(4.65242e-3/rsn)*3600*2;
sp->s_mag = -26.8;
sp->s_hlong = lsn-PI; /* geo- to helio- centric */
range (&sp->s_hlong, 2*PI);
sp->s_hlat = 0.0;
ecl_eq (mjd, 0.0, lsn, &sp->s_ra, &sp->s_dec);
if (epoch != EOD)
precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
new = 1;
}
now_lst (np, &lst);
hadec_aa (lat, hrrad(lst) - sp->s_ra, sp->s_dec, &alt, &az);
refract (pressure, temp, alt, &alt);
sp->s_alt = alt;
sp->s_az = az;
last_sky = *sp;
return (new);
}
/* find moon's circumstances now.
* as is the desired accuracy, in arc seconds; use 0.0 for best possible.
* return 0 if only alt/az changes, else 1 if all other stuff updated too.
*/
moon_cir (as, np, sp)
double as;
Now *np;
Sky *sp;
{
static Sky last_sky;
static Now last_now;
static double ehp;
double lst, alt, az;
double ha, dec;
int new;
if (same_cir (np, &last_now) && about_now (np, &last_now, as*.000021)) {
*sp = last_sky;
new = 0;
} else {
double lam, bet;
double deps, dpsi;
double lsn, rsn; /* sun long in rads, earth-sun dist in au */
double edistau; /* earth-moon dist, in au */
double el; /* elongation, rads east */
last_now = *np;
moon (mjd, &lam, &bet, &ehp); /* moon's true ecliptic loc */
nutation (mjd, &deps, &dpsi); /* correct for nutation */
lam += dpsi;
range (&lam, 2*PI);
sp->s_edist = 6378.14/sin(ehp); /* earth-moon dist, want km */
sp->s_size = 3600*31.22512*sin(ehp);/* moon angular dia, seconds */
ecl_eq (mjd, bet, lam, &sp->s_ra, &sp->s_dec);
if (epoch != EOD)
precess (mjd, epoch, &sp->s_ra, &sp->s_dec);
sunpos (mjd, &lsn, &rsn);
range (&lsn, 2*PI);
elongation (lam, bet, lsn, &el);
/* solve triangle of earth, sun, and elongation for moon-sun dist */
edistau = sp->s_edist/1.495979e8; /* km -> au */
sp->s_sdist =
sqrt (edistau*edistau + rsn*rsn - 2.0*edistau*rsn*cos(el));
/* TODO: improve mag; this is based on a flat moon model. */
sp->s_mag = -12 + 2.5*(log10(PI) - log10(PI/2*(1+1.e-6-cos(el))));
sp->s_elong = raddeg(el); /* want degrees */
sp->s_phase = fabs(el)/PI*100.0; /* want non-negative % */
sp->s_hlong = sp->s_hlat = 0.0;
new = 1;
}
/* show topocentric alt/az by correcting ra/dec for parallax
* as well as refraction.
*/
now_lst (np, &lst);
ha = hrrad(lst) - sp->s_ra;
ta_par (ha, sp->s_dec, lat, height, ehp, &ha, &dec);
hadec_aa (lat, ha, dec, &alt, &az);
refract (pressure, temp, alt, &alt);
sp->s_alt = alt;
sp->s_az = az;
last_sky = *sp;
return (new);
}
/* given geocentric ecliptic longitude and latitude, lam and bet, of some object
* and the longitude of the sun, lsn, find the elongation, el. this is the
* actual angular separation of the object from the sun, not just the difference
* in the longitude. the sign, however, IS set simply as a test on longitude
* such that el will be >0 for an evening object <0 for a morning object.
* to understand the test for el sign, draw a graph with lam going from 0-2*PI
* down the vertical axis, lsn going from 0-2*PI across the hor axis. then
* define the diagonal regions bounded by the lines lam=lsn+PI, lam=lsn and
* lam=lsn-PI. the "morning" regions are any values to the lower left of the
* first line and bounded within the second pair of lines.
* all angles in radians.
*/
elongation (lam, bet, lsn, el)
double lam, bet, lsn;
double *el;
{
*el = acos(cos(bet)*cos(lam-lsn));
if (lam>lsn+PI || lam>lsn-PI && lam<lsn) *el = - *el;
}
/* return whether the two Nows are for the same observing circumstances. */
same_cir (n1, n2)
register Now *n1, *n2;
{
return (n1->n_lat == n2->n_lat
&& n1->n_lng == n2->n_lng
&& n1->n_temp == n2->n_temp
&& n1->n_pressure == n2->n_pressure
&& n1->n_height == n2->n_height
&& n1->n_tz == n2->n_tz
&& n1->n_epoch == n2->n_epoch);
}
/* return whether the two Nows are for the same LOCAL day */
same_lday (n1, n2)
Now *n1, *n2;
{
return (mjd_day(n1->n_mjd - n1->n_tz/24.0) ==
mjd_day(n2->n_mjd - n2->n_tz/24.0));
}
/* return whether the mjd of the two Nows are within dt */
static
about_now (n1, n2, dt)
Now *n1, *n2;
double dt;
{
return (fabs (n1->n_mjd - n2->n_mjd) <= dt/2.0);
}
now_lst (np, lst)
Now *np;
double *lst;
{
utc_gst (mjd_day(mjd), mjd_hr(mjd), lst);
*lst += radhr(lng);
range (lst, 24.0);
}
/* round a time in days, *t, to the nearest second, IN PLACE. */
rnd_second (t)
double *t;
{
*t = floor(*t*SPD+0.5)/SPD;
}
double mjd_day(jd)
double jd;
{
return (floor(jd-0.5)+0.5);
}
double mjd_hr(jd)
double jd;
{
return ((jd-mjd_day(jd))*24.0);
}
xXx
echo x compiler.c
cat > compiler.c << 'xXx'
/* module to compile and execute a c-style arithmetic expression.
* public entry points are compile_expr() and execute_expr().
*
* one reason this is so nice and tight is that all opcodes are the same size
* (an int) and the tokens the parser returns are directly usable as opcodes,
* for the most part. constants and variables are compiled as an opcode
* with an offset into the auxiliary opcode tape, opx.
*/
#include <math.h>
#include "screen.h"
/* parser tokens and opcodes, as necessary */
#define HALT 0 /* good value for HALT since program is inited to 0 */
/* binary operators (precedences in table, below) */
#define ADD 1
#define SUB 2
#define MULT 3
#define DIV 4
#define AND 5
#define OR 6
#define GT 7
#define GE 8
#define EQ 9
#define NE 10
#define LT 11
#define LE 12
/* unary op, precedence in NEG_PREC #define, below */
#define NEG 13
/* symantically operands, ie, constants, variables and all functions */
#define CONST 14
#define VAR 15
#define ABS 16 /* add functions if desired just like this is done */
/* purely tokens - never get compiled as such */
#define LPAREN 255
#define RPAREN 254
#define ERR (-1)
/* precedence of each of the binary operators.
* in case of a tie, compiler associates left-to-right.
* N.B. each entry's index must correspond to its #define!
*/
static int precedence[] = {0,5,5,6,6,2,1,4,4,3,3,4,4};
#define NEG_PREC 7 /* negation is highest */
/* execute-time operand stack */
#define MAX_STACK 16
static double stack[MAX_STACK], *sp;
/* space for compiled opcodes - the "program".
* opcodes go in lower 8 bits.
* when an opcode has an operand (as CONST and VAR) it is really in opx[] and
* the index is in the remaining upper bits.
*/
#define MAX_PROG 32
static int program[MAX_PROG], *pc;
#define OP_SHIFT 8
#define OP_MASK 0xff
/* auxiliary operand info.
* the operands (all but lower 8 bits) of CONST and VAR are really indeces
* into this array. thus, no point in making this any longer than you have
* bits more than 8 in your machine's int to index into it, ie, make
* MAX_OPX <= 1 << ((sizeof(int)-1)*8)
* also, the fld's must refer to ones being flog'd, so not point in more
* of these then that might be used for plotting and srching combined.
*/
#define MAX_OPX 16
typedef union {
double opu_f; /* value when opcode is CONST */
int opu_fld; /* rcfpack() of field when opcode is VAR */
} OpX;
static OpX opx[MAX_OPX];
static int opxidx;
/* these are global just for easy/rapid access */
static int parens_nest; /* to check that parens end up nested */
static char *err_msg; /* caller provides storage; we point at it with this */
static char *cexpr, *lcexpr; /* pointers that move along caller's expression */
static int good_prog; /* != 0 when program appears to be good */
/* compile the given c-style expression.
* return 0 and set good_prog if ok,
* else return -1 and a reason message in errbuf.
*/
compile_expr (ex, errbuf)
char *ex;
char *errbuf;
{
int instr;
/* init the globals.
* also delete any flogs used in the previous program.
*/
cexpr = ex;
err_msg = errbuf;
pc = program;
opxidx = 0;
parens_nest = 0;
do {
instr = *pc++;
if ((instr & OP_MASK) == VAR)
flog_delete (opx[instr >> OP_SHIFT].opu_fld);
} while (instr != HALT);
pc = program;
if (compile(0) == ERR) {
sprintf (err_msg + strlen(err_msg), " at \"%.10s\"", lcexpr);
good_prog = 0;
return (-1);
}
*pc++ = HALT;
good_prog = 1;
return (0);
}
/* execute the expression previously compiled with compile_expr().
* return 0 with *vp set to the answer if ok, else return -1 with a reason
* why not message in errbuf.
*/
execute_expr (vp, errbuf)
double *vp;
char *errbuf;
{
int s;
err_msg = errbuf;
sp = stack + MAX_STACK; /* grows towards lower addresses */
pc = program;
s = execute(vp);
if (s < 0)
good_prog = 0;
return (s);
}
/* this is a way for the outside world to ask whether there is currently a
* reasonable program compiled and able to execute.
*/
prog_isgood()
{
return (good_prog);
}
/* get and return the opcode corresponding to the next token.
* leave with lcexpr pointing at the new token, cexpr just after it.
* also watch for mismatches parens and proper operator/operand alternation.
*/
static
next_token ()
{
static char toomt[] = "More than %d terms";
static char badop[] = "Illegal operator";
int tok = ERR; /* just something illegal */
char c;
while ((c = *cexpr) == ' ')
cexpr++;
lcexpr = cexpr++;
/* mainly check for a binary operator */
switch (c) {
case '\0': --cexpr; tok = HALT; break; /* keep returning HALT */
case '+': tok = ADD; break; /* compiler knows when it's really unary */
case '-': tok = SUB; break; /* compiler knows when it's really negate */
case '*': tok = MULT; break;
case '/': tok = DIV; break;
case '(': parens_nest++; tok = LPAREN; break;
case ')':
if (--parens_nest < 0) {
sprintf (err_msg, "Too many right parens");
return (ERR);
} else
tok = RPAREN;
break;
case '|':
if (*cexpr == '|') { cexpr++; tok = OR; }
else { sprintf (err_msg, badop); return (ERR); }
break;
case '&':
if (*cexpr == '&') { cexpr++; tok = AND; }
else { sprintf (err_msg, badop); return (ERR); }
break;
case '=':
if (*cexpr == '=') { cexpr++; tok = EQ; }
else { sprintf (err_msg, badop); return (ERR); }
break;
case '!':
if (*cexpr == '=') { cexpr++; tok = NE; }
else { sprintf (err_msg, badop); return (ERR); }
break;
case '<':
if (*cexpr == '=') { cexpr++; tok = LE; }
else tok = LT;
break;
case '>':
if (*cexpr == '=') { cexpr++; tok = GE; }
else tok = GT;
break;
}
if (tok != ERR)
return (tok);
/* not op so check for a constant, variable or function */
if (isdigit(c) || c == '.') {
if (opxidx > MAX_OPX) {
sprintf (err_msg, toomt, MAX_OPX);
return (ERR);
}
opx[opxidx].opu_f = atof (lcexpr);
tok = CONST | (opxidx++ << OP_SHIFT);
skip_double();
} else if (isalpha(c)) {
/* check list of functions */
if (strncmp (lcexpr, "abs", 3) == 0) {
cexpr += 2;
tok = ABS;
} else {
/* not a function, so assume it's a variable */
int fld;
if (opxidx > MAX_OPX) {
sprintf (err_msg, toomt, MAX_OPX);
return (ERR);
}
fld = parse_fieldname ();
if (fld < 0) {
sprintf (err_msg, "Unknown field");
return (ERR);
} else {
if (flog_add (fld) < 0) { /* register with field logger */
sprintf (err_msg, "Sorry; too many fields");
return (ERR);
}
opx[opxidx].opu_fld = fld;
tok = VAR | (opxidx++ << OP_SHIFT);
}
}
}
return (tok);
}
/* move cexpr on past a double.
* allow sci notation.
* no need to worry about a leading '-' or '+' but allow them after an 'e'.
* TODO: this handles all the desired cases, but also admits a bit too much
* such as things like 1eee2...3. geeze; to skip a double right you almost
* have to go ahead and crack it!
*/
static
skip_double()
{
int sawe = 0; /* so we can allow '-' or '+' right after an 'e' */
while (1) {
char c = *cexpr;
if (isdigit(c) || c=='.' || (sawe && (c=='-' || c=='+'))) {
sawe = 0;
cexpr++;
} else if (c == 'e') {
sawe = 1;
cexpr++;
} else
break;
}
}
/* call this whenever you want to dig out the next (sub)expression.
* keep compiling instructions as long as the operators are higher precedence
* than prec, then return that "look-ahead" token that wasn't (higher prec).
* if error, fill in a message in err_msg[] and return ERR.
*/
static
compile (prec)
int prec;
{
int expect_binop = 0; /* set after we have seen any operand.
* used by SUB so it can tell if it really
* should be taken to be a NEG instead.
*/
int tok = next_token ();
while (1) {
int p;
if (tok == ERR)
return (ERR);
if (pc - program >= MAX_PROG) {
sprintf (err_msg, "Program is too long");
return (ERR);
}
/* check for special things like functions, constants and parens */
switch (tok & OP_MASK) {
case HALT: return (tok);
case ADD:
if (expect_binop)
break; /* procede with binary addition */
/* just skip a unary positive(?) */
tok = next_token();
continue;
case SUB:
if (expect_binop)
break; /* procede with binary subtract */
tok = compile (NEG_PREC);
*pc++ = NEG;
expect_binop = 1;
continue;
case ABS: /* other funcs would be handled the same too ... */
/* eat up the function parenthesized argument */
if (next_token() != LPAREN || compile (0) != RPAREN) {
sprintf (err_msg, "Function arglist error");
return (ERR);
}
/* then handled same as ... */
case CONST: /* handled same as... */
case VAR:
*pc++ = tok;
tok = next_token();
expect_binop = 1;
continue;
case LPAREN:
if (compile (0) != RPAREN) {
sprintf (err_msg, "Unmatched left paren");
return (ERR);
}
tok = next_token();
expect_binop = 1;
continue;
case RPAREN:
return (RPAREN);
}
/* everything else is a binary operator */
p = precedence[tok];
if (p > prec) {
int newtok = compile (p);
if (newtok == ERR)
return (ERR);
*pc++ = tok;
expect_binop = 1;
tok = newtok;
} else
return (tok);
}
}
/* "run" the program[] compiled with compile().
* if ok, return 0 and the final result,
* else return -1 with a reason why not message in err_msg.
*/
static
execute(result)
double *result;
{
int instr;
do {
instr = *pc++;
switch (instr & OP_MASK) {
/* put these in numberic order so hopefully even the dumbest
* compiler will choose to use a jump table, not a cascade of ifs.
*/
case HALT: break; /* outer loop will stop us */
case ADD: sp[1] = sp[1] + sp[0]; sp++; break;
case SUB: sp[1] = sp[1] - sp[0]; sp++; break;
case MULT: sp[1] = sp[1] * sp[0]; sp++; break;
case DIV: sp[1] = sp[1] / sp[0]; sp++; break;
case AND: sp[1] = sp[1] && sp[0] ? 1 : 0; sp++; break;
case OR: sp[1] = sp[1] || sp[0] ? 1 : 0; sp++; break;
case GT: sp[1] = sp[1] > sp[0] ? 1 : 0; sp++; break;
case GE: sp[1] = sp[1] >= sp[0] ? 1 : 0; sp++; break;
case EQ: sp[1] = sp[1] == sp[0] ? 1 : 0; sp++; break;
case NE: sp[1] = sp[1] != sp[0] ? 1 : 0; sp++; break;
case LT: sp[1] = sp[1] < sp[0] ? 1 : 0; sp++; break;
case LE: sp[1] = sp[1] <= sp[0] ? 1 : 0; sp++; break;
case NEG: *sp = -*sp; break;
case CONST: *--sp = opx[instr >> OP_SHIFT].opu_f; break;
case VAR:
if (flog_get (opx[instr >> OP_SHIFT].opu_fld, --sp) < 0) {
sprintf (err_msg, "Bug! VAR field not logged");
return (-1);
}
break;
case ABS: *sp = fabs (*sp); break;
default:
sprintf (err_msg, "Bug! bad opcode: 0x%x", instr);
return (-1);
}
if (sp < stack) {
sprintf (err_msg, "Runtime stack overflow");
return (-1);
} else if (sp - stack > MAX_STACK) {
sprintf (err_msg, "Bug! runtime stack underflow");
return (-1);
}
} while (instr != HALT);
/* result should now be on top of stack */
if (sp != &stack[MAX_STACK - 1]) {
sprintf (err_msg, "Bug! stack has %d items",MAX_STACK-(sp-stack));
return (-1);
}
*result = *sp;
return (0);
}
static
isdigit(c)
char c;
{
return (c >= '0' && c <= '9');
}
static
isalpha (c)
char c;
{
return ((c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z'));
}
/* starting with lcexpr pointing at a string expected to be a field name,
* return an rcfpack(r,c,0) of the field else -1 if bad.
* when return, leave lcexpr alone but move cexpr to just after the name.
*/
static
parse_fieldname ()
{
int r = -1, c = -1; /* anything illegal */
char *fn = lcexpr; /* likely faster than using the global */
char f0, f1;
char *dp;
/* search for first thing not an alpha char.
* leave it in f0 and leave dp pointing to it.
*/
dp = fn;
while (isalpha(f0 = *dp))
dp++;
/* crack the new field name.
* when done trying, leave dp pointing at first char just after it.
* set r and c if we recognized it.
*/
if (f0 == '.') {
/* planet.column pair.
* first crack the planet portion (pointed to by fn): set r.
* then the second portion (pointed to by dp+1): set c.
*/
char xn[MAXOBJXNM+1];
f0 = fn[0];
f1 = fn[1];
/* if there is an object-x we insist on a perfect match */
objx_get ((double *)0, (double *)0, (double *)0, xn);
if (xn[0] && f0 == xn[0] && (!xn[1] || f1 == xn[1]))
r = R_OBJX;
else
switch (f0) {
case 'j':
r = R_JUPITER;
break;
case 'm':
if (f1 == 'a') r = R_MARS;
else if (f1 == 'e') r = R_MERCURY;
else if (f1 == 'o') r = R_MOON;
break;
case 'n':
r = R_NEPTUNE;
break;
case 'p':
r = R_PLUTO;
break;
case 's':
if (f1 == 'a') r = R_SATURN;
else if (f1 == 'u') r = R_SUN;
break;
case 'u':
r = R_URANUS;
break;
case 'v':
r = R_VENUS;
break;
}
/* now crack the column (stuff after the dp) */
dp++; /* point at good stuff just after the decimal pt */
f0 = dp[0];
f1 = dp[1];
switch (f0) {
case 'a':
if (f1 == 'l') c = C_ALT;
else if (f1 == 'z') c = C_AZ;
break;
case 'd':
c = C_DEC;
break;
case 'e':
if (f1 == 'd') c = C_EDIST;
else if (f1 == 'l') c = C_ELONG;
break;
case 'h':
if (f1 == 'l') {
if (dp[2] == 'a') c = C_HLAT;
else if (dp[2] == 'o') c = C_HLONG;
} else if (f1 == 'r' || f1 == 'u') c = C_TUP;
break;
case 'j':
c = C_JUPITER;
break;
case 'm':
if (f1 == 'a') c = C_MARS;
else if (f1 == 'e') c = C_MERCURY;
else if (f1 == 'o') c = C_MOON;
break;
case 'n':
c = C_NEPTUNE;
break;
case 'p':
if (f1 == 'h') c = C_PHASE;
else if (f1 == 'l') c = C_PLUTO;
break;
case 'r':
if (f1 == 'a') {
if (dp[2] == 'z') c = C_RISEAZ;
else c = C_RA;
} else if (f1 == 't') c = C_RISETM;
break;
case 's':
if (f1 == 'a') {
if (dp[2] == 'z') c = C_SETAZ;
else c = C_SATURN;
} else if (f1 == 'd') c = C_SDIST;
else if (f1 == 'i') c = C_SIZE;
else if (f1 == 't') c = C_SETTM;
else if (f1 == 'u') c = C_SUN;
break;
case 't':
if (f1 == 'a') c = C_TRANSALT;
else if (f1 == 't') c = C_TRANSTM;
break;
case 'u':
c = C_URANUS;
break;
case 'v':
if (f1 == 'e') c = C_VENUS;
else if (f1 == 'm') c = C_MAG;
break;
}
/* now skip dp on past the column stuff */
while (isalpha(*dp))
dp++;
} else {
/* no decimal point; some field in the top of the screen */
f0 = fn[0];
f1 = fn[1];
switch (f0) {
case 'd':
if (f1 == 'a') r = R_DAWN, c = C_DAWN;
else if (f1 == 'u') r = R_DUSK, c = C_DUSK;
break;
case 'n':
r = R_LON, c = C_LON;
break;
}
}
cexpr = dp;
if (r <= 0 || c <= 0) return (-1);
return (rcfpack (r, c, 0));
}
xXx
echo x eq_ecl.c
cat > eq_ecl.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
#define EQtoECL 1
#define ECLtoEQ (-1)
/* given the modified Julian date, mjd, and an equitorial ra and dec, each in
* radians, find the corresponding geocentric ecliptic latitude, *lat, and
* longititude, *lng, also each in radians.
* the effects of nutation and the angle of the obliquity are included.
*/
eq_ecl (mjd, ra, dec, lat, lng)
double mjd, ra, dec;
double *lat, *lng;
{
ecleq_aux (EQtoECL, mjd, ra, dec, lng, lat);
}
/* given the modified Julian date, mjd, and a geocentric ecliptic latitude,
* *lat, and longititude, *lng, each in radians, find the corresponding
* equitorial ra and dec, also each in radians.
* the effects of nutation and the angle of the obliquity are included.
*/
ecl_eq (mjd, lat, lng, ra, dec)
double mjd, lat, lng;
double *ra, *dec;
{
ecleq_aux (ECLtoEQ, mjd, lng, lat, ra, dec);
}
static
ecleq_aux (sw, mjd, x, y, p, q)
int sw; /* +1 for eq to ecliptic, -1 for vv. */
double mjd, x, y; /* sw==1: x==ra, y==dec. sw==-1: x==lng, y==lat. */
double *p, *q; /* sw==1: p==lng, q==lat. sw==-1: p==ra, q==dec. */
{
static double lastmjd; /* last mjd calculated */
static double seps, ceps; /* sin and cos of mean obliquity */
double sx, cx, sy, cy, ty;
if (mjd != lastmjd) {
double deps, dpsi, eps;
obliquity (mjd, &eps); /* mean obliquity for date */
#ifdef NONUTATION
#define NONUTATION
deps = 0; /* checkout handler did not
* include nutation correction.
*/
#else
nutation (mjd, &deps, &dpsi); /* correctin due to nutation */
#endif
eps += deps; /* true obliquity for date */
seps = sin(eps);
ceps = cos(eps);
lastmjd = mjd;
}
sy = sin(y);
cy = cos(y); /* always non-negative */
if (fabs(cy)<1e-20) cy = 1e-20; /* insure > 0 */
ty = sy/cy;
cx = cos(x);
sx = sin(x);
*q = asin((sy*ceps)-(cy*seps*sx*sw));
*p = atan(((sx*ceps)+(ty*seps*sw))/cx);
if (cx<0) *p += PI; /* account for atan quad ambiguity */
range (p, 2*PI);
}
xXx
echo x flog.c
cat > flog.c << 'xXx'
/* this is a simple little package to manage the saving and retrieving of
* field values, which we call field logging or "flogs". a flog consists of a
* field location, ala rcfpack(), and its value as a double. you can reset the
* list of flogs, add to and remove from the list of registered fields and log
* a field if it has been registered.
*
* this is used by the plotting and searching facilities of ephem to maintain
* the values of the fields that are being plotted or used in search
* expressions.
*
* a field can be in use for more than one
* thing at a time (eg, all the X plot values may the same time field, or
* searching and plotting might be on at one time using the same field) so
* we consider the field to be in use as long a usage count is > 0.
*/
#include "screen.h"
#define NFLOGS 32
typedef struct {
int fl_usagecnt; /* number of "users" logging to this field */
int fl_fld; /* an rcfpack(r,c,0) */
double fl_val;
} FLog;
static FLog flog[NFLOGS];
/* add fld to the list. if already there, just increment usage count.
* return 0 if ok, else -1 if no more room.
*/
flog_add (fld)
int fld;
{
FLog *flp, *unusedflp = 0;
/* scan for fld already in list, or find an unused one along the way */
for (flp = &flog[NFLOGS]; --flp >= flog; ) {
if (flp->fl_usagecnt > 0) {
if (flp->fl_fld == fld) {
flp->fl_usagecnt++;
return (0);
}
} else
unusedflp = flp;
}
if (unusedflp) {
unusedflp->fl_fld = fld;
unusedflp->fl_usagecnt = 1;
return (0);
}
return (-1);
}
/* decrement usage count for flog for fld. if goes to 0 take it out of list.
* ok if not in list i guess...
*/
flog_delete (fld)
int fld;
{
FLog *flp;
for (flp = &flog[NFLOGS]; --flp >= flog; )
if (flp->fl_fld == fld && flp->fl_usagecnt > 0) {
if (--flp->fl_usagecnt <= 0) {
flp->fl_usagecnt = 0;
}
break;
}
}
/* if plotting or searching is active then
* if rcfpack(r,c,0) is in the fld list, set its value to val.
* return 0 if ok, else -1 if not in list.
*/
flog_log (r, c, val)
int r, c;
double val;
{
if (plot_ison() || srch_ison()) {
FLog *flp;
int fld = rcfpack (r, c, 0);
for (flp = &flog[NFLOGS]; --flp >= flog; )
if (flp->fl_fld == fld && flp->fl_usagecnt > 0) {
flp->fl_val = val;
return(0);
}
return (-1);
} else
return (0);
}
/* search for fld in list. if find it return its value.
* return 0 if found it, else -1 if not in list.
*/
flog_get (fld, vp)
int fld;
double *vp;
{
FLog *flp;
for (flp = &flog[NFLOGS]; --flp >= flog; )
if (flp->fl_fld == fld && flp->fl_usagecnt > 0) {
*vp = flp->fl_val;
return (0);
}
return (-1);
}
xXx
echo x formats.c
cat > formats.c << 'xXx'
/* basic formating routines.
* all the screen oriented printing should go through here.
*/
#include <stdio.h>
#include <math.h>
#include "astro.h"
#include "screen.h"
/* suppress screen io if this is true, but always flog stuff.
*/
static int f_scrnoff;
f_on ()
{
f_scrnoff = 0;
}
f_off ()
{
f_scrnoff = 1;
}
/* draw n blanks at the given cursor position. */
f_blanks (r, c, n)
int r, c, n;
{
if (f_scrnoff)
return;
c_pos (r, c);
while (--n >= 0)
putchar (' ');
}
/* print the given value, v, in "sexadecimal" format at [r,c]
* ie, in the form A:m.P, where A is a digits wide, P is p digits.
* if p == 0, then no decimal point either.
*/
f_sexad (r, c, a, p, mod, v)
int r, c;
int a, p; /* left space, min precision */
int mod; /* don't let whole portion get this big */
double v;
{
char astr[32], str[32];
int dec;
double frac;
int visneg;
(void) flog_log (r, c, v);
if (f_scrnoff)
return;
if (v >= 0.0)
visneg = 0;
else {
visneg = 1;
v = -v;
}
dec = v;
frac = (v - dec)*60.0;
sprintf (str, "59.%.*s5", p, "999999999");
if (frac >= atof (str)) {
dec += 1;
frac = 0.0;
}
dec %= mod;
if (dec == 0 && visneg)
strcpy (str, "-0");
else
sprintf (str, "%d", visneg ? -dec : dec);
/* would just do this if Turbo-C 2.0 %?.0f" worked:
* sprintf (astr, "%*s:%0*.*f", a, str, p == 0 ? 2 : p+3, p, frac);
*/
if (p == 0)
sprintf (astr, "%*s:%02d", a, str, (int)(frac+0.5));
else
sprintf (astr, "%*s:%0*.*f", a, str, p+3, p, frac);
f_string (r, c, astr);
}
/* print the given value, t, in sexagesimal format at [r,c]
* ie, in the form T:mm:ss, where T is nd digits wide.
* N.B. we assume nd >= 2.
*/
f_sexag (r, c, nd, t)
int r, c, nd;
double t;
{
char tstr[32];
int h, m, s;
int tisneg;
(void) flog_log (r, c, t);
if (f_scrnoff)
return;
dec_sex (t, &h, &m, &s, &tisneg);
if (h == 0 && tisneg)
sprintf (tstr, "%*s-0:%02d:%02d", nd-2, "", m, s);
else
sprintf (tstr, "%*d:%02d:%02d", nd, tisneg ? -h : h, m, s);
f_string (r, c, tstr);
}
/* print angle ra, in radians, in ra hours as hh:mm.m at [r,c]
* N.B. we assume ra is >= 0.
*/
f_ra (r, c, ra)
int r, c;
double ra;
{
f_sexad (r, c, 2, 1, 24, radhr(ra));
}
/* print time, t, as hh:mm:ss */
f_time (r, c, t)
int r, c;
double t;
{
f_sexag (r, c, 2, t);
}
/* print time, t, as +/-hh:mm:ss (don't show leading +) */
f_signtime (r, c, t)
int r, c;
double t;
{
f_sexag (r, c, 3, t);
}
/* print time, t, as hh:mm */
f_mtime (r, c, t)
int r, c;
double t;
{
f_sexad (r, c, 2, 0, 24, t);
}
/* print angle, a, in rads, as degress at [r,c] in form ddd:mm */
f_angle(r, c, a)
int r, c;
double a;
{
f_sexad (r, c, 3, 0, 360, raddeg(a));
}
/* print angle, a, in rads, as degress at [r,c] in form dddd:mm:ss */
f_gangle(r, c, a)
int r, c;
double a;
{
f_sexag (r, c, 4, raddeg(a));
}
/* print the given modified Julian date, jd, as the starting date at [r,c]
* in the form mm/dd/yyyy.
*/
f_date (r, c, jd)
int r, c;
double jd;
{
char dstr[32];
int m, y;
double d, tmp;
/* shadow to the plot subsystem as years. */
mjd_year (jd, &tmp);
(void) flog_log (r, c, tmp);
if (f_scrnoff)
return;
mjd_cal (jd, &m, &d, &y);
sprintf (dstr, "%2d/%02d/%04d", m, (int)(d), y);
f_string (r, c, dstr);
}
f_char (row, col, c)
int row, col;
char c;
{
if (f_scrnoff)
return;
c_pos (row, col);
putchar (c);
}
f_string (r, c, s)
int r, c;
char *s;
{
if (f_scrnoff)
return;
c_pos (r, c);
fputs (s, stdout);
}
f_double (r, c, fmt, f)
int r, c;
char *fmt;
double f;
{
char str[80];
(void) flog_log (r, c, f);
sprintf (str, fmt, f);
f_string (r, c, str);
}
/* print prompt line */
f_prompt (p)
char *p;
{
c_pos (R_PROMPT, C_PROMPT);
c_eol ();
c_pos (R_PROMPT, C_PROMPT);
fputs (p, stdout);
}
/* print a message and wait for op to hit any key */
f_msg (m)
char *m;
{
f_prompt (m);
(void) read_char();
}
/* crack a line of the form X?X?X into its components,
* where X is an integer and ? can be any character except '0-9' or '-',
* such as ':' or '/'.
* only change those fields that are specified:
* eg: ::10 only changes *s
* 10 only changes *d
* 10:0 changes *d and *m
* if see '-' anywhere, first non-zero component will be made negative.
*/
f_sscansex (bp, d, m, s)
char *bp;
int *d, *m, *s;
{
char c;
int *p = d;
int *nonzp = 0;
int sawneg = 0;
int innum = 0;
while (c = *bp++)
if (c >= '0' && c <= '9') {
if (!innum) {
*p = 0;
innum = 1;
}
*p = *p*10 + (c - '0');
if (*p && !nonzp)
nonzp = p;
} else if (c == '-') {
sawneg = 1;
} else if (c != ' ') {
/* advance to next component */
p = (p == d) ? m : s;
innum = 0;
}
if (sawneg && nonzp)
*nonzp = -*nonzp;
}
/* just like dec_sex() but makes the first non-zero element negative if
* x is negative (instead of returning a sign flag).
*/
f_dec_sexsign (x, h, m, s)
double x;
int *h, *m, *s;
{
int n;
dec_sex (x, h, m, s, &n);
if (n) {
if (*h)
*h = -*h;
else if (*m)
*m = -*m;
else
*s = -*s;
}
}
xXx