downey@cs.umn.edu@dimed1.UUCP (08/31/90)
Posting-number: Volume 14, Issue 79
Submitted-by: downey@cs.umn.edu@dimed1.UUCP
Archive-name: ephem-4.21/part04
#! /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 4 (of 6)."
# Contents: compiler.c plans.c watch.c
# Wrapped by allbery@uunet on Thu Aug 30 20:46:35 1990
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'compiler.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'compiler.c'\"
else
echo shar: Extracting \"'compiler.c'\" \(15331 characters\)
sed "s/^X//" >'compiler.c' <<'END_OF_FILE'
X/* module to compile and execute a c-style arithmetic expression.
X * public entry points are compile_expr() and execute_expr().
X *
X * one reason this is so nice and tight is that all opcodes are the same size
X * (an int) and the tokens the parser returns are directly usable as opcodes,
X * for the most part. constants and variables are compiled as an opcode
X * with an offset into the auxiliary opcode tape, opx.
X */
X
X#include <math.h>
X#ifdef VMS
X#include <stdlib.h>
X#endif
X#include "screen.h"
X
X/* parser tokens and opcodes, as necessary */
X#define HALT 0 /* good value for HALT since program is inited to 0 */
X/* binary operators (precedences in table, below) */
X#define ADD 1
X#define SUB 2
X#define MULT 3
X#define DIV 4
X#define AND 5
X#define OR 6
X#define GT 7
X#define GE 8
X#define EQ 9
X#define NE 10
X#define LT 11
X#define LE 12
X/* unary op, precedence in NEG_PREC #define, below */
X#define NEG 13
X/* symantically operands, ie, constants, variables and all functions */
X#define CONST 14
X#define VAR 15
X#define ABS 16 /* add functions if desired just like this is done */
X/* purely tokens - never get compiled as such */
X#define LPAREN 255
X#define RPAREN 254
X#define ERR (-1)
X
X/* precedence of each of the binary operators.
X * in case of a tie, compiler associates left-to-right.
X * N.B. each entry's index must correspond to its #define!
X */
Xstatic int precedence[] = {0,5,5,6,6,2,1,4,4,3,3,4,4};
X#define NEG_PREC 7 /* negation is highest */
X
X/* execute-time operand stack */
X#define MAX_STACK 16
Xstatic double stack[MAX_STACK], *sp;
X
X/* space for compiled opcodes - the "program".
X * opcodes go in lower 8 bits.
X * when an opcode has an operand (as CONST and VAR) it is really in opx[] and
X * the index is in the remaining upper bits.
X */
X#define MAX_PROG 32
Xstatic int program[MAX_PROG], *pc;
X#define OP_SHIFT 8
X#define OP_MASK 0xff
X
X/* auxiliary operand info.
X * the operands (all but lower 8 bits) of CONST and VAR are really indeces
X * into this array. thus, no point in making this any longer than you have
X * bits more than 8 in your machine's int to index into it, ie, make
X * MAX_OPX <= 1 << ((sizeof(int)-1)*8)
X * also, the fld's must refer to ones being flog'd, so not point in more
X * of these then that might be used for plotting and srching combined.
X */
X#define MAX_OPX 16
Xtypedef union {
X double opu_f; /* value when opcode is CONST */
X int opu_fld; /* rcfpack() of field when opcode is VAR */
X} OpX;
Xstatic OpX opx[MAX_OPX];
Xstatic int opxidx;
X
X/* these are global just for easy/rapid access */
Xstatic int parens_nest; /* to check that parens end up nested */
Xstatic char *err_msg; /* caller provides storage; we point at it with this */
Xstatic char *cexpr, *lcexpr; /* pointers that move along caller's expression */
Xstatic int good_prog; /* != 0 when program appears to be good */
X
X/* compile the given c-style expression.
X * return 0 and set good_prog if ok,
X * else return -1 and a reason message in errbuf.
X */
Xcompile_expr (ex, errbuf)
Xchar *ex;
Xchar *errbuf;
X{
X int instr;
X
X /* init the globals.
X * also delete any flogs used in the previous program.
X */
X cexpr = ex;
X err_msg = errbuf;
X pc = program;
X opxidx = 0;
X parens_nest = 0;
X do {
X instr = *pc++;
X if ((instr & OP_MASK) == VAR)
X flog_delete (opx[instr >> OP_SHIFT].opu_fld);
X } while (instr != HALT);
X
X pc = program;
X if (compile(0) == ERR) {
X (void) sprintf (err_msg + strlen(err_msg), " at \"%.10s\"", lcexpr);
X good_prog = 0;
X return (-1);
X }
X *pc++ = HALT;
X good_prog = 1;
X return (0);
X}
X
X/* execute the expression previously compiled with compile_expr().
X * return 0 with *vp set to the answer if ok, else return -1 with a reason
X * why not message in errbuf.
X */
Xexecute_expr (vp, errbuf)
Xdouble *vp;
Xchar *errbuf;
X{
X int s;
X
X err_msg = errbuf;
X sp = stack + MAX_STACK; /* grows towards lower addresses */
X pc = program;
X s = execute(vp);
X if (s < 0)
X good_prog = 0;
X return (s);
X}
X
X/* this is a way for the outside world to ask whether there is currently a
X * reasonable program compiled and able to execute.
X */
Xprog_isgood()
X{
X return (good_prog);
X}
X
X/* get and return the opcode corresponding to the next token.
X * leave with lcexpr pointing at the new token, cexpr just after it.
X * also watch for mismatches parens and proper operator/operand alternation.
X */
Xstatic
Xnext_token ()
X{
X static char toomt[] = "More than %d terms";
X static char badop[] = "Illegal operator";
X int tok = ERR; /* just something illegal */
X char c;
X
X while ((c = *cexpr) == ' ')
X cexpr++;
X lcexpr = cexpr++;
X
X /* mainly check for a binary operator */
X switch (c) {
X case '\0': --cexpr; tok = HALT; break; /* keep returning HALT */
X case '+': tok = ADD; break; /* compiler knows when it's really unary */
X case '-': tok = SUB; break; /* compiler knows when it's really negate */
X case '*': tok = MULT; break;
X case '/': tok = DIV; break;
X case '(': parens_nest++; tok = LPAREN; break;
X case ')':
X if (--parens_nest < 0) {
X (void) sprintf (err_msg, "Too many right parens");
X return (ERR);
X } else
X tok = RPAREN;
X break;
X case '|':
X if (*cexpr == '|') { cexpr++; tok = OR; }
X else { (void) sprintf (err_msg, badop); return (ERR); }
X break;
X case '&':
X if (*cexpr == '&') { cexpr++; tok = AND; }
X else { (void) sprintf (err_msg, badop); return (ERR); }
X break;
X case '=':
X if (*cexpr == '=') { cexpr++; tok = EQ; }
X else { (void) sprintf (err_msg, badop); return (ERR); }
X break;
X case '!':
X if (*cexpr == '=') { cexpr++; tok = NE; }
X else { (void) sprintf (err_msg, badop); return (ERR); }
X break;
X case '<':
X if (*cexpr == '=') { cexpr++; tok = LE; }
X else tok = LT;
X break;
X case '>':
X if (*cexpr == '=') { cexpr++; tok = GE; }
X else tok = GT;
X break;
X }
X
X if (tok != ERR)
X return (tok);
X
X /* not op so check for a constant, variable or function */
X if (isdigit(c) || c == '.') {
X if (opxidx > MAX_OPX) {
X (void) sprintf (err_msg, toomt, MAX_OPX);
X return (ERR);
X }
X opx[opxidx].opu_f = atof (lcexpr);
X tok = CONST | (opxidx++ << OP_SHIFT);
X skip_double();
X } else if (isalpha(c)) {
X /* check list of functions */
X if (strncmp (lcexpr, "abs", 3) == 0) {
X cexpr += 2;
X tok = ABS;
X } else {
X /* not a function, so assume it's a variable */
X int fld;
X if (opxidx > MAX_OPX) {
X (void) sprintf (err_msg, toomt, MAX_OPX);
X return (ERR);
X }
X fld = parse_fieldname ();
X if (fld < 0) {
X (void) sprintf (err_msg, "Unknown field");
X return (ERR);
X } else {
X if (flog_add (fld) < 0) { /* register with field logger */
X (void) sprintf (err_msg, "Sorry; too many fields");
X return (ERR);
X }
X opx[opxidx].opu_fld = fld;
X tok = VAR | (opxidx++ << OP_SHIFT);
X }
X }
X }
X
X return (tok);
X}
X
X/* move cexpr on past a double.
X * allow sci notation.
X * no need to worry about a leading '-' or '+' but allow them after an 'e'.
X * TODO: this handles all the desired cases, but also admits a bit too much
X * such as things like 1eee2...3. geeze; to skip a double right you almost
X * have to go ahead and crack it!
X */
Xstatic
Xskip_double()
X{
X int sawe = 0; /* so we can allow '-' or '+' right after an 'e' */
X
X while (1) {
X char c = *cexpr;
X if (isdigit(c) || c=='.' || (sawe && (c=='-' || c=='+'))) {
X sawe = 0;
X cexpr++;
X } else if (c == 'e') {
X sawe = 1;
X cexpr++;
X } else
X break;
X }
X}
X
X/* call this whenever you want to dig out the next (sub)expression.
X * keep compiling instructions as long as the operators are higher precedence
X * than prec, then return that "look-ahead" token that wasn't (higher prec).
X * if error, fill in a message in err_msg[] and return ERR.
X */
Xstatic
Xcompile (prec)
Xint prec;
X{
X int expect_binop = 0; /* set after we have seen any operand.
X * used by SUB so it can tell if it really
X * should be taken to be a NEG instead.
X */
X int tok = next_token ();
X
X while (1) {
X int p;
X if (tok == ERR)
X return (ERR);
X if (pc - program >= MAX_PROG) {
X (void) sprintf (err_msg, "Program is too long");
X return (ERR);
X }
X
X /* check for special things like functions, constants and parens */
X switch (tok & OP_MASK) {
X case HALT: return (tok);
X case ADD:
X if (expect_binop)
X break; /* procede with binary addition */
X /* just skip a unary positive(?) */
X tok = next_token();
X continue;
X case SUB:
X if (expect_binop)
X break; /* procede with binary subtract */
X tok = compile (NEG_PREC);
X *pc++ = NEG;
X expect_binop = 1;
X continue;
X case ABS: /* other funcs would be handled the same too ... */
X /* eat up the function parenthesized argument */
X if (next_token() != LPAREN || compile (0) != RPAREN) {
X (void) sprintf (err_msg, "Function arglist error");
X return (ERR);
X }
X /* then handled same as ... */
X case CONST: /* handled same as... */
X case VAR:
X *pc++ = tok;
X tok = next_token();
X expect_binop = 1;
X continue;
X case LPAREN:
X if (compile (0) != RPAREN) {
X (void) sprintf (err_msg, "Unmatched left paren");
X return (ERR);
X }
X tok = next_token();
X expect_binop = 1;
X continue;
X case RPAREN:
X return (RPAREN);
X }
X
X /* everything else is a binary operator */
X p = precedence[tok];
X if (p > prec) {
X int newtok = compile (p);
X if (newtok == ERR)
X return (ERR);
X *pc++ = tok;
X expect_binop = 1;
X tok = newtok;
X } else
X return (tok);
X }
X}
X
X/* "run" the program[] compiled with compile().
X * if ok, return 0 and the final result,
X * else return -1 with a reason why not message in err_msg.
X */
Xstatic
Xexecute(result)
Xdouble *result;
X{
X int instr;
X
X do {
X instr = *pc++;
X switch (instr & OP_MASK) {
X /* put these in numberic order so hopefully even the dumbest
X * compiler will choose to use a jump table, not a cascade of ifs.
X */
X case HALT: break; /* outer loop will stop us */
X case ADD: sp[1] = sp[1] + sp[0]; sp++; break;
X case SUB: sp[1] = sp[1] - sp[0]; sp++; break;
X case MULT: sp[1] = sp[1] * sp[0]; sp++; break;
X case DIV: sp[1] = sp[1] / sp[0]; sp++; break;
X case AND: sp[1] = sp[1] && sp[0] ? 1 : 0; sp++; break;
X case OR: sp[1] = sp[1] || sp[0] ? 1 : 0; sp++; break;
X case GT: sp[1] = sp[1] > sp[0] ? 1 : 0; sp++; break;
X case GE: sp[1] = sp[1] >= sp[0] ? 1 : 0; sp++; break;
X case EQ: sp[1] = sp[1] == sp[0] ? 1 : 0; sp++; break;
X case NE: sp[1] = sp[1] != sp[0] ? 1 : 0; sp++; break;
X case LT: sp[1] = sp[1] < sp[0] ? 1 : 0; sp++; break;
X case LE: sp[1] = sp[1] <= sp[0] ? 1 : 0; sp++; break;
X case NEG: *sp = -*sp; break;
X case CONST: *--sp = opx[instr >> OP_SHIFT].opu_f; break;
X case VAR:
X if (flog_get(opx[instr>>OP_SHIFT].opu_fld, --sp, (char *)0)<0) {
X (void) sprintf (err_msg, "Bug! VAR field not logged");
X return (-1);
X }
X break;
X case ABS: *sp = fabs (*sp); break;
X default:
X (void) sprintf (err_msg, "Bug! bad opcode: 0x%x", instr);
X return (-1);
X }
X if (sp < stack) {
X (void) sprintf (err_msg, "Runtime stack overflow");
X return (-1);
X } else if (sp - stack > MAX_STACK) {
X (void) sprintf (err_msg, "Bug! runtime stack underflow");
X return (-1);
X }
X } while (instr != HALT);
X
X /* result should now be on top of stack */
X if (sp != &stack[MAX_STACK - 1]) {
X (void) sprintf (err_msg, "Bug! stack has %d items",
X MAX_STACK - (sp-stack));
X return (-1);
X }
X *result = *sp;
X return (0);
X}
X
Xstatic
Xisdigit(c)
Xchar c;
X{
X return (c >= '0' && c <= '9');
X}
X
Xstatic
Xisalpha (c)
Xchar c;
X{
X return ((c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z'));
X}
X
X/* starting with lcexpr pointing at a string expected to be a field name,
X * return an rcfpack(r,c,0) of the field else -1 if bad.
X * when return, leave lcexpr alone but move cexpr to just after the name.
X */
Xstatic
Xparse_fieldname ()
X{
X int r = -1, c = -1; /* anything illegal */
X char *fn = lcexpr; /* likely faster than using the global */
X char f0, f1;
X char *dp;
X
X /* search for first thing not an alpha char.
X * leave it in f0 and leave dp pointing to it.
X */
X dp = fn;
X while (isalpha(f0 = *dp))
X dp++;
X
X /* crack the new field name.
X * when done trying, leave dp pointing at first char just after it.
X * set r and c if we recognized it.
X */
X if (f0 == '.') {
X /* planet.column pair.
X * first crack the planet portion (pointed to by fn): set r.
X * then the second portion (pointed to by dp+1): set c.
X */
X f0 = fn[0];
X f1 = fn[1];
X switch (f0) {
X case 'j':
X r = R_JUPITER;
X break;
X case 'm':
X if (f1 == 'a') r = R_MARS;
X else if (f1 == 'e') r = R_MERCURY;
X else if (f1 == 'o') r = R_MOON;
X break;
X case 'n':
X r = R_NEPTUNE;
X break;
X case 'p':
X r = R_PLUTO;
X break;
X case 's':
X if (f1 == 'a') r = R_SATURN;
X else if (f1 == 'u') r = R_SUN;
X break;
X case 'u':
X r = R_URANUS;
X break;
X case 'x':
X r = R_OBJX;
X break;
X case 'y':
X r = R_OBJY;
X break;
X case 'v':
X r = R_VENUS;
X break;
X }
X
X /* now crack the column (stuff after the dp) */
X dp++; /* point at good stuff just after the decimal pt */
X f0 = dp[0];
X f1 = dp[1];
X switch (f0) {
X case 'a':
X if (f1 == 'l') c = C_ALT;
X else if (f1 == 'z') c = C_AZ;
X break;
X case 'd':
X c = C_DEC;
X break;
X case 'e':
X if (f1 == 'd') c = C_EDIST;
X else if (f1 == 'l') c = C_ELONG;
X break;
X case 'h':
X if (f1 == 'l') {
X if (dp[2] == 'a') c = C_HLAT;
X else if (dp[2] == 'o') c = C_HLONG;
X } else if (f1 == 'r' || f1 == 'u') c = C_TUP;
X break;
X case 'j':
X c = C_JUPITER;
X break;
X case 'm':
X if (f1 == 'a') c = C_MARS;
X else if (f1 == 'e') c = C_MERCURY;
X else if (f1 == 'o') c = C_MOON;
X break;
X case 'n':
X c = C_NEPTUNE;
X break;
X case 'p':
X if (f1 == 'h') c = C_PHASE;
X else if (f1 == 'l') c = C_PLUTO;
X break;
X case 'r':
X if (f1 == 'a') {
X if (dp[2] == 'z') c = C_RISEAZ;
X else c = C_RA;
X } else if (f1 == 't') c = C_RISETM;
X break;
X case 's':
X if (f1 == 'a') {
X if (dp[2] == 'z') c = C_SETAZ;
X else c = C_SATURN;
X } else if (f1 == 'd') c = C_SDIST;
X else if (f1 == 'i') c = C_SIZE;
X else if (f1 == 't') c = C_SETTM;
X else if (f1 == 'u') c = C_SUN;
X break;
X case 't':
X if (f1 == 'a') c = C_TRANSALT;
X else if (f1 == 't') c = C_TRANSTM;
X break;
X case 'u':
X c = C_URANUS;
X break;
X case 'x':
X c = C_OBJX;
X break;
X case 'y':
X c = C_OBJY;
X break;
X case 'v':
X if (f1 == 'e') c = C_VENUS;
X else if (f1 == 'm') c = C_MAG;
X break;
X }
X
X /* now skip dp on past the column stuff */
X while (isalpha(*dp))
X dp++;
X } else {
X /* no decimal point; some field in the top of the screen */
X f0 = fn[0];
X f1 = fn[1];
X switch (f0) {
X case 'd':
X if (f1 == 'a') r = R_DAWN, c = C_DAWNV;
X else if (f1 == 'u') r = R_DUSK, c = C_DUSKV;
X break;
X case 'n':
X r = R_LON, c = C_LONV;
X break;
X }
X }
X
X cexpr = dp;
X if (r <= 0 || c <= 0) return (-1);
X return (rcfpack (r, c, 0));
X}
END_OF_FILE
if test 15331 -ne `wc -c <'compiler.c'`; then
echo shar: \"'compiler.c'\" unpacked with wrong size!
fi
# end of 'compiler.c'
fi
if test -f 'plans.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'plans.c'\"
else
echo shar: Extracting \"'plans.c'\" \(17647 characters\)
sed "s/^X//" >'plans.c' <<'END_OF_FILE'
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X
X#define TWOPI (2*PI)
X#define mod2PI(x) ((x) - (long)((x)/TWOPI)*TWOPI)
X
X/* given a modified Julian date, mjd, and a planet, p, find:
X * lpd0: heliocentric longitude,
X * psi0: heliocentric latitude,
X * rp0: distance from the sun to the planet,
X * rho0: distance from the Earth to the planet,
X * none corrected for light time, ie, they are the true values for the
X * given instant.
X * lam: geocentric ecliptic longitude,
X * bet: geocentric ecliptic latitude,
X * each corrected for light time, ie, they are the apparent values as
X * seen from the center of the Earth for the given instant.
X * dia: angular diameter in arcsec at 1 AU,
X * mag: visual magnitude when 1 AU from sun and earth at 0 phase angle.
X *
X * all angles are in radians, all distances in AU.
X * the mean orbital elements are found by calling pelement(), then mutual
X * perturbation corrections are applied as necessary.
X *
X * corrections for nutation and abberation must be made by the caller. The RA
X * and DEC calculated from the fully-corrected ecliptic coordinates are then
X * the apparent geocentric coordinates. Further corrections can be made, if
X * required, for atmospheric refraction and geocentric parallax although the
X * intrinsic error herein of about 10 arcseconds is usually the dominant
X * error at this stage.
X * TODO: combine the several intermediate expressions when get a good compiler.
X */
Xplans (mjd, p, lpd0, psi0, rp0, rho0, lam, bet, dia, mag)
Xdouble mjd;
Xint p;
Xdouble *lpd0, *psi0, *rp0, *rho0, *lam, *bet, *dia, *mag;
X{
X static double plan[8][9];
X static double lastmjd = -10000;
X double dl; /* perturbation correction for longitude */
X double dr; /* " orbital radius */
X double dml; /* " mean longitude */
X double ds; /* " eccentricity */
X double dm; /* " mean anomaly */
X double da; /* " semi-major axis */
X double dhl; /* " heliocentric longitude */
X double lsn, rsn;/* true geocentric longitude of sun and sun-earth rad */
X double mas; /* mean anomaly of the sun */
X double re; /* radius of earth's orbit */
X double lg; /* longitude of earth */
X double map[8]; /* array of mean anomalies for each planet */
X double lpd, psi, rp, rho;
X double ll, sll, cll;
X double t;
X double dt;
X int pass;
X int j;
X double s, ma;
X double nu, ea;
X double lp, om;
X double lo, slo, clo;
X double inc, y;
X double spsi, cpsi;
X double rpd;
X
X /* only need to fill in plan[] once for a given mjd */
X if (mjd != lastmjd) {
X pelement (mjd, plan);
X lastmjd = mjd;
X }
X
X dt = 0;
X t = mjd/36525.;
X sunpos (mjd, &lsn, &rsn);
X masun (mjd, &mas);
X re = rsn;
X lg = lsn+PI;
X
X /* first find the true position of the planet at mjd.
X * then repeat a second time for a slightly different time based
X * on the position found in the first pass to account for light-travel
X * time.
X */
X for (pass = 0; pass < 2; pass++) {
X
X for (j = 0; j < 8; j++)
X map[j] = degrad(plan[j][0]-plan[j][2]-dt*plan[j][1]);
X
X /* set initial corrections to 0.
X * then modify as necessary for the planet of interest.
X */
X dl = 0;
X dr = 0;
X dml = 0;
X ds = 0;
X dm = 0;
X da = 0;
X dhl = 0;
X
X switch (p) {
X
X case MERCURY:
X p_mercury (map, &dl, &dr);
X break;
X
X case VENUS:
X p_venus (t, mas, map, &dl, &dr, &dml, &dm);
X break;
X
X case MARS:
X p_mars (mas, map, &dl, &dr, &dml, &dm);
X break;
X
X case JUPITER:
X p_jupiter (t, plan[p][3], &dml, &ds, &dm, &da);
X break;
X
X case SATURN:
X p_saturn (t, plan[p][3], &dml, &ds, &dm, &da, &dhl);
X break;
X
X case URANUS:
X p_uranus (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
X break;
X
X case NEPTUNE:
X p_neptune (t, plan[p][3], &dl, &dr, &dml, &ds, &dm, &da, &dhl);
X break;
X
X case PLUTO:
X /* no perturbation theory for pluto */
X break;
X }
X
X s = plan[p][3]+ds;
X ma = map[p]+dm;
X anomaly (ma, s, &nu, &ea);
X rp = (plan[p][6]+da)*(1-s*s)/(1+s*cos(nu));
X lp = raddeg(nu)+plan[p][2]+raddeg(dml-dm);
X lp = degrad(lp);
X om = degrad(plan[p][5]);
X lo = lp-om;
X slo = sin(lo);
X clo = cos(lo);
X inc = degrad(plan[p][4]);
X rp = rp+dr;
X spsi = slo*sin(inc);
X y = slo*cos(inc);
X psi = asin(spsi)+dhl;
X spsi = sin(psi);
X lpd = atan(y/clo)+om+degrad(dl);
X if (clo<0) lpd += PI;
X range (&lpd, TWOPI);
X cpsi = cos(psi);
X rpd = rp*cpsi;
X ll = lpd-lg;
X rho = sqrt(re*re+rp*rp-2*re*rp*cpsi*cos(ll));
X
X /* when we view a planet we see it in the position it occupied
X * dt days ago, where rho is the distance between it and earth,
X * in AU. use this as the new time for the next pass.
X */
X dt = rho*5.775518e-3;
X
X if (pass == 0) {
X /* save heliocentric coordinates after first pass since, being
X * true, they are NOT to be corrected for light-travel time.
X */
X *lpd0 = lpd;
X range (lpd0, TWOPI);
X *psi0 = psi;
X *rp0 = rp;
X *rho0 = rho;
X }
X }
X
X sll = sin(ll);
X cll = cos(ll);
X if (p < MARS)
X *lam = atan(-1*rpd*sll/(re-rpd*cll))+lg+PI;
X else
X *lam = atan(re*sll/(rpd-re*cll))+lpd;
X range (lam, TWOPI);
X *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*re*sll));
X *dia = plan[p][7];
X *mag = plan[p][8];
X}
X
X/* set auxilliary variables used for jupiter, saturn, uranus, and neptune */
Xstatic
Xaux_jsun (t, x1, x2, x3, x4, x5, x6)
Xdouble t;
Xdouble *x1, *x2, *x3, *x4, *x5, *x6;
X{
X *x1 = t/5+0.1;
X *x2 = mod2PI(4.14473+5.29691e1*t);
X *x3 = mod2PI(4.641118+2.132991e1*t);
X *x4 = mod2PI(4.250177+7.478172*t);
X *x5 = 5 * *x3 - 2 * *x2;
X *x6 = 2 * *x2 - 6 * *x3 + 3 * *x4;
X}
X
X/* find the mean anomaly of the sun at mjd.
X * this is the same as that used in sun() but when it was converted to C it
X * was not known it would be required outside that routine.
X * TODO: add an argument to sun() to return mas and eliminate this routine.
X */
Xstatic
Xmasun (mjd, mas)
Xdouble mjd;
Xdouble *mas;
X{
X double t, t2;
X double a, b;
X
X t = mjd/36525;
X t2 = t*t;
X a = 9.999736042e1*t;
X b = 360.*(a-(long)a);
X *mas = degrad (3.5847583e2-(1.5e-4+3.3e-6*t)*t2+b);
X}
X
X/* perturbations for mercury */
Xstatic
Xp_mercury (map, dl, dr)
Xdouble map[];
Xdouble *dl, *dr;
X{
X *dl = 2.04e-3*cos(5*map[2-1]-2*map[1-1]+2.1328e-1)+
X 1.03e-3*cos(2*map[2-1]-map[1-1]-2.8046)+
X 9.1e-4*cos(2*map[3]-map[1-1]-6.4582e-1)+
X 7.8e-4*cos(5*map[2-1]-3*map[1-1]+1.7692e-1);
X
X *dr = 7.525e-6*cos(2*map[3]-map[1-1]+9.25251e-1)+
X 6.802e-6*cos(5*map[2-1]-3*map[1-1]-4.53642)+
X 5.457e-6*cos(2*map[2-1]-2*map[1-1]-1.24246)+
X 3.569e-6*cos(5*map[2-1]-map[1-1]-1.35699);
X}
X
X/* ....venus */
Xstatic
Xp_venus (t, mas, map, dl, dr, dml, dm)
Xdouble t, mas, map[];
Xdouble *dl, *dr, *dml, *dm;
X{
X *dml = degrad (7.7e-4*sin(4.1406+t*2.6227));
X *dm = *dml;
X
X *dl = 3.13e-3*cos(2*mas-2*map[2-1]-2.587)+
X 1.98e-3*cos(3*mas-3*map[2-1]+4.4768e-2)+
X 1.36e-3*cos(mas-map[2-1]-2.0788)+
X 9.6e-4*cos(3*mas-2*map[2-1]-2.3721)+
X 8.2e-4*cos(map[3]-map[2-1]-3.6318);
X
X *dr = 2.2501e-5*cos(2*mas-2*map[2-1]-1.01592)+
X 1.9045e-5*cos(3*mas-3*map[2-1]+1.61577)+
X 6.887e-6*cos(map[3]-map[2-1]-2.06106)+
X 5.172e-6*cos(mas-map[2-1]-5.08065e-1)+
X 3.62e-6*cos(5*mas-4*map[2-1]-1.81877)+
X 3.283e-6*cos(4*mas-4*map[2-1]+1.10851)+
X 3.074e-6*cos(2*map[3]-2*map[2-1]-9.62846e-1);
X}
X
X/* ....mars */
Xstatic
Xp_mars (mas, map, dl, dr, dml, dm)
Xdouble mas, map[];
Xdouble *dl, *dr, *dml, *dm;
X{
X double a;
X
X a = 3*map[3]-8*map[2]+4*mas;
X *dml = degrad (-1*(1.133e-2*sin(a)+9.33e-3*cos(a)));
X *dm = *dml;
X
X *dl = 7.05e-3*cos(map[3]-map[2]-8.5448e-1)+
X 6.07e-3*cos(2*map[3]-map[2]-3.2873)+
X 4.45e-3*cos(2*map[3]-2*map[2]-3.3492)+
X 3.88e-3*cos(mas-2*map[2]+3.5771e-1)+
X 2.38e-3*cos(mas-map[2]+6.1256e-1)+
X 2.04e-3*cos(2*mas-3*map[2]+2.7688)+
X 1.77e-3*cos(3*map[2]-map[2-1]-1.0053)+
X 1.36e-3*cos(2*mas-4*map[2]+2.6894)+
X 1.04e-3*cos(map[3]+3.0749e-1);
X
X *dr = 5.3227e-5*cos(map[3]-map[2]+7.17864e-1)+
X 5.0989e-5*cos(2*map[3]-2*map[2]-1.77997)+
X 3.8278e-5*cos(2*map[3]-map[2]-1.71617)+
X 1.5996e-5*cos(mas-map[2]-9.69618e-1)+
X 1.4764e-5*cos(2*mas-3*map[2]+1.19768)+
X 8.966e-6*cos(map[3]-2*map[2]+7.61225e-1);
X *dr += 7.914e-6*cos(3*map[3]-2*map[2]-2.43887)+
X 7.004e-6*cos(2*map[3]-3*map[2]-1.79573)+
X 6.62e-6*cos(mas-2*map[2]+1.97575)+
X 4.93e-6*cos(3*map[3]-3*map[2]-1.33069)+
X 4.693e-6*cos(3*mas-5*map[2]+3.32665)+
X 4.571e-6*cos(2*mas-4*map[2]+4.27086)+
X 4.409e-6*cos(3*map[3]-map[2]-2.02158);
X}
X
X/* ....jupiter */
Xstatic
Xp_jupiter (t, s, dml, ds, dm, da)
Xdouble t, s;
Xdouble *dml, *ds, *dm, *da;
X{
X double dp;
X double x1, x2, x3, x4, x5, x6, x7;
X double sx3, cx3, s2x3, c2x3;
X double sx5, cx5, s2x5;
X double sx6;
X double sx7, cx7, s2x7, c2x7, s3x7, c3x7, s4x7, c4x7, c5x7;
X
X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
X x7 = x3-x2;
X sx3 = sin(x3);
X cx3 = cos(x3);
X s2x3 = sin(2*x3);
X c2x3 = cos(2*x3);
X sx5 = sin(x5);
X cx5 = cos(x5);
X s2x5 = sin(2*x5);
X sx6 = sin(x6);
X sx7 = sin(x7);
X cx7 = cos(x7);
X s2x7 = sin(2*x7);
X c2x7 = cos(2*x7);
X s3x7 = sin(3*x7);
X c3x7 = cos(3*x7);
X s4x7 = sin(4*x7);
X c4x7 = cos(4*x7);
X c5x7 = cos(5*x7);
X
X *dml = (3.31364e-1-(1.0281e-2+4.692e-3*x1)*x1)*sx5+
X (3.228e-3-(6.4436e-2-2.075e-3*x1)*x1)*cx5-
X (3.083e-3+(2.75e-4-4.89e-4*x1)*x1)*s2x5+
X 2.472e-3*sx6+1.3619e-2*sx7+1.8472e-2*s2x7+6.717e-3*s3x7+
X 2.775e-3*s4x7+6.417e-3*s2x7*sx3+
X (7.275e-3-1.253e-3*x1)*sx7*sx3+
X 2.439e-3*s3x7*sx3-(3.5681e-2+1.208e-3*x1)*sx7*cx3;
X *dml += -3.767e-3*c2x7*sx3-(3.3839e-2+1.125e-3*x1)*cx7*sx3-
X 4.261e-3*s2x7*cx3+
X (1.161e-3*x1-6.333e-3)*cx7*cx3+
X 2.178e-3*cx3-6.675e-3*c2x7*cx3-2.664e-3*c3x7*cx3-
X 2.572e-3*sx7*s2x3-3.567e-3*s2x7*s2x3+2.094e-3*cx7*c2x3+
X 3.342e-3*c2x7*c2x3;
X *dml = degrad(*dml);
X
X *ds = (3606+(130-43*x1)*x1)*sx5+(1289-580*x1)*cx5-6764*sx7*sx3-
X 1110*s2x7*sx3-224*s3x7*sx3-204*sx3+(1284+116*x1)*cx7*sx3+
X 188*c2x7*sx3+(1460+130*x1)*sx7*cx3+224*s2x7*cx3-817*cx3+
X 6074*cx3*cx7+992*c2x7*cx3+
X 508*c3x7*cx3+230*c4x7*cx3+108*c5x7*cx3;
X *ds += -(956+73*x1)*sx7*s2x3+448*s2x7*s2x3+137*s3x7*s2x3+
X (108*x1-997)*cx7*s2x3+480*c2x7*s2x3+148*c3x7*s2x3+
X (99*x1-956)*sx7*c2x3+490*s2x7*c2x3+
X 158*s3x7*c2x3+179*c2x3+(1024+75*x1)*cx7*c2x3-
X 437*c2x7*c2x3-132*c3x7*c2x3;
X *ds *= 1e-7;
X
X dp = (7.192e-3-3.147e-3*x1)*sx5-4.344e-3*sx3+
X (x1*(1.97e-4*x1-6.75e-4)-2.0428e-2)*cx5+
X 3.4036e-2*cx7*sx3+(7.269e-3+6.72e-4*x1)*sx7*sx3+
X 5.614e-3*c2x7*sx3+2.964e-3*c3x7*sx3+3.7761e-2*sx7*cx3+
X 6.158e-3*s2x7*cx3-
X 6.603e-3*cx7*cx3-5.356e-3*sx7*s2x3+2.722e-3*s2x7*s2x3+
X 4.483e-3*cx7*s2x3-2.642e-3*c2x7*s2x3+4.403e-3*sx7*c2x3-
X 2.536e-3*s2x7*c2x3+5.547e-3*cx7*c2x3-2.689e-3*c2x7*c2x3;
X
X *dm = *dml-(degrad(dp)/s);
X
X *da = 205*cx7-263*cx5+693*c2x7+312*c3x7+147*c4x7+299*sx7*sx3+
X 181*c2x7*sx3+204*s2x7*cx3+111*s3x7*cx3-337*cx7*cx3-
X 111*c2x7*cx3;
X *da *= 1e-6;
X}
X
X/* ....saturn */
Xstatic
Xp_saturn (t, s, dml, ds, dm, da, dhl)
Xdouble t, s;
Xdouble *dml, *ds, *dm, *da, *dhl;
X{
X double dp;
X double x1, x2, x3, x4, x5, x6, x7, x8;
X double sx3, cx3, s2x3, c2x3, s3x3, c3x3, s4x3, c4x3;
X double sx5, cx5, s2x5, c2x5;
X double sx6;
X double sx7, cx7, s2x7, c2x7, s3x7, c3x7, s4x7, c4x7, c5x7, s5x7;
X double s2x8, c2x8, s3x8, c3x8;
X
X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
X x7 = x3-x2;
X sx3 = sin(x3);
X cx3 = cos(x3);
X s2x3 = sin(2*x3);
X c2x3 = cos(2*x3);
X sx5 = sin(x5);
X cx5 = cos(x5);
X s2x5 = sin(2*x5);
X sx6 = sin(x6);
X sx7 = sin(x7);
X cx7 = cos(x7);
X s2x7 = sin(2*x7);
X c2x7 = cos(2*x7);
X s3x7 = sin(3*x7);
X c3x7 = cos(3*x7);
X s4x7 = sin(4*x7);
X c4x7 = cos(4*x7);
X c5x7 = cos(5*x7);
X
X s3x3 = sin(3*x3);
X c3x3 = cos(3*x3);
X s4x3 = sin(4*x3);
X c4x3 = cos(4*x3);
X c2x5 = cos(2*x5);
X s5x7 = sin(5*x7);
X x8 = x4-x3;
X s2x8 = sin(2*x8);
X c2x8 = cos(2*x8);
X s3x8 = sin(3*x8);
X c3x8 = cos(3*x8);
X
X *dml = 7.581e-3*s2x5-7.986e-3*sx6-1.48811e-1*sx7-4.0786e-2*s2x7-
X (8.14181e-1-(1.815e-2-1.6714e-2*x1)*x1)*sx5-
X (1.0497e-2-(1.60906e-1-4.1e-3*x1)*x1)*cx5-1.5208e-2*s3x7-
X 6.339e-3*s4x7-6.244e-3*sx3-1.65e-2*s2x7*sx3+
X (8.931e-3+2.728e-3*x1)*sx7*sx3-5.775e-3*s3x7*sx3+
X (8.1344e-2+3.206e-3*x1)*cx7*sx3+1.5019e-2*c2x7*sx3;
X *dml += (8.5581e-2+2.494e-3*x1)*sx7*cx3+1.4394e-2*c2x7*cx3+
X (2.5328e-2-3.117e-3*x1)*cx7*cx3+
X 6.319e-3*c3x7*cx3+6.369e-3*sx7*s2x3+9.156e-3*s2x7*s2x3+
X 7.525e-3*s3x8*s2x3-5.236e-3*cx7*c2x3-7.736e-3*c2x7*c2x3-
X 7.528e-3*c3x8*c2x3;
X *dml = degrad(*dml);
X
X *ds = (-7927+(2548+91*x1)*x1)*sx5+(13381+(1226-253*x1)*x1)*cx5+
X (248-121*x1)*s2x5-(305+91*x1)*c2x5+412*s2x7+12415*sx3+
X (390-617*x1)*sx7*sx3+(165-204*x1)*s2x7*sx3+26599*cx7*sx3-
X 4687*c2x7*sx3-1870*c3x7*sx3-821*c4x7*sx3-
X 377*c5x7*sx3+497*c2x8*sx3+(163-611*x1)*cx3;
X *ds += -12696*sx7*cx3-4200*s2x7*cx3-1503*s3x7*cx3-619*s4x7*cx3-
X 268*s5x7*cx3-(282+1306*x1)*cx7*cx3+(-86+230*x1)*c2x7*cx3+
X 461*s2x8*cx3-350*s2x3+(2211-286*x1)*sx7*s2x3-
X 2208*s2x7*s2x3-568*s3x7*s2x3-346*s4x7*s2x3-
X (2780+222*x1)*cx7*s2x3+(2022+263*x1)*c2x7*s2x3+248*c3x7*s2x3+
X 242*s3x8*s2x3+467*c3x8*s2x3-490*c2x3-(2842+279*x1)*sx7*c2x3;
X *ds += (128+226*x1)*s2x7*c2x3+224*s3x7*c2x3+
X (-1594+282*x1)*cx7*c2x3+(2162-207*x1)*c2x7*c2x3+
X 561*c3x7*c2x3+343*c4x7*c2x3+469*s3x8*c2x3-242*c3x8*c2x3-
X 205*sx7*s3x3+262*s3x7*s3x3+208*cx7*c3x3-271*c3x7*c3x3-
X 382*c3x7*s4x3-376*s3x7*c4x3;
X *ds *= 1e-7;
X
X dp = (7.7108e-2+(7.186e-3-1.533e-3*x1)*x1)*sx5-7.075e-3*sx7+
X (4.5803e-2-(1.4766e-2+5.36e-4*x1)*x1)*cx5-7.2586e-2*cx3-
X 7.5825e-2*sx7*sx3-2.4839e-2*s2x7*sx3-8.631e-3*s3x7*sx3-
X 1.50383e-1*cx7*cx3+2.6897e-2*c2x7*cx3+1.0053e-2*c3x7*cx3-
X (1.3597e-2+1.719e-3*x1)*sx7*s2x3+1.1981e-2*s2x7*c2x3;
X dp += -(7.742e-3-1.517e-3*x1)*cx7*s2x3+
X (1.3586e-2-1.375e-3*x1)*c2x7*c2x3-
X (1.3667e-2-1.239e-3*x1)*sx7*c2x3+
X (1.4861e-2+1.136e-3*x1)*cx7*c2x3-
X (1.3064e-2+1.628e-3*x1)*c2x7*c2x3;
X
X *dm = *dml-(degrad(dp)/s);
X
X *da = 572*sx5-1590*s2x7*cx3+2933*cx5-647*s3x7*cx3+33629*cx7-
X 344*s4x7*cx3-3081*c2x7+2885*cx7*cx3-1423*c3x7+
X (2172+102*x1)*c2x7*cx3-671*c4x7+296*c3x7*cx3-320*c5x7-
X 267*s2x7*s2x3+1098*sx3-778*cx7*s2x3-2812*sx7*sx3;
X *da += 495*c2x7*s2x3+688*s2x7*sx3+250*c3x7*s2x3-393*s3x7*sx3-
X 856*sx7*c2x3-228*s4x7*sx3+441*s2x7*c2x3+2138*cx7*sx3+
X 296*c2x7*c2x3-999*c2x7*sx3+211*c3x7*c2x3-642*c3x7*sx3-
X 427*sx7*s3x3-325*c4x7*sx3+398*s3x7*s3x3-890*cx3+
X 344*cx7*c3x3+2206*sx7*cx3-427*c3x7*c3x3;
X *da *= 1e-6;
X
X *dhl = 7.47e-4*cx7*sx3+1.069e-3*cx7*cx3+2.108e-3*s2x7*s2x3+
X 1.261e-3*c2x7*s2x3+1.236e-3*s2x7*c2x3-2.075e-3*c2x7*c2x3;
X *dhl = degrad(*dhl);
X}
X
X/* ....uranus */
Xstatic
Xp_uranus (t, s, dl, dr, dml, ds, dm, da, dhl)
Xdouble t, s;
Xdouble *dl, *dr, *dml, *ds, *dm, *da, *dhl;
X{
X double dp;
X double x1, x2, x3, x4, x5, x6;
X double x8, x9, x10, x11, x12;
X double sx4, cx4, s2x4, c2x4;
X double sx9, cx9, s2x9, c2x9;
X double sx11, cx11;
X
X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
X
X x8 = mod2PI(1.46205+3.81337*t);
X x9 = 2*x8-x4;
X sx9 = sin(x9);
X cx9 = cos(x9);
X s2x9 = sin(2*x9);
X c2x9 = cos(2*x9);
X
X x10 = x4-x2;
X x11 = x4-x3;
X x12 = x8-x4;
X
X *dml = (8.64319e-1-1.583e-3*x1)*sx9+(8.2222e-2-6.833e-3*x1)*cx9+
X 3.6017e-2*s2x9-3.019e-3*c2x9+8.122e-3*sin(x6);
X *dml = degrad(*dml);
X
X dp = 1.20303e-1*sx9+6.197e-3*s2x9+(1.9472e-2-9.47e-4*x1)*cx9;
X *dm = *dml-(degrad(dp)/s);
X
X *ds = (163*x1-3349)*sx9+20981*cx9+1311*c2x9;
X *ds *= 1e-7;
X
X *da = -3.825e-3*cx9;
X
X *dl = (1.0122e-2-9.88e-4*x1)*sin(x4+x11)+
X (-3.8581e-2+(2.031e-3-1.91e-3*x1)*x1)*cos(x4+x11)+
X (3.4964e-2-(1.038e-3-8.68e-4*x1)*x1)*cos(2*x4+x11)+
X 5.594e-3*sin(x4+3*x12)-1.4808e-2*sin(x10)-
X 5.794e-3*sin(x11)+2.347e-3*cos(x11)+9.872e-3*sin(x12)+
X 8.803e-3*sin(2*x12)-4.308e-3*sin(3*x12);
X
X sx11 = sin(x11);
X cx11 = cos(x11);
X sx4 = sin(x4);
X cx4 = cos(x4);
X s2x4 = sin(2*x4);
X c2x4 = cos(2*x4);
X *dhl = (4.58e-4*sx11-6.42e-4*cx11-5.17e-4*cos(4*x12))*sx4-
X (3.47e-4*sx11+8.53e-4*cx11+5.17e-4*sin(4*x11))*cx4+
X 4.03e-4*(cos(2*x12)*s2x4+sin(2*x12)*c2x4);
X *dhl = degrad(*dhl);
X
X *dr = -25948+4985*cos(x10)-1230*cx4+3354*cos(x11)+904*cos(2*x12)+
X 894*(cos(x12)-cos(3*x12))+(5795*cx4-1165*sx4+1388*c2x4)*sx11+
X (1351*cx4+5702*sx4+1388*s2x4)*cos(x11);
X *dr *= 1e-6;
X}
X
X/* ....neptune */
Xstatic
Xp_neptune (t, s, dl, dr, dml, ds, dm, da, dhl)
Xdouble t, s;
Xdouble *dl, *dr, *dml, *ds, *dm, *da, *dhl;
X{
X double dp;
X double x1, x2, x3, x4, x5, x6;
X double x8, x9, x10, x11, x12;
X double sx8, cx8;
X double sx9, cx9, s2x9, c2x9;
X double s2x12, c2x12;
X
X aux_jsun (t, &x1, &x2, &x3, &x4, &x5, &x6);
X
X x8 = mod2PI(1.46205+3.81337*t);
X x9 = 2*x8-x4;
X sx9 = sin(x9);
X cx9 = cos(x9);
X s2x9 = sin(2*x9);
X c2x9 = cos(2*x9);
X
X x10 = x8-x2;
X x11 = x8-x3;
X x12 = x8-x4;
X
X *dml = (1.089e-3*x1-5.89833e-1)*sx9+(4.658e-3*x1-5.6094e-2)*cx9-
X 2.4286e-2*s2x9;
X *dml = degrad(*dml);
X
X dp = 2.4039e-2*sx9-2.5303e-2*cx9+6.206e-3*s2x9-5.992e-3*c2x9;
X
X *dm = *dml-(degrad(dp)/s);
X
X *ds = 4389*sx9+1129*s2x9+4262*cx9+1089*c2x9;
X *ds *= 1e-7;
X
X *da = 8189*cx9-817*sx9+781*c2x9;
X *da *= 1e-6;
X
X s2x12 = sin(2*x12);
X c2x12 = cos(2*x12);
X sx8 = sin(x8);
X cx8 = cos(x8);
X *dl = -9.556e-3*sin(x10)-5.178e-3*sin(x11)+2.572e-3*s2x12-
X 2.972e-3*c2x12*sx8-2.833e-3*s2x12*cx8;
X
X *dhl = 3.36e-4*c2x12*sx8+3.64e-4*s2x12*cx8;
X *dhl = degrad(*dhl);
X
X *dr = -40596+4992*cos(x10)+2744*cos(x11)+2044*cos(x12)+1051*c2x12;
X *dr *= 1e-6;
X}
X
END_OF_FILE
if test 17647 -ne `wc -c <'plans.c'`; then
echo shar: \"'plans.c'\" unpacked with wrong size!
fi
# end of 'plans.c'
fi
if test -f 'watch.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'watch.c'\"
else
echo shar: Extracting \"'watch.c'\" \(11959 characters\)
sed "s/^X//" >'watch.c' <<'END_OF_FILE'
X/* these functions allow you to watch the sky or the solar system via a
X * simple character-graphics representation on the screen.
X * the interaction starts by using the current time. then control with
X * END returns to table form; or
X * RETURN advances time by one StpSz; or
X * h advances once by 1 hour; or
X * d advances once by 24 hours (1 day); or
X * w advances once by 7 days (1 week); or
X * any other key free runs by StpSz until any key is hit.
X */
X
X#include <stdio.h>
X#include <math.h>
X#include "astro.h"
X#include "circum.h"
X#include "screen.h"
X
X#define SSZCOL 1 /* column to show solar system z coords */
X
X#define DOMESKY 0 /* flags for watch_sky() */
X#define ALTAZSKY 1 /* flags for watch_sky() */
X
X#define SKYACC 3600. /* desired sky plot accuracy, in arc seconds */
X#define SSACC 3600. /* desired solar system plot accuracy, in arc secs */
X
X/* macros to convert row(col) in range 1..NR(1..NC) to fraction in range 0..1 */
X#define r2fr(r) (((r)-1)/(NR-1))
X#define c2fc(c) (((c)-1)/(NC-1))
X#define fr2r(fr) ((int)((fr)*(NR-1)+.5)+1)
X#define fc2c(fc) ((int)((fc)*(NC-1)+.5)+1)
X
X/* single-character tag for each body.
X * order must match the #defines in astro.h and screen.h additions.
X */
Xstatic char body_tags[] = "evmjsunpSMxy";
X
X/* multiple and single loop prompts */
Xstatic char frprompt[] = "Running... press any key to stop.";
Xstatic char qprompt[] =
X"q to quit, RETURN/h/d/w to step by StpSz/hr/day/wk, or any other to freerun";
X
X/* used to locate, record and then erase last plotted chars */
Xtypedef struct {
X double l_fr, l_fc; /* 2d coords as 0..1 (lower left corner is [0,0]) */
X int l_r, l_c; /* screen 2d coords (upper left corner is [1,1]) */
X char l_tag; /* char to use to print on screen */
X} LastDraw;
X
Xstatic int trails; /* !0 if want to leave trails */
X
Xwatch (np, tminc, wbodies)
XNow *np; /* time now and on each step */
Xdouble tminc; /* hrs to increment time by each step */
Xint wbodies; /* each bit is !=0 if want that body */
X{
X static char *flds[4] = {
X "Sky dome", "Alt/az sky", "Solar system"
X };
X static int fn; /* begin with 0, then remember for next time */
X int didone = 0;
X
X while (1) {
X int nf;
X flds[3] = trails ? "Leave trails" : "No trails";
X if ((nf = popup (flds, fn, 4)) < 0)
X break;
X fn = nf;
X switch (nf) {
X case 0: watch_sky (DOMESKY, np, tminc, wbodies); didone = 1; break;
X case 1: watch_sky (ALTAZSKY, np, tminc, wbodies); didone = 1; break;
X case 2: watch_solarsystem (np, tminc, wbodies); didone = 1; break;
X case 3: trails ^= 1; break;
X }
X }
X
X if (didone)
X redraw_screen(2);
X}
X
X/* full alt/az or dome sky view (like the popular astro mags).
X * alt/az: north is at left and right of screen, south at center.
X * 0 elevation is at bottom of screen, zenith at the top.
X * dome: east is left, north is up.
X */
Xstatic
Xwatch_sky (style, np, tminc, wbodies)
Xint style; /* DOMESKY or ALTAZSKY */
XNow *np; /* time now and on each step */
Xdouble tminc; /* hrs to increment time by each step */
Xint wbodies; /* each bit is !=0 if want */
X{
X static char east[] = "East";
X static char west[] = "West";
X static char north[] = "North";
X static char south[] = "South";
X double tminc0 = tminc; /* remember the original */
X /* two draw buffers so we can leave old up while calc new then
X * erase and draw in one quick operation. always calc new in newp
X * buffer and erase previous from lastp. buffers alternate roles.
X */
X LastDraw ld0[NOBJ], ld1[NOBJ], *lp, *lastp = ld0, *newp = ld1;
X int nlast = 0, nnew;
X int once = 1;
X double lmjd, tmp;
X Sky s;
X int p;
X
X /* clear screen and put up the permanent labels */
X c_erase();
X if (style == DOMESKY) {
X f_string (fr2r(.5), fc2c(.5-.5/ASPECT)-7, "East |");
X f_char (fr2r(.5+.5*.707), fc2c(.5-.5*.707/ASPECT)-1, '\\');
X f_string (fr2r(1.), fc2c(.5)-2, south);
X f_char (fr2r(.5+.5*.707), fc2c(.5+.5*.707/ASPECT)+1, '/');
X f_string (fr2r(.5), fc2c(.5+.5/ASPECT)+1, "| West");
X f_char (fr2r(.5-.5*.707), fc2c(.5+.5*.707/ASPECT)+1, '\\');
X f_string (2, NC/2-2, north);
X f_char (fr2r(.5-.5*.707), fc2c(.5-.5*.707/ASPECT)-1, '/');
X } else {
X f_string (NR, 1, north);
X f_string (NR, NC/4, east);
X f_string (NR, NC/2, south);
X f_string (NR, 3*NC/4, west);
X f_string (NR, NC-5, north); /* -1 more to avoid scrolling */
X f_string (2, NC/2-3, "Zenith");
X }
X f_string (2, 1, tznm);
X f_string (3, 1, "LST");
X
X while (1) {
X if (once)
X print_updating();
X
X /* calculate desired stuff into newp[] */
X nnew = 0;
X for (p = nxtbody(-1); p != -1; p = nxtbody(p))
X if (wbodies & (1<<p)) {
X (void) body_cir (p, SKYACC, np, &s);
X if (s.s_alt >= 0) {
X LastDraw *lnp = newp + nnew;
X if (style == DOMESKY) {
X tmp = 0.5 - s.s_alt/PI;
X lnp->l_fr = 0.5 - tmp*cos(s.s_az);
X lnp->l_fc = 0.5 - tmp*sin(s.s_az)/ASPECT;
X } else {
X lnp->l_fr = 1.0 - s.s_alt/(PI/2);
X lnp->l_fc = s.s_az/(2*PI);
X }
X lnp->l_tag = body_tags[p];
X nnew++;
X }
X }
X set_screencoords (newp, nnew);
X
X /* unless we want trails,
X * erase any previous tags (in same order as written) from lastp[].
X */
X if (!trails)
X for (lp = lastp; --nlast >= 0; lp++)
X f_char (lp->l_r, lp->l_c, ' ');
X
X /* print LOCAL time and date we will be using */
X lmjd = mjd - tz/24.0;
X f_time (2, 5, mjd_hr(lmjd));
X f_date (2, 14, mjd_day(lmjd));
X now_lst (np, &tmp);
X f_time (3, 5, tmp);
X
X /* now draw new stuff from newp[] */
X for (lp = newp; lp < newp + nnew; lp++)
X f_char (lp->l_r, lp->l_c, lp->l_tag);
X
X /* swap new and last roles and save new count */
X if (newp == ld0)
X newp = ld1, lastp = ld0;
X else
X newp = ld0, lastp = ld1;
X nlast = nnew;
X
X if (!once)
X slp_sync();
X
X if (once || (chk_char()==0 && read_char()!=0)) {
X if (readwcmd (tminc0, &tminc, &once) < 0)
X break;
X }
X
X /* advance time */
X inc_mjd (np, tminc);
X }
X}
X
X/* solar system view, "down from the top", first point of aries to the right.
X * always include earth.
X */
Xstatic
Xwatch_solarsystem (np, tminc, wbodies)
XNow *np; /* time now and on each step */
Xdouble tminc; /* hrs to increment time by each step */
Xint wbodies;
X{
X /* max au of each planet from sun; in astro.h #defines order */
X static double auscale[] = {.38, .75, 1.7, 5.2, 11., 20., 31., 50.};
X double tminc0 = tminc; /* remember the original */
X /* two draw buffers so we can leave old up while calc new then
X * erase and draw in one quick operation. always calc new in newp
X * buffer and erase previous from lastp. buffers alternate roles.
X */
X LastDraw ld0[2*NOBJ], ld1[2*NOBJ], *lp, *lastp = ld0, *newp = ld1;
X int nlast = 0, nnew;
X int once = 1;
X double lmjd;
X double scale;
X Sky s;
X int p;
X
X /* set screen scale: largest au we will have to plot.
X * never make it less than 1 au since we always show earth.
X */
X scale = 1.0;
X for (p = MARS; p <= PLUTO; p++)
X if ((wbodies & (1<<p)) && auscale[p] > scale)
X scale = auscale[p];
X
X /* clear screen and put up the permanent labels */
X c_erase();
X f_string (2, 1, tznm);
X
X while (1) {
X if (once)
X print_updating();
X
X /* calculate desired stuff into newp[].
X * fake a sun at center and add earth first.
X * (we get earth's loc when ask for sun)
X */
X nnew = 0;
X set_ss (newp+nnew, 0.0, 0.0, 0.0, 'S');
X nnew += 2;
X (void) body_cir (SUN, SSACC, np, &s);
X set_ss (newp+nnew, s.s_edist/scale, s.s_hlong, 0.0, 'E');
X nnew += 2;
X for (p = MERCURY; p <= PLUTO; p++)
X if (p != MOON && (wbodies & (1<<p))) {
X (void) body_cir (p, SSACC, np, &s);
X set_ss (newp+nnew, s.s_sdist/scale, s.s_hlong, s.s_hlat,
X body_tags[p]);
X nnew += 2;
X }
X for (p = OBJX; p != -1; p = (p == OBJX) ? OBJY : -1)
X if (wbodies & (1<<p)) {
X (void) body_cir (p, SSACC, np, &s);
X if (s.s_hlong != NOHELIO && s.s_sdist <= scale) {
X set_ss (newp+nnew, s.s_sdist/scale, s.s_hlong, s.s_hlat,
X body_tags[p]);
X nnew += 2;
X }
X }
X
X set_screencoords (newp, nnew);
X
X /* unless we want trails,
X * erase any previous tags (in same order as written) from lastp[].
X */
X if (!trails)
X for (lp = lastp; --nlast >= 0; lp++)
X safe_f_char (lp->l_r, lp->l_c, ' ');
X
X /* print LOCAL time and date we will be using */
X lmjd = mjd - tz/24.0;
X f_time (2, 5, mjd_hr(lmjd));
X f_date (2, 14, mjd_day(lmjd));
X
X /* now draw new stuff from newp[] */
X for (lp = newp; lp < newp + nnew; lp++)
X safe_f_char (lp->l_r, lp->l_c, lp->l_tag);
X
X /* swap new and last roles and save new count */
X if (newp == ld0)
X newp = ld1, lastp = ld0;
X else
X newp = ld0, lastp = ld1;
X nlast = nnew;
X
X if (!once)
X slp_sync();
X
X if (once || (chk_char()==0 && read_char()!=0)) {
X if (readwcmd (tminc0, &tminc, &once) < 0)
X break;
X }
X
X /* advance time */
X inc_mjd (np, tminc);
X }
X}
X
X/* fill in two LastDraw solar system entries,
X * one for the x/y display, one for the z.
X */
Xstatic
Xset_ss (lp, dist, lg, lt, tag)
XLastDraw *lp;
Xdouble dist, lg, lt; /* scaled heliocentric distance, longitude and lat */
Xchar tag;
X{
X lp->l_fr = 0.5 - dist*sin(lg)*0.5;
X lp->l_fc = 0.5 + dist*cos(lg)*0.5/ASPECT;
X lp->l_tag = tag;
X lp++;
X /* row is to show course helio altitude but since we resolve collisions
X * by adjusting columns we can get more detail by smaller variations
X * within one column.
X */
X lp->l_fr = 0.5 - dist*sin(lt)*0.5;
X lp->l_fc = c2fc(SSZCOL) + (1 - lp->l_fr)/NC;
X lp->l_tag = tag;
X}
X
X/* given a list of LastDraw structs with their l_{fr,fc} filled in,
X * fill in their l_{r,c}.
X * TODO: better collision avoidance.
X */
Xstatic
Xset_screencoords (lp, np)
XLastDraw lp[];
Xint np;
X{
X LastDraw *lpi; /* the current basis for comparison */
X LastDraw *lpj; /* the sweep over other existing cells */
X int i; /* index of the current basis cell, lpi */
X int j; /* index of sweep cell, lpj */
X int n; /* total cells placed so far (ie, # to check) */
X
X /* idea is to place each new item onto the screen.
X * after each placement, look for collisions.
X * if find a colliding pair, move the one with the greater l_fc to
X * the right one cell, then rescan for more collisions.
X * this will yield a result that is sorted by columns by l_fc.
X * TODO: don't just move to the right, try up too for true 2d adjusts.
X */
X for (n = 0; n < np; n++) {
X lpi = lp + n;
X i = n;
X lpi->l_r = fr2r(lpi->l_fr);
X lpi->l_c = fc2c(lpi->l_fc);
X chk:
X for (j = 0; j < n; j++) {
X lpj = lp + j;
X if (i!=j && lpi->l_r == lpj->l_r && lpi->l_c == lpj->l_c) {
X if (lpj->l_fc > lpi->l_fc) {
X /* move lpj and use it as basis for checks now */
X lpi = lpj;
X i = j;
X }
X if (++lpi->l_c > NC)
X lpi->l_c = 1;
X goto chk;
X }
X }
X }
X}
X
X/* since the solar system scaling is only approximate, and doesn't include
X * object x/y at all, characters might get mapped off screen. this funtion
X * guards against drawing chars off screen. it also moves a char being drawn
X * on the lower right corner of the screem left one to avoid scrolling.
X */
Xstatic
Xsafe_f_char (r, c, tag)
Xint r, c;
Xchar tag;
X{
X if (r >= 1 && r <= NR && c >= 1 && c <= NC) {
X if (r == NR && c == NC)
X c -= 1;
X f_char (r, c, tag);
X }
X}
X
X/* see what the op wants to do now and update prompt/times accordingly.
X * return -1 if we are finished, else 0.
X */
Xstatic int
Xreadwcmd (tminc0, tminc, once)
Xdouble tminc0;
Xdouble *tminc;
Xint *once;
X{
X f_prompt (qprompt);
X
X switch (read_char()) {
X case END: /* back to table */
X return (-1);
X case '\r': /* one StpSz step */
X *tminc = tminc0;
X *once = 1;
X break;
X case 'h': /* one 1-hour step */
X *tminc = 1.0;
X *once = 1;
X break;
X case 'd': /* one 24-hr step */
X *tminc = 24.0;
X *once = 1;
X break;
X case 'w': /* 7 day step */
X *tminc = 7*24.0;
X *once = 1;
X break;
X default: /* free-run */
X *once = 0;
X f_prompt (frprompt);
X }
X return (0);
X}
END_OF_FILE
if test 11959 -ne `wc -c <'watch.c'`; then
echo shar: \"'watch.c'\" unpacked with wrong size!
fi
# end of 'watch.c'
fi
echo shar: End of archive 4 \(of 6\).
cp /dev/null ark4isdone
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