allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc) (06/04/89)
Posting-number: Volume 7, Issue 9
Submitted-by: ecd@ncs-med.UUCP (Elwood C. Downey)
Archive-name: ephem/part01
[Manual uses -ms. May I suggest that the only formatter macros you can be
*sure* exist are the ones in -man? -ms and -mm are Berkeleyisms, -me is
an AT&Tism. Maybe you should use -mn, which comes with News 2.11. ++bsa]
Here is an ephemeris program I have written I thought I'd toss to the wind.
I started it a few years ago and now it is pretty stable and I know of
no bugs in it. It uses termcap on unix but all the terminal i/o is based
on a very few simple functions in io.c, the intent is that it could be
easily ported to DOS or other non-unix environment; even the unix version
uses no other unix-ism.
I have gathered everything together into three sh archive files; this
is the first. Cut this chit-chat off and run the script to create the first
set of files. Do the same for the others. When it's done, you should have
these files:
Makefile # makes ephem (given lib/lib.a is done)
Man.ms # nroff -ms manual
TODO # stuff I've thought to add someday
ephem.cfg # sample configuration file
borders.c
circum.c
circum.h
io.c
main.c
plot.c
pr.c
pr0.c
screen.h
sel_lfd.c
time.c
version.c
lib/Makefile # makes the pure astronomy stuff, lib.a
lib/aa_hadec.c
lib/anomaly.c
lib/astro.h
lib/cal_mjd.c
lib/eq_ecl.c
lib/moon.c
lib/moonnf.c
lib/moonrs.c
lib/nutation.c
lib/obliq.c
lib/parallax.c
lib/pelement.c
lib/plans.c
lib/precess.c
lib/refract.c
lib/riset.c
lib/sex_dec.c
lib/sun.c
lib/sunrs.c
lib/utc_gst.c
(the third script makes the lib directory).
When it's done, make it by making the library, then the main program::
cd lib
make
cd ..
make
then try to run it by just typing:
ephem
Let me know about suggestions, bugs, problems, etc. I don't guarantee I'll
maintain it particularly regularly but I'm enough of an astro nut that I will
likely fool with it some more without much provoking. Enjoy!
Elwood Downey
umn-cs!ncs-med!ecd
---CUT---CUT---CUT---CUT---CUT---CUT---CUT---CUT---CUT---CUT---CUT---CUT---CUT
#!/bin/sh
echo extracting ephem.cfg
cat > ephem.cfg << 'xXx'
UT NOW
LONG 93:42:8
LAT 44:50:37
HEIGHT 800
TEMP 40
PRES 29.5
STPSZ RTC
PROPTS TSMevajsunp
EPOCH EOD
SITE The Hill Observatory
OBJN Orion
OBJRA 6
OBJDEC 0
NSTEP 1
xXx
echo extracting Man.ms
cat > Man.ms << 'xXx'
.LP
.ND
.po .5i
.nr LL 7.4i \" line length
.nr LT \n(LL \" title length
.nr FL \n(LL \" footnote length
.nr PO .5i \" page offset
.nh
.nr Pd \n(PD \" save paragraph spacing for changing .IP
.na
.LP
.DS C
Ephem
V3.11
by
Elwood Downey
8860 Abbywood Road
Chaska, MN 55343
.DE
.SH
Introduction
.LP
Ephem.exe is a program that displays observing circumstances
for all the planets and any one extra object you wish to enter. Included are
RA and Dec, heliocentric coordinates, azimuth and altitude, distance from
sun and earth, phase, visual magnitude, solar elongation, and angular size.
It also displays local ephemeris information. Included are UTC and local
date and time, local sidereal time, local sun and moon rise and set times,
times of astronomical twilight, length of day and night,
and a monthly calendar.
RA/Dec calculations are geocentric and include the effects of precession,
nutation and aberration.
Alt/az calculations are topocentric and also include effects of refraction and
parallax.
A running plot file of selected field values
may be generated as the program runs.
Ephem.exe includes a very crude
quick-look capability to view these plot files or they may
be plotted by other programs.
The program is written in C for unix or DOS. It uses only a very simple
set of io routines and should be easily ported to any ASCII display.
The DOS version requires the ANSI.SYS screen driver.
The planetary data and correction algorithms are taken, with permission,
from "Astronomy With Your Personal Computer",
by Peter Duffett-Smith, Cambridge University Press, 1985.
.bp
.SH
Sample Display
.LP
Here is a typical screen from
.I ephem:
.DS L
The Hill Observatory
Move to another field, RETURN to change this field, ? for help, or ESC to run
May 1989
JD 2447664.08208 LST 23:24:22 TZ 5:00:00 Su Mo Tu We Th Fr Sa
UTC 13:58:12 5/17/1989 Lat 44:50:37 Long 93:42:08 1 2 3 4 NM 6
CDT 8:58:12 5/17/1989 Dawn 3:44 Dusk 22:40 7 8 9 10 11 12 13
SRis 5:44 @ 61:10 StpSz RT CLOCK NStep 1 14 15 16 17 18 19 FM
SSet 20:39 @ 299:03 Elev 800 ft Temp 40 F 21 22 23 24 25 26 27
MRis 17:58 @ 109:36 AtmPr 29.50 in DayLn 14:55 28 29 30 31
MSet 3:54 @ 254:42 Plot off NiteLn 5:04
Ob R.A. Dec Helio Helio Az Alt Ea Dst Sn Dst Elong Size VMag Phs
(Epoch of date)Long Lat Deg E Deg Up AU(mi) AU Deg E ArcS %
Su 3:37.3 19:24 236:38 0:00 94:20 32:24 1.0114 1898 -27
Mo 13:12.1 -11:48 316:42 -50:18 251450 1.0136 144.3 1772 -12 80
Me 4:15.7 21:37 224:39 0:26 85:39 27:04 0.5804 0.4484 9.2 11.6 0.3 3
Ve 4:23.8 21:47 83:30 0:25 84:08 25:44 1.6854 0.7200 11.1 10.0 -4.0 98
Ma 6:50.0 24:18 127:03 1:48 58:11 3:21 2.2010 1.6450 44.8 4.3 1.3 95
Ju 4:47.6 21:57 76:32 -0:32 80:03 21:40 6.0333 5.0725 16.6 32.6 -2.0 100
Sa 18:58.1 -22:06 279:14 0:37 238:13 0:24 9.3194 10.039 -133.2 17.8 1.0 100
Ur 18:20.7 -23:37 272:54 -0:15 243:35 -7:09 18.533 19.339 -141.9 3.6 5.6 100
Ne 18:52.2 -21:58 280:43 0:55 239:21 -0:09 29.498 30.216 -134.5 2.1 7.9 100
Pl 15:03.1 -0:16 223:51 15:46 296:22 -24:25 28.710 29.657 159.1 0.3 13.7 100
Or 6:00.0 0:00 83:42 -6:18
.DE
.SH
Program Operation
.LP
When
.I ephem
is started, it first displays a disclaimer banner.
Then, after any key is depressed,
it reads a configuration file
to set the initial values of several fields.
This file, ephem.cfg by default, is described below.
It then draws all fields on the screen with their initial values.
The program then loops advancing time each step, by some amount you may
control, and updating all fields accordingly.
There are two fields that control this looping behavior: NStep and StpSz.
These control the number of steps and the amount of time
to add each step, respectively.
When the number of steps, NStep, goes to 0 or any key is depressed,
the looping stops and you enter
.I "command mode."
Command mode allows you to modify most of the fields.
The idea is that you move to each field on the screen you wish to change
and change them.
When you have changed everything you want to,
type ESC to update the other fields on the screen.
To change a field: move the cursor to the field;
type RETURN; then
type in the new value along the command line at the top according to
the format indicated in the prompt. To accept the new value type RETURN,
or to leave it unchanged after all type ESC.
A few fields don't require you to type anything; just typing RETURN does
all the work.
If you can't move to it, you can't change it.
The arrow keys on most systems move the cursor around.
If these do not function or function incorrectly,
the h/j/k/l keys also move the cursor left/down/up/right, respectively.
When you have changed a field that would invalidate any of the other fields
the message NEW CIRCUMSTANCES appears in the upper right of the screen.
This will remain until you
type ESC to allow at least one screen update loop to occur.
If you change any time field, the StpSz value is not added to the first loop.
Note also that after a series of loops, NStep is automatically reset to 1
so ESC will do exactly one loop and return you to command mode.
To quit the program, type control-x from command mode.
To show a simple help line, type ? any time.
.bp
.SH
Screen Fields
.LP
These are the fields displayed by the program.
Following each name, in parentheses, might be
"c" to mean the field may be picked to be changed or
"p" to mean the field may be picked to be included in a plot (see below).
.nr PD 0
.IP JD 14
the current Julian date, to about 1-second accuracy.
.IP LST(p)
the current local sidereal time.
.IP LT(c)
the local timezone name.
This field can be changed to any three-character mnemonic for the current time
zone name.
.IP UTC(cp)
universally coordinated time. set to NOW to set from computer clock.
.IP TZ(cp)
hours local time is behind utc, ie, >0 west, <0 east of Greenwich.
.IP Lat(cp)
location latitude, positive degrees north of equator.
.IP Long(cp)
location longitude, positive degrees west of greenwich meridian.
.IP Temp(cp)
local surface air temperature, in degrees F.
.IP AtmPr(cp)
local surface air pressure, in inches of mercury.
.IP NStep(c)
The number of times the display with be updated (time advanced by StpSz each
step) before entering command mode.
.IP StpSz(c)
the amount of time UTC (and its derivatives) is incremented each loop.
set this to RTC to use real-time based on the computer clock.
you may also set it in terms of days by appending a D (or d)
after the number when you set it.
.IP Elev(cp)
local elevation of the ground above sea level, in feet.
.IP Dusk(cp)
local time when the sun is about 15 degrees below the horizon after sunset.
.IP SRis(cp)
local time when the sun upper limb appears, ie, sunrise.
.IP MRis(cp)
Local time when the moon upper limb appears, ie, moonrise.
.IP Dawn(cp)
local time when the sun is about 18 degrees below the horizon before sunrise
.IP SSet(cp)
local time when the sun upper limb disappears, ie, sunset.
.IP MSet(cp)
local time when the moon upper limb disappears, ie, moonset.
.IP DayLn(cp)
length of time sun is above horizon, ie, SSet - SRis.
.IP NiteLn(cp)
length of astronical night, ie, Dawn - Dusk.
.IP Plot(c)
controls plotting; see separate discussion, below.
.nr PD \(Pd
In the upper right of the screen is a calendar for the current local month.
Dates of new and full moons are marked NM and FM, respectively.
.LP
Some things may be turned off to reduce compute times.
Each planet may be turned on and off
by selecting the planet name field.
The calculation of dawn/dusk, sunrise/set, and moonrise/set may each be
turned on and off by selecting their fields (note
that even when they are turned
on the software only recalculates them when the local date changes.)
.bp
.LP
The planets are displayed in a table in the bottom portion of the screen.
There is one row per planet, and several columns, described next.
.LP
One object may also be added to the display by putting its name and location
in the bottom row of the screen. This is done by moving the cursor to
the bottom row and selecting the field under the Ob, R.A., and Dec columns.
.nr PD 0
.IP Ob 14
name of object.
.IP R.A.(p)
right ascension of object, precessed to given epoch.
.IP Dec(p)
declination of object, precessed to given epoch.
.IP Helio Long(p)
heliocentric longitude. the earth's is displayed on the Sun's line.
.IP Helio Lat(p)
heliocentric latitude.
.IP Az(p)
degrees eastward of true north for object.
.IP Alt(p)
degrees up from a horizontal plane Elev feet above sea level.
.IP "Ea Dst(p)"
distance from earth center to object center, in AU, except distance
to moon is in miles.
.IP "Sn Dst(p)"
distance from sun center to object center, in AU.
.IP Elong(p)
spherical angular separation between sun and object. note this
is not just difference in ecliptic longitude. the sign, however, is
simply sign(obj_long - sun_long), ie, degrees east. thus, a positive
elongation means the object rises after the sun.
.IP Size(p)
angular size of object, in arc seconds.
.IP VMag(p)
visual magnitude of object.
.IP Phs(p)
percent of visible surface in sunlight, ie, the phase. note the
moon phase is calculated simplistically as just abs(elongation)/180*100
which can be a few degrees off... this means that because of how
elongation is defined it doesn't say 0 during new moon (or 100 during full)
except during close eclipses (maybe that's a "feature"?).
.nr PD \n(Pd
The precession epoch is displayed beneath the RA/DEC column headings.
.bp
.SH
Date and Time Formats.
.LP
Times are displayed and entered in h:m:s format.
Any of the h, m, and s components that are not specified are left unchanged
from their current value.
For example, 0:5:0 set hours to 0, minutes to 5, seconds to 0, whereas :5
sets minutes to 5 but leaves hours and seconds unchanged.
A negative time is indicated by
a minus sign (-) anywhere before the first digit.
Dates are displayed and entered in American m:d:y format.
As with time,
components omitted remain the current value.
For example, if the current date is 10/20/1988 and you type 20/20 the new
date will become 20/20/1988. Note you must type the full year since the
program is accurate over several centuries either side of 1900.
.SH
Configuration File
.LP
The ephem.cfg configuration file allows you to set the initial values of
many of the screen fields. You can still change any field while the program
is running too; this file just sets the initial conditions.
You can have several different configuration files if you wish. By
default,
.I ephem
looks for one named ephem.cfg. You can tell it to use
an alternate file by using the -c switch as follows:
.DS L
ephem -c <filespec>
.DE
The format of the file uses the form KEYWORD=VALUE, where the possible
KEYWORDS and the types of VALUES for each are described below. Any KEYWORDS
not in the file will take on some sort of default.
Note:
because of the way unspecified time and date components are left unchanged
(see section on Date and Time Formats) always specify the complete time and
date for all entries in the configuration file. For example, to initialize
the longitude to zero degrees, say 0:0:0, not just 0.
.nr PD 0
.IP UD 10
initial UTC date, such as 10/20/1988 or NOW to use the computer clock.
.IP TZONE
hours the local time is behind utc, such as 5:0:0.
you need not set this if you use NOW for UT or UD.
.IP TZNAME
name of the local time zone, such as CDT. 3 chars max.
you need not set this if you use NOW for UT or UD.
.IP LONG
longitude, in degrees west of greenwich, in the form d:m:s.
.IP LAT
latitude, in degrees north of the equator, in the form d:m:s.
.IP HEIGHT
height above sea level, in feet, such as 800
.IP TEMP
air temperature, in degrees F, such as 50
.IP PRES
air pressure, in inches of Mercury, such as 29
.IP STPSZ
the time increment between screen updates, such as "1" to give one
hour updates. this can be a specific amount or RTC to use the system
clock as a real-time source. You may also specify a time in days, by
appending a D (or d) after the number.
.IP PROPTS
this selects what you want included in the display. since IBM-PC math
is not very fast, you can reduce the time to update the screen
by only printing those fields of interest. the VALUE is a collection
of letters to turn on each item from the following set:
.DS L
T twilight (dawn-dusk)
S sun rise/set and circumstance for the sun
M moon rise/set and circumstance for the moon
e circumstances for mercury
v circumstances for venus
a circumstances for mars
j circumstances for jupiter
s circumstances for saturn
u circumstances for uranus
n circumstances for neptune
p circumstances for pluto
.DE
For example, to just display sun rise/set and track the sun and
saturn, say PROPTS Ss
.IP NSTEP
number of times program will loop before entering command mode.
see the discussion under Program Operation.
.IP EPOCH
this sets the desired ra/dec precession epoch. you can put any date
here or EOD to use the current instant ("Epoch of Date").
.IP SITE
the name of the observing site. it places the characters found here
across the top of the screen as a title.
.IP OBJN
name of the extra object to track.
.IP OBJRA
right ascension of the extra object to track.
.IP OBJDEC
declination of the extra object to track.
.nr PD \n(Pd
.bp
.LP
Example ephem.cfg files:
.DS L
Create an essentially free-running real-time update:
UT NOW
LONG 90:10:8
LAT 40:50:20
HEIGHT 800
TEMP 50
PRES 29
STPSZ RTC
PROPTS TSMevajsunp
NSTEP 10000000
EPOCH EOD
SITE The Observatory
.DE
.DS L
Initialize things to investigate the 1991 Hawaiian solar eclipse:
UD 7/11/1991
UT 19:10:0
LAT 20:30:0
LONG 157:0:0
TZONE 10
TZNAME HST
EPOCH 2000
PRES 30
HEIGHT 0
TEMP 80
SETPSZ :10
NSTEP 1
PROPTS SM
SITE 1991 Hawaii Eclipse
.DE
.bp
.SH Plotting
.LP
Each time a field is drawn on the screen its full-precision
value may be written to a file.
Each line in the file consists of a tag character followed by two or three
floating point variables, all separated by commas. If there are two
values, they should be interpreted to be x and y (or perhaps r and theta).
If there is a third, it is a z or trace value.
.LP
The "Plot" field controls plotting.
Whether plotting is currently active is indicated by "on" or "off" immediately
to its right.
.LP
To initiate a plot file, pick the "off" field. You will be asked for the
name of the file to use and, if it already exists, whether to overwrite it or
append to it. Once you have chosen a file, plotting is on and the field changes
to "on".
.LP
You should then define which fields should be plotted. Select the "Plot"
field. You will be asked
for a tag character, then asked
to move the cursor to the field you want to use as the x coordinate (abscissa),
then asked to choose the y coordinate (ordinate), then asked to choose an
optional z trace variable.
You may choose up to four of these sets for any given plot run. Type ESC to
quick making choices and begin plotting. The values are written to the plot
file each time they are updated on the screen until you select "on" to turn
plotting back off.
.LP
You may use
.I ephem
to make a crude plot of the plot files. When plotting is
off, select the "Plot" field and give a filename. The entries will be
drawn with their tag characters; the plot remains on the screen until you
type any character.
.bp
.SH
Implementation Notes
.LP
All rise/set times are for the current local date.
However, if the event occurred just before midnight this morning
the time reported might be for the previous day, and similarly after tonights
midnight. If in doubt, set the time and check the altitude.
The calendar is for current local month.
The program uses a horizontal plane tangent to the earth as the
horizon for all altitude calculations, rise/set events, etc.
This is
.I not
the same as the angle up from the local horizon unless the observer is
directly on the ground due to earth's curvature.
The effect can be found from:
.DS L
sin(a)**2 = (h**2 + 2Rh) / (R+h)**2
where:
R = radius of earth
h = height above ground (same units as R)
a = increase in altitude
.DE
For example, the effect is more than two arc minutes at a height of 5 feet.
visual magnitudes are not very accurate at all... haven't bother to fix.
internally, the program is good to a few tens of arc seconds on the planets;
displaying to nearest arc minute should be quite reliable.
a "negative 0" is displayed when a value is negative but less than half the
display precision.
The sun-moon distance is the solution for the third side of a planar triangle
whose two other sides are the earth-moon distance and earth-sun
distance separated by the angle of elongation.
Sun
rise and set times and azimuths could be found from an iteration directly on
altitude. However, the "standard" algorithm is based on a fixed, average
refraction correction at sea level (elevation makes a slight effect on moon
times). This matches published values pretty well.
The "adaptive" algorithm takes into
account the air pressure and temperture and sun diameter but I have
no independent means to check (eg, an ocean).
Anyway, I figure it isn't worth too much effort.
.bp
.SH
DOS Installation Procedure
.LP
Summary:
You must be running DOS V2.0 or later.
A 8087 floating point chip will be used if present.
The distribution floppy contains two files, ephem.exe and ephem.cfg.
Ephem.exe is the executable program; ephem.cfg is a sample configuration
file. To run the program, make working copies of these two files
in a directory and run "ephem" from that directory. The program uses the
ANSI.SYS terminal driver for screen control. It also uses an
environment variable, TZ, to establish the local timezone.
Details:
.DS L
1) The ANSI.SYS screen driver is required for this program. Edit the
CONFIG.SYS file, if necessary, so it contains the following line:
ANSI.SYS
If it wasn't already there and you had to add it, note it will not
take effect until you reboot DOS.
2) Set a DOS environment variable, TZ, in the following form:
set TZ=SSSnDDD
"SSS" is the 3-letter abbreviation for the local standard timezone;
"n" is a number between -23 to 24 indicating the number of hours
that are subtracted from GMT to obtain local standard time;
"DDD" is an optional 3-letter abbreviation for the local daylight savings
time zone name. Leave it off if you do not have savings time in your
area. If the changeover dates differ from the internal algorithm,
just use SSS and n directly.
For example, in the midwestern United States with savings times:
set TZ=CST6CDT
You can put this in your AUTOEXEC.BAT file so it gets set each time
you boot DOS.
2) place the distribution floppy into drive a:.
3) copy the ephem.* files to a working diskette:
copy ephem.* b:*.*/v
or hard disk:
copy ephem.* c:*.*/v
4) run using the sample configuration file by just running
ephem
To run with a different configuration file, use the -c switch:
ephem -c <filespec>
.DE
xXx
echo extracting TODO
cat > TODO << 'xXx'
allow for plotting backwards: times down, ra left
add facility for finding max alt on a given day. useful for plotting
planet sky positions over many days.
rise, set times for each planet?
new scheme for showing more info/planet than can fit on one line.
for example, make a new pick that allows selection of what you DO want
displayed, up to the maximum allowable columns. add all rise/set, times,
the new max alt, etc
a real-time plot mode. select fields to plot and update.
could watch planets revolve around sun, or watch them go across the sky.
xXx
echo extracting Makefile
cat > Makefile << 'xXx'
CLNFLAGS=$(CLNF)
LNFLAGS=$(CLNFLAGS) $(LNF)
CFLAGS=$(CLNFLAGS) $(CF) -Ilib
LIBRARIES=lib/lib.a /usr/lib/libtermcap.a
LINTFLAGS=$(CFLAGS) $(LINTF)
LINTLIBS=
EPHEM= borders.o \
circum.o \
io.o \
main.o \
plot.o \
pr.o \
pr0.o \
sel_fld.o \
time.o \
version.o \
$(LIBRARIES)
ephem: $(EPHEM)
$(CC) -o $@ $(LNFLAGS) $(EPHEM) $(LNTAIL)
xXx
echo extracting borders.c
cat > borders.c << 'xXx'
#include "screen.h"
borders()
{
/*
register i;
for (i = 1; i <= 80; i++)
pr_char (R_PLANTAB-1, i, '-');
for (i = R_LST; i < R_PLANTAB-1; i++)
pr_char (i, C_LST-1, '|');
for (i = R_TZONE; i < R_PLANTAB-1; i++)
pr_char (i, C_TZONE-1, '|');
for (i = R_CAL+1; i < R_PLANTAB-1; i++)
pr_char (i, C_CAL-1, '|');
*/
}
xXx
echo extracting circum.c
cat > circum.c << 'xXx'
#include <stdio.h>
#include <math.h>
#include "astro.h"
#include "circum.h"
/* shorthands into np */
#define mjd np->n_mjd
#define lat np->n_lat
#define lng np->n_lng
#define tz np->n_tz
#define temp np->n_temp
#define pressure np->n_pressure
#define height np->n_height
/* find sun's circumstances now */
sun_cir (np, sp)
Now *np;
Sky *sp;
{
static Sky last_sky;
static Now last_now;
float lst, alt, az;
if (same_cir (np, &last_now) && about_now (np, &last_now, 7e-4))
*sp = last_sky;
else {
float lsn, rsn;
float deps, dpsi;
last_now = *np;
sun ((float)mjd,&lsn,&rsn); /* sun's true ecliptic long and dist */
nutation ((float)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_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 ((float)mjd, 0.0, lsn, &sp->s_ra, &sp->s_dec);
}
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;
}
/* find moon's circumstances now */
moon_cir (np, sp)
Now *np;
Sky *sp;
{
static Sky last_sky;
static Now last_now;
static float ehp;
float lst, alt, az;
float ha, dec;
if (same_cir (np, &last_now) && about_now (np, &last_now, 2e-4))
*sp = last_sky;
else {
float lam, bet;
float deps, dpsi;
float lsn, rsn; /* sun long in rads, earth-sun dist in au */
float edistau; /* earth-moon dist, in au */
float el; /* elongation, rads east */
last_now = *np;
moon ((float)mjd,&lam,&bet,&ehp); /* moon's true ecliptic loc */
nutation ((float)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 ((float)mjd, bet, lam, &sp->s_ra, &sp->s_dec);
sun ((float)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;
}
/* 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;
}
/* find planet p's circumstances now */
planet_cir (p, np, sp)
int p;
Now *np;
Sky *sp;
{
typedef struct {
float l_every; /* recalc l_sky every these days */
Now l_now; /* when l_sky was found */
Sky l_sky;
} Last;
/* set to number of days planet takes to move about 1 arcsec */
static Last last[8] =
{{2e-4},{4e-4},{1e-3},{.006},{.02},{.06},{.1},{.2}};
float lst, alt, az;
register Last *lp = last + p;
/* if less than l_every days from last time for this planet
* just redo alt/az.
*/
if (same_cir(np, &lp->l_now) && about_now (np, &lp->l_now, lp->l_every))
*sp = lp->l_sky;
else {
float lpd0, psi0; /* heliocentric ecliptic longitude and latitude */
float rp0; /* dist from sun */
float rho0; /* dist from earth */
float lam, bet; /* geocentric ecliptic long and lat */
float dia, mag; /* angular diameter at 1 AU and magnitude */
float lsn, rsn; /* true geoc lng of sun, dist from sn to earth*/
float deps, dpsi;
float a, ca, sa;
float el; /* elongation */
float f; /* phase from earth */
lp->l_now = *np;
plans ((float)mjd,p,&lpd0,&psi0,&rp0,&rho0,&lam,&bet,&dia,&mag);
nutation ((float)mjd, &deps, &dpsi); /* correct for nutation */
lam += dpsi;
sun ((float)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 ((float)mjd, bet, lam, &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;
}
/* 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;
}
/* find s_ra/dec/alt/az @ EOD for object at given loc */
obj_cir (ra, dec, e, np, sp)
float ra, dec, e; /* objects location and epoch of coords */
Now *np;
Sky *sp;
{
float lst, alt, az;
/* always want EOD ra/dec in Sky */
sp->s_ra = ra;
sp->s_dec = dec;
if (e != mjd)
precess (e, (float)mjd, &sp->s_ra, &sp->s_dec);
/* find alt/az based on EOD */
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;
}
/* find times when sun is 18 degrees below horizon */
twilight_cir (np, dawn, dusk, status)
Now *np;
float *dawn, *dusk;
int *status;
{
static Now last_now;
static float last_dawn, last_dusk;
static int last_status;
if (same_cir (np, &last_now) && same_lday (np, &last_now)) {
*dawn = last_dawn;
*dusk = last_dusk;
*status = last_status;
} else {
float tmp;
sunrs ((float)(mjd-tz/24.0), lat, lng, degrad(18.), dawn, dusk,
&tmp, &tmp, status);
last_dawn = *dawn;
last_dusk = *dusk;
last_status = *status;
last_now = *np;
}
}
sunrs_cir (np, dis, utcr, utcs, azr, azs, status)
Now *np;
float dis;
float *utcr, *utcs, *azr, *azs;
int *status;
{
static Now last_now;
static float last_r, last_s, last_azr, last_azs, last_dis;
static int last_status;
if (same_cir (np, &last_now) && same_lday (np, &last_now)
&& last_dis == dis) {
*utcr = last_r;
*utcs = last_s;
*azr = last_azr;
*azs = last_azs;
*status = last_status;
} else {
last_dis = dis; /* save before we modify it in place if ADPREF */
if (dis == ADPREF) {
/* use the real sun diameter and current refraction conditions.
* unrefract the sun upper limb, then subtract sun semi-diam.
*/
Sky sk;
unrefract (pressure, temp, 0.0, &dis);
sun_cir (np, &sk);
dis -= degrad(sk.s_size/3600./2.0);
dis = -dis;
}
sunrs ((float)(mjd-tz/24.0), lat, lng, dis, utcr, utcs, azr, azs,
status);
last_r = *utcr;
last_s = *utcs;
last_azr = *azr;
last_azs = *azs;
last_status = *status;
last_now = *np;
}
}
moonrs_cir (np, utcr, utcs, azr, azs, status)
Now *np;
float *utcr, *utcs, *azr, *azs;
int *status;
{
static Now last_now;
static float last_r, last_s, last_azr, last_azs;
static int last_status;
if (same_cir (np, &last_now) && same_lday (np, &last_now)) {
*utcr = last_r;
*utcs = last_s;
*azr = last_azr;
*azs = last_azs;
*status = last_status;
} else {
moonrs ((float)(mjd-tz/24.0), lat, lng, utcr, utcs, azr, azs,
status);
last_r = *utcr;
last_s = *utcs;
last_azr = *azr;
last_azs = *azs;
last_status = *status;
last_now = *np;
}
}
/* 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.
*/
static
elongation (lam, bet, lsn, el)
float lam, bet, lsn;
float *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. */
static
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);
}
/* return whether the two Nows are for the same LOCAL day */
static
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;
float dt;
{
return (fabs (n1->n_mjd - n2->n_mjd) < dt);
}
now_lst (np, lst)
Now *np;
float *lst;
{
utc_gst ((float)mjd_day(mjd), (float)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 extracting circum.h
cat > circum.h << 'xXx'
#define SPD (24.0*3600.0) /* seconds per day */
#define EOD (-9876) /* special epoch flag: use epoch of date */
#define RTC (-1234) /* special tminc flag: use rt clock */
#define ADPREF (-100) /* special sunrs dis flag: use pres/temp */
/* info about our local observing circumstances */
typedef struct {
double n_mjd; /* modified Julian date, ie, days since
* Jan 0.5 1900 (== 12 noon, Dec 30, 1899), utc.
* enough precision to get well better than 1 second.
*/
float n_lat; /* latitude, >0 north, rads */
float n_lng; /* longitude, >0 east, rads */
float n_tz; /* time zone, hrs behind UTC */
float n_temp; /* atmospheric temp, degrees C */
float n_pressure; /* atmospheric pressure, mBar */
float n_height; /* height above sea level, earth radii */
char n_tznm[4]; /* time zone name; 3 chars or less, always 0 at end */
} Now;
extern double mjd_day(), mjd_hr();
/* info about where and how we see something in the sky */
typedef struct {
float s_ra; /* ra, rads (equinox of date) */
float s_dec; /* dec, rads (equinox of date) */
float s_az; /* azimuth, >0 e of n, rads */
float s_alt; /* altitude above topocentric horizon, rads */
float s_sdist; /* dist from object to sun, au */
float s_edist; /* dist from object to earth, au */
float s_elong; /* angular sep between object and sun, >0 if east */
float s_hlong; /* heliocentric longitude, rads */
float s_hlat; /* heliocentric latitude, rads */
float s_size; /* angular size, arc secs */
float s_phase; /* phase, % */
float s_mag; /* visual magnitude */
} Sky;
xXx
echo extracting io.c
cat > io.c << 'xXx'
/* this file (in principle) contains all the device-dependent code for
* handling screen movement and reading the keyboard. public routines are:
* c_pos(r,c), c_erase(), c_eol();
* chk_char(), read_char(), read_line (buf, max); and
* byetty().
* N.B. we assume output may be performed by printf() and putchar(). output
* buffering has been disabled elsewhere as setbuf((char *)0).
* #define UNIX or LATTICE_C to give two popular versions.
* UNIX uses termcap; LATTICE_C uses ANSI and the IBM-PC keyboard codes.
*/
#define UNIX
/* #define LATTICE_C */
#include <stdio.h>
#define CNTRL(c) ((c)&037)
#define ESC CNTRL('[')
#ifdef UNIX
#include <sgtty.h>
#include <signal.h>
extern char *tgoto();
static char *cm, *ce, *cl, *kl, *kr, *ku, *kd; /* curses sequences */
static int tloaded;
static int ttysetup;
static struct sgttyb orig_sgtty;
/* move cursor to row, col, 1-based.
* we assume this also moves a visible cursor to this location.
*/
c_pos (r, c)
int r, c;
{
if (!tloaded) tload();
fputs (tgoto (cm, c-1, r-1), stdout);
}
/* erase entire screen. */
c_erase()
{
if (!tloaded) tload();
fputs (cl, stdout);
}
/* erase to end of line */
c_eol()
{
if (!tloaded) tload();
fputs (ce, stdout);
}
/* return 0 if there is a char that may be read without blocking, else -1 */
chk_char()
{
long n = 0;
if (!ttysetup) setuptty();
ioctl (0, FIONREAD, &n);
return (n > 0 ? 0 : -1);
}
/* read the next char, blocking if necessary, and return it. don't echo.
* map the arrow keys if we can too into hjkl
*/
read_char()
{
char c;
if (!ttysetup) setuptty();
read (0, &c, 1);
return (chk_arrow (c & 0177)); /* just ASCII, please */
}
/* used to time out of a read */
static got_alrm;
static
on_alrm()
{
got_alrm = 1;
}
/* see if c is the first of any of the curses arrow key sequences.
* if it is, read the rest of the sequence, and return the hjkl code
* that corresponds.
* if no match, just return c.
*/
static
chk_arrow (c)
register char c;
{
register char *seq;
if (c == *(seq = kl) || c == *(seq = kd) || c == *(seq = ku)
|| c == *(seq = kr)) {
char seqa[10]; /* maximum arrow escape sequence ever expected */
int l = strlen(seq);
seqa[0] = c;
/* most arrow keys generate sequences starting with ESC. if so
* c might just be a lone ESC; time out if so.
*/
got_alrm=0;
if (c == ESC) {
signal (SIGALRM, on_alrm);
alarm(1);
}
read (0, seqa+1, l-1);
if (got_alrm == 0) {
if (c == ESC)
alarm(0);
seqa[l] = '\0';
if (strcmp (seqa, kl) == 0)
return ('h');
if (strcmp (seqa, kd) == 0)
return ('j');
if (strcmp (seqa, ku) == 0)
return ('k');
if (strcmp (seqa, kr) == 0)
return ('l');
}
}
return (c);
}
/* do whatever might be necessary to get the screen and/or tty back into shape.
*/
byetty()
{
ioctl (0, TIOCSETP, &orig_sgtty);
}
static
tload()
{
extern char *getenv(), *tgetstr();
extern char *UP, *BC;
char *egetstr();
static char tbuf[512];
char rawtbuf[1024];
char *tp;
char *ptr;
if (!(tp = getenv ("TERM"))) {
printf ("Must have addressable cursor\n");
exit(1);
}
if (!ttysetup) setuptty();
if (tgetent (rawtbuf, tp) != 1) {
printf ("can't find termcap for %s\n", tp);
exit (1);
}
ptr = tbuf;
ku = egetstr ("ku", &ptr);
kd = egetstr ("kd", &ptr);
kl = egetstr ("kl", &ptr);
kr = egetstr ("kr", &ptr);
cm = egetstr ("cm", &ptr);
ce = egetstr ("ce", &ptr);
cl = egetstr ("cl", &ptr);
UP = egetstr ("up", &ptr);
if (!tgetflag ("bs"))
BC = egetstr ("bc", &ptr);
tloaded = 1;
}
/* like tgetstr() but discard curses delay codes, for now anyways */
static char *
egetstr (name, sptr)
char *name;
char **sptr;
{
extern char *tgetstr();
register char c, *s;
s = tgetstr (name, sptr);
while ((c = *s) >= '0' && c <= '9')
s += 1;
return (s);
}
static
setuptty()
{
extern ospeed;
struct sgttyb sg;
ioctl (0, TIOCGETP, &orig_sgtty);
sg = orig_sgtty;
ospeed = sg.sg_ospeed;
sg.sg_flags &= ~ECHO; /* do our own echoing */
sg.sg_flags &= ~CRMOD; /* leave CR and LF unchanged */
sg.sg_flags |= XTABS; /* no tabs with termcap */
sg.sg_flags |= CBREAK; /* wake up on each char but can still kill */
ioctl (0, TIOCSETP, &sg);
ttysetup = 1;
}
#endif
#ifdef LATTICE_C
#include <dos.h>
/* (ANSI: ESC [ r ; c F) (r/c are numbers given in ASCII digits)
*/
c_pos (r, c)
int r, c;
{
printf ("%c[%d;%dF", ESC, r, c);
}
/* erase entire screen. (ANSI: ESC [ 2 j) */
c_erase()
{
printf ("%c[2j", ESC);
}
/* erase to end of line. (ANSI: ESC [ K) */
c_eol()
{
printf ("%c[K", ESC);
}
/* return 0 if there is a char that may be read without blocking, else -1 */
chk_char()
{
return (kbhit() == 0 ? -1 : 0);
}
/* read the next char, blocking if necessary, and return it. don't echo.
* map the arrow keys if we can too into hjkl
*/
read_char()
{
int c;
c = getch();
if (c == 0) {
/* get scan code; convert to direction hjkl if possible */
c = getch();
switch (c) {
case 0x4b: c = 'h'; break;
case 0x50: c = 'j'; break;
case 0x48: c = 'k'; break;
case 0x4d: c = 'l'; break;
}
}
return (c);
}
/* do whatever might be necessary to get the screen and/or tty back into shape.
*/
byetty()
{
}
#endif
/* read up to max chars into buf, with cannonization.
* add trailing '\0' (buf is really max+1 chars long).
* return count of chars read (not counting '\0').
* assume cursor is already positioned as desired.
* ESC: return -1 immediately.
*/
read_line (buf, max)
char buf[];
int max;
{
static char erase[] = "\b \b";
int n, c;
int done;
#ifdef UNIX
if (!ttysetup) setuptty();
#endif
for (done = 0, n = 0; !done; )
switch (c = read_char()) { /* does not echo */
case 0177: case CNTRL('h'): /* char erase */
if (n > 0) {
fputs (erase, stdout);
n -= 1;
}
break;
case CNTRL('u'): /* line erase */
while (n > 0) {
fputs (erase, stdout);
n -= 1;
}
break;
case '\r': /* EOL */
done++;
break;
case ESC: /* ESC to abort */
return (-1);
default: /* echo and store */
if (n >= max)
putchar (CNTRL('g'));
else {
putchar (c);
buf[n++] = c;
}
}
buf[n] = '\0';
return (n);
}
xXx