AstroCalc - Solar System Ephemeris for the hp 39g+

Introduction

As an amateur astronomer it is often useful to know the approximate positions of the Sun, Moon and planets to decide whether it will be worth going outside to observe them. Knowing their altitudes and azimuths can be particularly helpful in determining whether they will be accessible from your observing site.

There are several programs available for the PC and other platforms which will calculate data for solar system objects and produce detailed sky displays, complete with stars, constellation boundaries, etc. Though these programs can be very useful, it takes time to wait for the PC to boot up, find the CD-ROM, start the program, sit through the introduction animation, choose a planet, view its coordinates... and after all that Red Shift 3 which I use does not seem to have any way of reporting a planet's altitude in degrees so I cannot tell if it will clear the rooftops across the road.

On the other hand a fairly small program running on a programmable calculator can find all the relevant data to a sufficient accuracy and is quicker to use.


AstroCalc on the hp 39g+ Graphing Calculator

AstroCalc is a set of routines that run on a Hewlett-Packard 39g+ calculator and will calculate the following properties for the main solar system objects (i.e. Sun, Moon, planets Mercury to Neptune, plus dwarf planets Ceres, Pluto and Eris.):

  • The right ascension and declination. These are the sky equivalents of longitude and latitude and can be looked up on a star atlas to locate a planet with regard to nearby stars.
  • The altitude and azimuth, i.e. how high the planet is above the horizon in degrees and its bearing, such that North = 0°, East = 90° etc. Gives you an idea of which part of the sky to look in to find the planet.
  • The time the planet rises and sets and its azimuth at rising and setting. Tells you the period during which it is potentially visible.
  • The visual magnitude, i.e. the brightness of the planet on the same scale as commonly used for stars whereby the faintest stars visible to the unaided eye are about magnitude +6 and the brightest are magnitude 0 to -1. Knowing the magnitude of a planet makes it easier to find it visually.
  • The angular diameter of the planet in seconds of arc. The closer planets (especially Mars) vary greatly in their distance from the Earth and hence in apparent size. If the angular diameter is much less than ten arc seconds you are unlikely to see significant surface detail on a planet, even through a large telescope.
  • The phase, i.e. the percentage of the side of the planet facing Earth which is illuminated. Obviously varies from 0 to 100% for the Moon, and also for the inner planets Mercury and Venus. Mars can show a noticeable phase but the further planets are always nearly at 100% phase.
  • The distance of the planet from the Earth in astronomical units.
  • The time when the planet culminates, i.e. is due south from the the observer's location (in the northern hemisphere), and its altitude at this time. Since at culmination the planet is at its highest in the sky there is least atmospheric disturbance so this is generally the best time to view it.
  • The heliocentric longitude. This is the direction of the planet, measured from the First Point of Aries, as seen from the Sun, rather than as seen from the Earth. This information was added specifically for setting up a mechanical orrery.


Loading AstroCalc

Copy the files HP39DIR.000, HP39DIR.CUR, and all those with names starting with AST... and ending in /000 to an empty directory on a PC, start the HP39G Connectivity program and point it to the directory holding the files.
Program AstroCalc consists of a main file, called AstroCalc, and sixteen subprograms called .AstCulm, .AstEast, .AstDisp, .AstDOW, .AstEc2E, .AstEq2A, .AstJDN, .AstJul, .AstMoRS, .AstPara, .AstPlan, .AstRSet, .AstSRS, .AstSun, .AstTrAn, .AstU2LS, which should all be downloaded from the computer into the PROGRAM catalog part of the hp 39g+'s memory using the RECV option on the calculator.
(Also included are subroutines .AstASep, .AstA2Eq, .AstPreC and .AstPreM but these are not used in the current version of AstroCalc so there is no need to transfer them to the calculator.)
Don't forget that you can ''check' multiple programs in the calculator's RECV screen and download them in one go. This is a lot quicker than transferring them one at a time but is only briefly mentioned in the hp39g+ manual.)


Running AstroCalc

Before running AstroCalc for the first time it is recommended that you edit the settings in the program listing to customise it for your own location. Choose AstroCalc in the program catalog then press softkey 1 to EDIT.
Near the start of the program four values are stored into List 0 (L0) elements 1 to 4. These represent:

  1.  L0(1) holds the geographic longitude of the observer's location in decimal degrees. Longitudes East of the Greenwich meridian are positive and those West of Greenwich should be entered as negative.
  2.  L0(2) represents the geographic latitude of the observing location in decimal degrees. Northern latitudes are positive and Southern latitudes are negative.
  3.  L0(3) is the observer's altitude above sea level in metres.
  4.  L0(4) is the time zone offset relative to Greenwich Mean Time (now more officially called Universal Time). If your local time is earlier than Greenwich Mean Time enter the number of hours as a negative number, while if your time is later than GMT (i.e if you are East of Greenwich) then enter a positive number of hours here.

Observing site location settings
Section of program to edit to suit your own observing location.

Note on time zones: Roughly speaking, for every 15 degrees of longitude East of the Greenwich meridian the time zone increases by one hour, and for every 15 degrees West of Greenwich it decreases by one hour. For example the east coast of the US is about 75 degrees west of Greenwich and the time zone is -5 so this is the number to store in L0(4). However not all countries keep to the rule of 15 degrees = ± 1 hour so you need to consult an atlas to be certain. India is especially unusual in that its time is GMT +5 hours 30 minutes.
If 'Summer Time' or 'Daylight Saving Time', whereby the clocks are moved forwards an hour during Summer, is currently in effect then you can add an extra hour to the value stored to L0(4) to allow for this, so that for Daylight Saving Time on the US east coast you would store -5 + 1 = -4 into L0(4). However if you do this remember that the extra hour will be added to all times displayed in the program, even if you change the date to one when Daylight Saving Time would not be in force.

Also, when using the 'current time' as set in the calculator's clock, it is assumed that this is the local time, with the same offset from GMT as stored in L0(4). Thus, if you have changed the calculator's built-in clock to Daylight Saving Time, (by storing to the variable TIME) you also need to add 1 to the value stored to L0(4), otherwise calculated positions for the current time will be wrong.

The longitude and latitude need to be reasonably accurate (say to within a tenth of a degree) or the sunrise and sunset times especially will be in error. An atlas should provide sufficient precision. You may need to use the calculator's Hours.MinutesSeconds to Decimal Degrees conversion before entering the values into the program.

Once the program has been configured, run AstroCalc from the PROGRAM catalog. Once it has compiled a choose box appears on the screen with seven options (scroll down to get the last three) as shown in figure 1.

Main choose menu

Second page of menu
Figure 1 - Main AstroCalc menu, both pages. Calculate and display option is selected as default.

The menu options have the following meanings:

  • Use Sys Date - When AstroCalc starts up it reads the current date and time from the calculator's internal clock and uses this as the date for ephemeris calculations.
    Selecting this option causes the program to re-read the system date and time, used for instance if you have previously entered a new time. Check that your calculator's date and time are set correctly; mine seems to jump an hour or a day occasionally.
  • New Date/Time - Prompts for you to enter a new date and time to use for the calculations. A standard input box is displayed for each of the year, month, date and time in turn, as shown in figure 2.
    Prompt for entry of year
    Figure 2 - Input screen for entering the year.
    The current values are shown as the defaults in each case. Months are entered as a number from 1 to 12. The time should be entered as local time at the observing location, taking into account the time zone as explained above. Format for entry of time is HH.MM, i.e. hours in 24-hour-clock form, then a decimal point, then the number of minutes.
  • Calculate/Disp - Performs all the calculations for the currently set date and time, then displays the results. While the calculations are being carried out (which takes about a minute) a list is shown to indicate what stage the calculations have reached, figure 3.
    Progress display
    Figure 3 - Calculation progress indicator.

    The display will then change to the ephemeris table with the calculated results. It starts as figure 4:
    Start of ephemeris table
    Figure 4 - Initial results display.
    The current day of the week, year, month, date, time, start and end of twilight, observer's longitude and latitude, and time zone are shown along the top line. If twilight lasts all night then this is stated instead of the time.
    On the second line are the abbreviated column headings:
    • Ri Asc = Right ascension in hours and minutes.
    • Declin = Declination north (+) or south (-) of the celestial equator in degrees and minutes.
    • Alti = Altitude above (+) or below (-) the observer's horizon in degrees and minutes.
    • Azim = Azimuth (compass bearing) where N is 0°, E is 90° etc.
    • Rises = Local time at which the object rises above the horizon on the current date, in hours and minutes. If it is circumpolar or never rises on this day at the observer's latitude then this is stated instead.
    • Sets = Local time at which the Sun/Moon/planet sets.
    • Rise Az = Azimuth at time of rising, in degrees and minutes, or blank if does not rise.
    • Set Az = Azimuth at setting.
    • Mag = Visual magnitude (brightness), lower means brighter.
    • Ang Dia = Angular diameter of the object in arc seconds.
    • Phase = Percentage of the visible face of the object illuminated. For the Moon a positive phase means it is waxing and a negative phase means it is waning.
    • Dist AU = Distance of the planet from the Earth in astronomical units. For the Moon the fraction of the Moon's mean distance of 384401km is given.
    • Culmin = Local time at which the object culminates or transits on the current day, in hours and minutes.
    • Max Alt = The maximum altitude above the horizon reached by the object. Will be negative if it does not rise.
    • Helio L = The Heliocentric Longitude of the planet as seen from the Sun, measured anti-clockwise from the first point of Aries. Note that since the heliocentric longitude of the Sun itself would be meaningless, the figure displayed in the Sun's row is actually the longitude of the Earth. Also, the figure given for the Moon is its longitude centred on the Earth. For the purpose of setting an orrery, set the Moon's direction relative to the Earth to the figure given, imagining that the heliocentric longitude is centred on Earth.

    Down the left hand side are the names of the Sun, Moon or planets. These remain in place as the table is scrolled.
    The right and left cursor keys scroll along the table, two columns at a time. See figures 5 to 7.
    Ephemeris table
    Figure 5 - After scrolling right twice.

    Ephemeris table
    Figure 6 - Scrolled to show observer's location along top line.

    End of ephemeris table
    Figure 7 - Scrolled to end of table.
    After you reach the right-hand end of the table ('Helio L'), the next scroll to the right will take you back to the start of the table.
    The up and down cursor keys move the results display up or down by five 'planets' at a time, i.e showing the Sun to Mars, Ceres to Neptune or Pluto and Eris as in figure 8.
    Ephemeris table p2

    Ephemeris table p3
    Figure 8 - After scrolling down once then twice to outer planets.
    To leave the results display and return to the main menu, press ENTER or HOME.

  • Exit program - Ends the program and returns to the program catalog. Note that you need to use this option to exit the program. Pressing 'ON' does not break into the program.
  • Location... - Prompts for a new observing location in terms of longitude, latitude and time zone, as in figure 9.
    New longitude prompt
    Figure 9 - Input screen for new longitude.
    Enter the longitude and latitude in degrees and minutes format, separated by a decimal point.
    This new location will then be used for future calculations until the program is exited.
  • Easter Sunday - Calculates the date of Easter Sunday in the currently set year and displays it as in figure 10.
    Easter Sunday display
    Figure 10 - Message box for Easter Sunday.
    Press OK or ENTER to return to the main menu.
  • Julian Day - Calculates the Julian Day Number of the currently set date and time and displays the results as in figure 11.
    Julian day display
    Figure 11 - Message box for Julian Day Number.
    As a reminder the current date is shown in YYYY-MM-DD format with the local time below.
    Press OK or ENTER to return to the main menu.

Hardware Requirements

AstroCalc runs on a Hewlett-Packard 39g+ calculator. It should also work on the 39g and 40g but will take about three times as long to calculate positions, due to the slower processor. It is written in HP-BASIC and may work on the 48g, 49g and 50g but has not been tried.

The program takes up around 34 kilobytes of RAM when run. It also uses 1.4 kilobytes to store matrix M9 and list L0 This means the calculator must have about 36 kilobytes of free RAM to run AstroCalc. This, plus use of the DISPXY command, rules it out for the 38g model.


Accuracy

In general calculated positions, i.e. right ascensions, declinations, altitudes and azimuths, should be correct to within a few minutes of arc for dates close to the epoch of the orbital data (year 2000 – see later). Accuracy decreases the further into the past or future we try to calculate positions, since many of the sources of error have a cumulative effect over time. E.g. if the figure for the orbital period of the Earth is in error by 1 part in 100 000, then after 1000 years the error will be 1 part in 100 or over 3 degrees of longitude. Also no allowance is made for precession of the equinoxes, nor for long-term changes in planetary orbits. The Moon is likely to show the greatest discrepancy since its orbit is subject to complicated perturbations. Generally the program is not sufficiently accurate more than a few centuries either side of year 2000.

However targets should be within the low magnification field of view of a telescope pointed at the reported coordinates.
Rising and setting times should also be within a few minutes of the correct times, though in practice the local topography (hills on the horizon) is likely to be a bigger effect than any errors in the calculation.

There are several minor corrections which are not made by AstroCalc, in order to keep the size of the program down and the speed up. These include:

  • Precession as already mentioned. Precession is caused by the gravity of the Sun and Moon acting on the non-spherical Earth. The direction of the Earth's rotation axis revolves in a circle with a radius of about 23.5° and a period of 25800 years. This causes the position of the zero line of right ascension to shift.
  • There is a small wobble superimposed on the precession, known as nutation, caused by the slightly varying distances of the Sun and Moon. It causes the direction of the rotation axis to vary by about 9 seconds of arc over a period of 18.6 years. There is also the Chandler wobble which is a movement of the Earth relative to the rotation axis, rather than the rotation axis relative to the stars. The Chandler wobble has an amplitude of 0.3 seconds of arc as the geographic location of the poles shifts by about 10 metres over a short period.
  • Aberration is due to the speed of light not being infinite. The Earth's orbital speed is a tiny but not zero fraction of the speed of light and has the effect of a vector addition to the direction of light from a distant planet. Unless the Earth happens to be travelling directly towards or away from a planet, aberration causes the apparent direction of the light from the planet to change slightly, by a maximum angle of 20.5 arcseconds. A second effect is that when we observe Jupiter, for example, the light we are seeing left the planet about 40 minutes ago so we are not seeing Jupiter where it is now but where it was 40 minutes earlier. Again the error is of the order of a few arcseconds.
  • Geocentric parallax. The calculated positions of the planets are really only correct if you were observing from the centre of the Earth. From a position on the surface you get a very slightly different view depending on the angle between you, the planet and the centre of the Earth. For Jupiter the difference is only about 2 arcseconds. For the much closer Moon parallax can make a difference of as much as a degree of arc and AstroCalc does make a suitable correction for the Moon.
  • The calculations assume that each planet purely orbits the Sun and do not take account of the gravitational influence of the other planets in the solar system. This can result in a maximum error of up to a few arcminutes in the calculated positions, but usually less.
  • Planetary orbits are not fixed but vary slowly, typically over timescales of thousands of years. The direction of perihelion gradually rotates around the Sun, and the degree of eccentricity and orbital inclination also change. (It has even been shown that over the longest timescales, hundreds of millions of years, orbits are mathematically 'chaotic' and is fundamentally impossible to predict them from their current parameters.)


Variables Used

The simple programming model of HP BASIC limits variable names to a single character and does not allow local variables, so care is needed when using variables in subroutines. AstroCalc uses all the HOME variables A-Z plus some others, and will therefore overwrite any existing information stored in them. Some variable names are used as a general-purpose temporary scratchpad but others have a fairly consistent purpose. The variables changed by the program are:

  • The plotting limit variables Xmin, Xmax, Ymin, Ymax
  • The Angle unit is set to degrees.
  • The number Format is set to floating decimal, and the number of Digits in fixed-decimal mode is set to 2.
  • List L0 holds details of the observer's location.
  • Expressions are stored into functions F8(X), F9(X) and F0(X) to perform common modulo operations (e.g. putting an angle into the range 0-360).
  • Real variables A to Z and q hold:
    A = altitude or azimuth at rising
    B = usually ecliptic latitude
    C = usually ecliptic longitude
    D = date in month
    E = declination
    H = usually rise time
    J = Julian day number
    L = local sidereal time
    M = month
    P = sometimes planet number
    R = right ascension
    S = usually setting time
    T = current (universal) time
    U = calculated universal time
    W = day of week
    Y = current year
    Z = azimuth, often at setting
    Other variables hold temporary values.
  • Matrix M9 stores the calculated results prior to display.


Known Bugs

There are some unwanted pixels set on the first line of the table display, near the first letter of the day of the week. Seems to be a mis-feature of the 39g+ due to wrap-around of letters off the screen.
Setting the latitude to exactly +90 or -90 causes an error, because tan(±90) is infinite. 89.9° is OK.
When calculating altitude at culmination no account is taken of atmospheric refraction so maximum altitudes within a few degrees of the horizon will be slightly wrong.


Orbital Elements

The orbits of the planets are calculated from their 'orbital elements' for the date 2000 January 1 12:00 UT ('January 1.5'), julian day number 2451545, known as the 'epoch' of the elements. The orbital elements, with their symbols, are as follows:

  • Tp — The orbital period of the planet, measured in Earth 'tropical years' of 365.242191 solar days. The orbital period, also known as the sidereal period, is the time taken for the planet to make one complete circuit of its orbit relative to the fixed stars. The tropical year of a planet is the time between successive passages of the planet through a fixed point in its orbit, such as its perihelion or ascending node. For the Earth a tropical year is the time between successive passages of the Sun through the first point of Aries, or equivalently the time between equinoxes, and is therefore the year that a calendar must reproduce if seasons are not to become out of step with dates. Note that the sidereal period of the Earth's orbit (365.2564 days) is slightly longer than one tropical year because of precession — it differs by 1 part in 25800. (There is also orbital precession where the direction of perihelion slowly rotates around the Sun, and this affects all the planets.) Thus Tp for the Earth is not exactly 1.00000. Strictly speaking the orbital period need not be separately tabulated since it can be calculated from the semi-major axis. Tp = a1.5
  • e (greek epsilon) — The longitude at epoch. Often the symbol L0 is used instead. The longitude at epoch is the longitude that a fictitious body moving in a circular orbit at constant angular speed but with the same period as the planet's true elliptical orbit would have at the epoch date. All longitudes of the orbital elements are heliocentric (sun-centred). Note that e is measured from the direction of perihelion.
  • W (greek capital omega) — The longitude of the ascending node. This is the direction of the point where the plane of the orbit of the planet crosses the Earth's orbital plane (the plane of the ecliptic), measured from the first point of Aries=0.
  • v (greek lowercase omega with a bar over) — The longitude of the perihelion. This is the angle of the point where the planet is closest to the Sun, measured from the first point of Aries. Note that some references use w (omega without a bar) to mean the angle of perihelion measured from the ascending node (the 'argument of the perihelion'.) v = w+W
  • e — The eccentricity of the orbit. A measure of how much the orbit differs from a circle. e=0 means a perfectly circular orbit, larger values up to a maximum of e=1 mean increasingly elliptical orbits.
  • a — The semi-major axis of the orbit. This is half the length of the orbit along its long axis and is the same as the radius of a circular orbit. a is normally measured in astronomical units, i.e. multiples of the semi-major axis of the Earth's orbit.
  • i — The inclination of the orbit. The angle the planet's orbit is tilted at relative to the plane of the Earth's orbit.
  • M0 — The 'mean anomaly' at the epoch. This is the position of a fictitious body following a circular orbit with the same period as the planet, at the date of the epoch, measured from the ascending node. M0 is sometimes used to define the position of the planet at the epoch instead of e, see above.
  • q0 (theta zero) — The angular diameter the planet would have, in seconds of arc, if it were at a distance of 1 astronomical unit (AU) from the Earth. Not a true orbital element but used to calculate the apparent sizes of the planets.
  • V0 — The visual magnitude the planet would have if it were 1 AU from the Sun and 1 AU from the Earth, and with a phase of 100%. Again not really an orbital element but used to calculate the apparent brightness of the planet. Note that the actual brightness falls off with the square of the planet's distance from the Sun and the square of its distance from the Earth.

The table below lists the orbital elements of the eight planets in the solar system, plus three dwarf planets (Ceres, Pluto, Eris) for the epoch 2000 January 1.5. Surprisingly, different sources quote slightly different values for these elements, such as an orbital period for Pluto ranging from 246.7 to 248.0 years, and some of the data for Eris will still have a fairly large uncertainty since this dwarf planet was only discovered a few years ago and has only completed a very small part of its orbit.
The figures below are what seem to be the most reliable.
PlanetTp
Period
(tropical years)		e
Long. at epoch
(degrees)		v or L0
Long. of peri-helion
(degrees)		M0
Mean anomaly
(degrees)		a
Semi-major axis of orbit
(AU)		e
Eccent. of orbit
(no units)		i
Inclin. of orbit
(degrees)		w
Arg. of peri-helion
(degrees)		W
Long. of the ascend. node
(degrees)		q0
Angular diam. at 1 AU
(arcseconds)		V0
Visual mag. at 1 AU
(mags.)
Mercury0.240852252.25177.456174.7950.387100.205637.00529.12548.3316.74-0.42
Venus0.615211181.980131.53350.4160.723330.006773.39554.85276.68116.92-4.40
Earth1.000038100.464102.947357.5291.000000.016710.000114.208348.73917.83-3.86
Mars1.880889355.453336.04119.3731.523660.093411.851286.46249.5799.36-1.52
Ceres4.599840161.741153.8007.9412.765800.0780010.60773.10080.7001.304.30
Jupiter11.86223634.40414.75420.0205.202600.048391.305274.290100.464196.74-9.40
Saturn29.45776949.94492.432317.0219.537070.054152.484338.717113.715165.60-8.88
Uranus84.013843313.232170.964141.05019.191260.047120.77096.73474.23065.80-7.19
Neptune164.792460304.88044.971256.22530.06896 0.008591.769273.249131.72262.20-6.87
Pluto247.687000238.929224.06714.88239.481690.2488117.142113.764110.3033.20-1.00
Eris557.20000021.649187.301194.34867.670000.4417744.187151.43135.8703.30-1.10


Component Programs

The individual subprograms which make up AstroCalc are listed below in a readable format, with comments on parameters passed to and from them and their purpose.
Programs are adapted from 'Practical Astronomy with your Calculator', 3rd edition, by Peter Duffet-Smith, published by Cambridge University Press in 1988. Planetary orbital elements from been updated from epoch 1990 as used in the book to epoch 2000.
Note: the character Þ represents the 'STO' arrow.
True line breaks in the program are shown by ¿, other line breaks are just to fit the listing into the table.


Program ListingsExplanation

Program AstroCalc


ERASE:DISP 1;" AstroCalc":¿
{0,0,0,0,0,0}ÞL0:¿
MAKEMAT(0,12,15)ÞM9:¿
-2.298ÞL0(1):¿
53.596ÞL0(2):¿
100ÞL0(3):¿
1ÞL0(4):¿
2000ÞL0(5):¿
QUOTE((X+L0(4)) MOD 24)ÞF0(X):
QUOTE(X MOD 360)ÞF9(X):
QUOTE(X MOD 24)ÞF8(X): ¿
1ÞAngle:2ÞDigits:¿
L0(5)ÞY:1ÞM:0ÞD:0ÞT:¿
RUN ".AstJDN":¿
(J-2451545)/36525ÞT:¿
23.439292-(46.815*T+0.0006*T²-0.00181*T*T²)/3600ÞL0(6):¿
1ÞO:¿
DO¿
IF O==1 THEN DISP 3;"Using current date":WAIT 1: ¿
HMS®(TIME)-L0(4)ÞT:¿
ROUND(FRAC(DATE*100)*10000,0)ÞY:¿
INT(DATE)ÞM:¿
INT(FRAC(DATE)*100)ÞD:¿
END:¿
ERASE:3ÞO:¿
CHOOSE O;"AstroCalc";"Use Sys Date";"
New Date/Time";"Calculate/Disp";"Exit program";
"Location...";"Easter Sunday";"Julian Day":¿
IF O==2 THEN¿
INPUT Y;"New Date";"Year";"Format '2006'";Y:¿
INPUT M;"New Date";"Month";"Format '5'";M:¿
INPUT D;"New Date";"Date";"Format '21'";D:¿
2ÞFormat:INPUT T;"New Time";"Local Time";
"Format 'HH.MM'";
ROUND(®HMS(F0(T)),2):
1ÞFormat:HMS®(T)-L0(4)ÞT:¿
DISP 3;"Using new time":WAIT 1: ¿
END: ¿
IF O==5 THEN 2ÞFormat: ¿
INPUT L;"Location";"Longitude";
"Format 'Degs.Mins' W -ve, E +ve";
ROUND(®HMS(L0(1)),2):
HMS®(L)ÞL0(1):¿
INPUT L;"Location";"Latitude";
"Format 'Degs.Mins' N +ve, S -ve";ROUND(®HMS(L0(2)),2):
HMS®(L)ÞL0(2):¿
1ÞFormat:¿
INPUT L;"Location";"Time Zone";
"Hours ahead + or behind - GMT";L0(4):
HMS®(L)ÞL0(4): ¿
DISP 3;"New location stored":WAIT 1:¿
END: ¿
IF O==6 THEN RUN ".AstEast":¿
END:¿
IF O==7 THEN RUN ".AstJul":¿
END:¿
IF O==3 THEN¿
ERASE:DISP 1;"Calculating":¿
DISP 3;" Planets":RUN ".AstPlan":¿
DISP 4;" Moon":RUN ".AstMoon":RÞM9(2,1):
EÞM9(2,2):SÞM9(2,9):NÞM9(2,10):PÞM9(2,11):
ZÞM9(2,12):R*15ÞM9(2,15): ¿
RUN ".AstMoRS":¿
DISP 5;" Sun":RUN ".AstSun":RÞM9(1,1):
EÞM9(1,2):-26.74+LOG(K²)/LOG(2.512)ÞM9(1,9):
NÞM9(1,10):1ÞM9(1,11):KÞM9(1,12):¿
RUN ".AstSRS":¿
DISP 6;" Alt/Az":RUN ".AstU2LS":¿
FOR P=1 TO 12;¿
M9(P,1)ÞR:M9(P,2)ÞE:¿
RUN ".AstEq2A":AÞM9(P,3):ZÞM9(P,4):¿
END:¿
RUN ".AstCulm":¿
RUN ".AstDOW":¿
1ÞC:1ÞR:¿
BEEP 1500;0.3:BEEP 1300;0.5:¿
DO RUN ".AstDisp":¿
GETKEY O:¿
R-5*(O==25.1)+5*(O==35.1)ÞR:¿
IF R<1 THEN 11ÞR:END:IF R>11 THEN 1ÞR:END: ¿
IF O==34.1 OR O==36.1 THEN
(C+2*(O-35.1)) MOD 14ÞC:END: ¿
UNTIL O==105.1 OR O==31.1 END:¿
3ÞO:END:¿
UNTIL O==4 END: 		This is the main routine and the one which is run directly to start AstroCalc. It calls the other routines as needed.
Matrix M9 is used to store the calculated results ready for display.
List L0 holds various constants used by the program. These are:
  • L0(1) is the geographic longitude of the default observing site in decimal degrees. Longitudes East of the Greenwich Meridian are taken to be positive and longitudes West of Greenwich are negative.
  • L0(2) is the geographic latitude of the default observing site in decimal degrees. Latitudes North of the equator are positive and latitudes South of the equator are negative.
  • L0(3) is the altitude in metres above sea level of the default observing site. This is only used in correcting for the parallax of the Moon and in practice makes virtually no difference unless the altitude is very large.
  • L0(4) is the time zone offset in hours relative to Greenwich Mean time. This is the number of hours that local civil time is ahead of Greenwich mean time. Places West of Greenwich have a negative offset and places East have a positive offset. It also includes any extra hour added for 'Daylight Saving Time'.
  • L0(5) is the epoch year used for calculations. It would be possible to convert calculated Right Ascensions and Declinations to a new epoch stored in L0(5) to correct for precession, but at present this is not done.
  • L0(6) stores the obliquity of the ecliptic for the epoch year in L0(5). This is not entered directly but is calculated from the epoch.
The current date and time stored in the calculator's built-in clock is used as the initial date for calculations.
Note that all times input by and displayed to the user are local civil times but internally the program uses Universal Time (essentially the same as Greenwich Mean Time), converting to and from civil time as needed.
A menu is then displayed which allows a new date and time or new location to be entered, the current date to be re-read or the positions of the Sun, Moon and Planets to be calculated and displayed in tabular form.
Once the results are being displayed they can be scrolled around using the cursor keys, or pressing ENTER or HOME returns to the menu.
Selecting 'Exit Program' from the menu returns to the calculator's program screen or home screen, depending on where AstroCalc was called from.

Three functions are also defined, to be used in the subroutines. These are:
  • F0 adds the civil time offset to a time and adjusts it to the range 0 to 24.
  • F9 converts an angle into the range 0 to 360 degrees.
  • F8 converts a time into the range 0 to 24 hours.

Note that if further menu options, with numbers above 7, are added to the main menu, the code to deal with them needs to be placed before that for option 3 (...O==3...), and the value of variable O (the selected option number) should not be changed by them.

Program .AstDOW


RUN ".AstJDN":¿
INT(J-(F0(T))/24+2.5) MOD 7ÞW:¿
 		Calculates the day of the week of the current date.
On entry Y=year, M=month, D=date in month, T=local time.
Returns W=day of week where 0=Sunday, 1=Monday etc.

Program .AstJDN


Y+(Y<0)ÞQ:MÞJ:¿
IF M<2.1 THEN Q-1ÞQ: J+12ÞJ:END:¿
(2-INT(Q/100)+INT(Q/400))*((Y*10000+M*100+D) ³15821015)+INT(365.25*Q-0.75*(Q<0))+ INT(30.6001*(J+1))+D+ 1720994.5+T/24ÞJ:¿
 		Calculates the Julian Day Number (JDN)of the current date and time.
JDN is defined as the number of days since Greenwich noon on January 1st 4713 B.C. This year was chosen because it predates all known astronomical observations so JDN is always positive.
The usefulness of Julian Day Numbers is that they make it easy to calculate the number of days between two dates, just by the difference between their JDNs.
This routine takes account of the change to the Gregorian calendar on 15th October 1582 and thus works for any date from 1 AD onwards.
On entry Y=year, M=month, D=date in month, T=Universal time.
Returns J=Julian Day Number.

Program .AstU2LS


TÞX:0ÞT:¿
RUN ".AstJDN":¿
(J-2451545)/36525ÞT:¿
F8(6.697374558+2400.051336*T+ 0.000025862*T²)ÞT:¿
F8(X*1.002737909+T)ÞT:¿
F8(T+L0(1)/15)ÞL:¿
XÞT:¿
 		Converts Universal Time to Local Sidereal Time.
A sidereal day is the time between successive transits of the line of the vernal equinox, approximately equal to the time the Earth takes to rotate once on its axis, 23 hours and 56 minutes and 4 seconds. (Because of precession a sidereal day is actually 0.0084 seconds shorter than the rotation period.)
Sidereal time is defined relative to the stars rather than the Sun such that there are 24 hours of sidereal time in one sidereal day. Local sidereal time is the same as the right ascension of stars currently on the observer's meridian, and thus depends on longitude as well as the time.
On entry Y=year, M=month, D=date in month, T=Universal time.
Returns L=local sidereal time.

Program .AstLS2U


F8(L-L0(1)/15)ÞU:¿
TÞX:0ÞT:¿
RUN ".AstJDN":¿
(J-2451545)/36525ÞT:¿
F8(6.697374558+2400.051336*T+ 0.000025862*T²)ÞT:¿
F8(U-T)*0.9972695663ÞU:¿
XÞT:¿
¿
 		Converts Universal Time to Local Sidereal Time, i.e. the reverse of the above routine.
On entry Y=year, M=month, D=date in month, L=local sidereal time.
Returns U=universal time.

Program .AstEq2A


F8(L-R)*15ÞX:¿
ASIN(SIN(E)*SIN(L0(2))+COS(E)*COS(L0(2))* COS(X))ÞA:¿
ACOS((SIN(E)-SIN(L0(2))*SIN(A))/ COS(L0(2))/COS(A))ÞZ:¿
IF SIN(X)>0 THEN 360-ZÞZ:END:¿
IF A>-0.75 THEN A+(159+19.6*A+0.02*A²)/ (1+0.505*A+0.0845*A²)/282ÞA:END:¿
 		Converts Equatorial coordinates (i.e. right ascension and declination) to Horizon coordinates (i.e. altitude and azimuth).
The altitude and azimuth of a planet are very useful to the observer in deciding whether or not it is currently visible, and depend on its right ascension and declination plus the current time and the observer's longitude and latitude.
The altitude is corrected for atmospheric refraction, which causes objects low down in the sky to appear higher than they really are, by a maximum of just over half a degree. (Strictly the correction depends on air pressure and temperature but since these are unlikely to be known typical values are used.)
On entry L=local sidereal time, R=right ascension (hours), E=declination (degrees).
Returns A=altitude above horizon (degrees), Z=azimuth (North=0, East=90, South=180, West=270)
Local sidereal time rather than the current date is used as the input to avoid the need to repeatedly calculate the same LST when calculating the positions of several objects.

Program .AstA2Eq


ASIN(SIN(A)*SIN(L0(2))+COS(A)*COS(L0(2))* COS(Z))ÞE:¿
ACOS((SIN(A)-SIN(L0(2))*SIN(E))/COS(L0(2))/ COS(E))ÞR:¿
IF SIN(Z)>0 THEN 360-RÞR: END:¿
F8(L-R/15)ÞR:¿
 		Converts Horizon coordinates to Equatorial coordinates (i.e. the reverse of the previous routine).
This routine is not currently used by AstroCalc but could be added as an extra function on the main menu.
On entry L=local sidereal time, A=altitude, Z=azimuth.
Returns R=right ascension (hours), E=declination (degrees).

Program .AstEc2E


ASIN(SIN(B)*COS(L0(6))+COS(B)*SIN(L0(6))* SIN(C))ÞE:¿
COS(C)ÞX:¿
SIN(C)*COS(L0(6))-TAN(B)*SIN(L0(6))ÞQ:¿
ARG((X,Q))ÞR:¿
F8(R/15)ÞR:¿
 		Converts Ecliptic coordinates to Equatorial coordinates.
Ecliptic coordinates use the 'First Point of Aries' (the direction of the Vernal Equinox) as the zero of longitude and the plane of the Earth's orbit around the Sun as the zero of latitude.
Equatorial coordinates also use the First Point of Aries as the zero of right ascension, and the zero of declination is the plane of the Earth's equator.
Conversion between the two depends on the current 'obliquity of the ecliptic', the angle of the Earth's axial tilt, which slowly changes over time due to perturbations by the Moon and planets and was calculated when AstroCalc started.
On entry B=ecliptic latitude (decimal degrees), C=ecliptic longitude, L0(6)=obliquity
Returns R=right ascension (hours), E=declination (degrees).

Program .AstASep


ACOS(SIN(E)*SIN(C)+COS(D)*COS(C)* COS((R-B)*15)ÞS:¿
 		Calculates the angular separation between two objects from their right ascensions and declinations.
Not currently used by AstroCalc but could be added as another menu option.
On entry R=right ascension of first object (decimal hours), E=declination of first object (decimal degrees), B=right ascension of second object, C=declination of second object.
Returns S=separation (decimal degrees).

Program .AstRSet


SIN(E)/COS(L0(2))ÞH:¿
IF ABS(H)>1 THEN IFTE(ABS(L0(2))> ABS(90*SIGN(L0(2))-E),99,-99)ÞH:¿
ELSE ACOS(SIN(L0(2))/ COS(E))ÞS:¿
ASIN(TAN(0.567)/TAN(S))ÞZ:¿
ACOS(H)-ZÞA: 360-AÞZ:¿
ACOS(-TAN(L0(2))*TAN(E))/15ÞH:¿
ASIN(SIN(0.567)/SIN(S))/15/ COS(E)Þq: ¿
F8(R+H+q)ÞL:¿
RUN ".AstLS2U":UÞS:¿
F8(24+R-H-q)ÞL:¿
RUN ".AstLS2U":UÞH:¿
END:¿
 		Calculates the times of rising and setting of an object on the current date from its right ascension and declination, and the observer's geographic location in L0(1) and L0(2). The times are corrected for atmospheric refraction which causes a planet to appear on the horizon while it is still really about 34 minutes of arc below the horizon.
The azimuth positions of rising and setting are also calculated, though these are not corrected for refraction and may be slightly out. On entry Y=year, M=month, D=date in month, R=right ascension (decimal hours), E=declination (decimal degrees)
Returns A=azimuth at rising (decimal degrees), Z=azimuth at setting, H=time of rising (local civil time, decimal hours), S=time of setting.
If the object is circumpolar then H=99. If it never rises then H=-99. In these cases A, Z and S are undefined.

Program .AstPreM


YÞG:MÞV: DÞX:TÞq:¿
OÞY:1ÞM: 0ÞD:0ÞT:¿
RUN ".AstJDN":¿
(J-2451545)/36525ÞT:¿
0.6406161*T+0.0000839*T²+ 0.000005*T*T²ÞY:¿
0.6406161*T+0.0003041*T²+ 0.0000051*T*T²ÞM:¿
0.556753*T-0.0001185*T²- 0.0000116*T*T²ÞD:¿
[[COS(Y)*COS(D)*COS(M)-SIN(Y)*SIN(M), COS(Y)*COS(D)*SIN(M)+SIN(Y)*COS(M), COS(Y)*SIN(D)], [-SIN(Y)*COS(D)*COS(M)-COS(Y)*SIN(M), -SIN(Y)*COS(D)*SIN(M)+COS(Y)*COS(M), -SIN(Y)*SIN(D)], [-SIN(D)*COS(M), -SIN(D)*SIN(M), COS(D)]]ÞM0:¿
L0(5)ÞY:1ÞM: 0ÞD:0ÞT:¿
RUN ".AstJDN":¿
(J-2451545)/36525ÞT:¿
0.6406161*T+0.0000839*T²+ 0.000005*T*T²ÞY:¿
0.6406161*T+0.0003041*T²+ 0.0000051*T*T²ÞM:¿
0.556753*T-0.0001185*T²- 0.0000116*T*T²ÞD:¿
[[COS(Y)*COS(D)*COS(M)-SIN(Y)*SIN(M), -SIN(Y)*COS(D)*COS(M)-COS(Y)*SIN(M), -SIN(D)*COS(M)], [COS(Y)*COS(D)*SIN(M)+SIN(Y)*COS(M), -SIN(Y)*COS(D)*SIN(M)+COS(Y)*COS(M), -SIN(D)*SIN(M)], [COS(Y)*SIN(D), -SIN(Y)*SIN(D), COS(D)]]*M0ÞM0:¿
GÞY:VÞM: XÞD:qÞT:¿
 		Due to the phenomenon of precession of the equinoxes (the Earth's rotation axis is not fixed in direction but 'wobbles' in a circle with a period of 25800 years) the direction of the vernal equinox, which defines the zero of right ascension, is not fixed but moves by 50.3 seconds of arc per year relative to the stars. Thus a star's right ascension (and declination) slowly changes over time. For this reason the positions shown in star atlases are only exactly correct at one particular date, known as the 'epoch'. Most current atlases are for epoch 2000.
The positions of the Sun, Moon and planets calculated by AstroCalc are also correct for epoch 2000 and therefore should agree with atlases. It is possible however to convert epoch 2000 coordinates to the corresponding coordinates in another epoch by means of matrix multiplication, which could be useful when determining the apparent position of a solar system object several decades or centuries before or after year 2000. This routine generates a suitable matrix.
Note that the current version of AstroCalc does not make use of these corrections since it would make the program significantly larger and slower for only a minor improvement in accuracy.
On entry O (the letter O)=the year of epoch 1, L0(5)=the year of epoch 2.
Returns M0= a 3x3 matrix to convert from epoch 1 to epoch 2.

Program .AstPreC


M0*[COS(R*15)*COS(E), SIN(R*15)*COS(E),SIN(E)]ÞM9:¿
ASIN(M9(3))ÞE:¿
ARG((M9(1),M9(2)))ÞR:¿
(R/15) MOD 24ÞR:¿
 		Once a matrix has been set up by routine .AstPreM, routine .AstPreC is used to actually perform the conversion of right ascension and declination from epoch 1 to epoch 2.
On entry R=right ascension (decimal hours) in epoch 1, E=declination (decimal degrees) in epoch 1.
Returns R=right ascension in epoch 2, E=declination in epoch 2.

Program .AstSun


RUN ".AstJDN":¿
F9(360*(J-2451545)/365.242191)ÞN:¿
N+280.464-282.947ÞN:¿
0.016713ÞK:¿
RUN ".AstTrAn":¿
F9(N+282.947)ÞC:¿
0ÞB:¿
RUN ".AstEc2E":¿
(1-0.016713²)/(1+0.016713*COS(N))ÞK:¿
1919/KÞN:¿
 		Calculates the right ascension and declination of the Sun on the current date and time. Also finds its distance from the Earth and its angular diameter.
On entry Y=year, M=month, D=date in month, T=universal time
Returns B=ecliptic latitude (decimal degrees), C=ecliptic longitude (decimal degrees), R=right ascension (decimal hours), E=declination (decimal degrees), K=distance from Earth (in Astronomical Units), N=angular diameter of sun (seconds of arc). One astronomical unit is (for practical purposes) the average radius of the Earth's orbit, 1AU=149597870 km.

Program .AstSRS


TÞO:¿
12-L0(1)/15ÞT:¿
RUN ".AstSun":¿
RUN ".AstRSet":¿
IF ABS(H)<90 THEN SÞW: HÞT:¿
RUN ".AstSun":EÞK:¿
RUN ".AstRSet":¿
HÞG:AÞV:¿
WÞT:¿
RUN ".AstSun":¿
RUN ".AstRSet":¿
VÞA:F0(G-N/108E3)ÞH: F0(S+N/108E3)ÞS:¿
END:¿
HÞM9(1,5):SÞM9(1,6): AÞM9(1,7):ZÞM9(1,8):¿
(COS(108)-SIN(L0(2))*SIN(E))/COS(L0(2))/COS(E)ÞK:¿
OÞT:¿
IF ABS(K)>1 THEN 99*-SIGN(K)ÞK: ELSE ACOS(K)/15ÞK:¿
RÞL:RUN ".AstLS2U":F0(U+K)ÞN: F0(U-K)ÞK:END:¿
 		Finds the times of sunrise and sunset on the current date, the azimuths of rising and setting, and the beginning and end of astronomical twilight.
Calculating the rising and setting times for the Sun is more complicated than for a star, for two reasons:
  1. The Sun's right ascension is not constant but changes through the day.
  2. The Sun is not a point but has a significant angular size and sunrise is when the top edge of the Sun first becomes visible, rather than its centre.

To find the rise time of the Sun we need to know its right ascension at the time it rises, which is what we are trying to find!
The solution adopted here is to calculate the Sun's position at local noon on the day in question, and calculate when it would rise and set if its position did not alter. (These times will be wrong by several minutes.) Next we recalculate its positions at these two times to obtain more accurate estimates of its positions at rising and setting, and finally recalculate the rising and setting times using the improved estimated positions.
The times are corrected by the routine .AstRSet for atmospheric refraction and further corrected for the Sun's radius of approx 0.267 degrees.

There are three stages of twilight:
  • During civil twilight the Sun is between 0 and 6 degrees below the horizon and it is still light enough to see one's way around outside without artificial light.
  • During nautical twilight the Sun sinks from 6 degrees to 12 degrees below the horizon and it is still possible to discern the horizon between sea and sky when at sea.
  • During astronomical twilight the Sun is between 12 and 18 degrees below the horizon and fainter stars gradually become visible. Once the Sun reaches 18 degrees below the horizon the sky does not become any darker and this is the end of astronomical twilight.

Thus astronomical twilight starts in the morning several hours before sunrise when the Sun reaches an altitude of -18 degrees moving upwards, and ends in the evening when the Sun drops below -18 degrees altitude.
On entry Y=year, M=month, D=date in month, L0(1)=observer's longitude, L0(2)=observer's latitude.
Returns H=time of sunrise (decimal hours, local civil time), S=time of sunset, A=azimuth at rising (decimal degrees), Z=azimuth at setting, K=time of start of morning twilight (decimal hours, local time), N=time of end of evening twilight.
The values of H, S, A and Z are also stored in matrix M9, row 1, columns 5 to 8.
Notes:For locations within the arctic circles there will be dates around midsummer and midwinter when the Sun does not set or does not rise. If the Sun is circumpolar then H>90 while if it never rises then H<-90.
For temperate latitudes there are periods in summer when the Sun is never less than 18 degrees below the horizon, which means the sky never becomes completely dark. This is indicated by K>90. For areas very close to the poles it is even possible in winter for the Sun to be always more than 18 degrees below the horizon, so there is never any twilight. In this situation K<-90.

Program .AstPlan


RUN ".AstJDN":¿
J-2451545ÞI:¿
F9(360/365.242191*I/1.00004+100.466-102.937)ÞN:¿
0.01671ÞK:¿
RUN ".AstTrAn":¿
F9(N+102.937)ÞF:FÞM9(1,15):¿
(1-K²)/(1+K*COS(N))ÞV:¿
RUN ".AstU2LS": ¿
FOR P=1 TO 10;¿
CASE¿
IF P==1 THEN {0.240852,252.251,77.456,0.20563,0.3871,7.005,48.331,6.74,-0.42}ÞL9:END¿
IF P==2 THEN {0.615211,181.98,131.533,0.00677,0.72333,3.395,76.681,16.92,-3.8}ÞL9:END¿
IF P==3 THEN {1.880889,355.453,336.041,0.09341,1.52366,1.851,49.579,9.36,-1.52}ÞL9:END¿
IF P==4 THEN {4.59984,161.741,153.8,0.078,2.7658,10.607,80.7,1.3,4.3}ÞL9:END¿
IF P==5 THEN {11.862236,34.404,14.754,0.04839,5.2026,1.305,100.464,196.74,-9.4}ÞL9:END¿
IF P==6 THEN {29.457769,49.944,92.432,0.05415,9.53707,2.484,113.715,165.6,-8.88}ÞL9:END¿
IF P==7 THEN {84.013843,313232,170.964,0.04712,19.19126,0.77,74.23,65.8,-7.19}ÞL9:END¿
IF P==8 THEN {164.79246,304.88,44.971,0.00859,30.06896,1.769,131.722,62.2,-6.87}ÞL9:END¿
IF P==9 THEN {247.687,238.929,224.067,0.24881,39.48169,17.142,110.303,3.2,-1}ÞL9:END¿
IF P==10 THEN {557.2,21.649,187.301,0.44177,67.67,44.187,35.87,3.3,-1.1}ÞL9:END ¿
END:¿
F9(360/365.242191*I/L9(1)+L9(2)-L9(3))ÞN:¿
L9(4)ÞK:¿
RUN ".AstTrAn":¿
F9(N+L9(3))ÞC:CÞM9(P+2,15):¿
L9(5)*(1-K²)/(1+K*COS(N))ÞK:¿
SIN(C-L9(7))*COS(L9(6))ÞQ:¿
COS(C-L9(7))ÞX:¿
F9(ARG((X,Q)))ÞX: ¿
X+L9(7)ÞX:¿
ASIN(SIN(C-L9(7))*SIN(L9(6))ÞQ:¿
Ö(V²+K²-2*V*K*COS(C-F)*COS(Q)ÞM9(P+2,12): ¿
K*COS(Q)ÞB:¿
IF P<3 THEN ¿
F9(180+F+ATAN(B*SIN(F-X)/(V-B*COS(F-X))))ÞC:¿
ATAN(B*TAN(Q)*SIN(C-X)/(V*SIN(X-F)))ÞB:¿
ELSE¿
F9(ATAN(V*SIN(X-F)/(B-V*COS(X-F)))+X)ÞC:¿
ATAN(B*TAN(Q*SIN(C-X)/(V*SIN(X-F)))ÞB:¿
END:¿
(1+COS(C-X))/2ÞM9(P+2,11):¿
M9(P+2,12)ÞX:¿
L9(8)/XÞM9(P+2,10):¿
5*LOG(K*X/ÖM9(P+2,11))+L9(9)ÞM9(P+2,9): ¿
RUN ".AstEc2E":¿
RÞM9(P+2,1):EÞM9(P+2,2):¿
RUN ".AstRSet":¿
IFTE(ABS(H)<90,F0(H),H)ÞM9(P+2,5):AÞM9(P+2,7):¿
F0(S)ÞM9(P+2,6):ZÞM9(P+2,8): ¿
END:¿
 		Calculates the positions at the current date and time of all the planets, from Mercury to Neptune, plus dwarf planets Ceres, Pluto and Eris. (The International Astronomical Union met in Prague in 2006 and controversially decided that Pluto is no longer a fully-fledged planet but instead is grouped with Ceres and the Kuiper Belt object Eris as a 'Dwarf Planet'.)
The parameters calculated for each planet are its right ascension, declination, magnitude, size, phase, distance from the Earth and current altitude and azimuth. The rising and setting times and azimuths are also calculated, though without any allowance for the movement of the planets over the day so these values will be slightly inaccurate for the inner planets.
Since there are too many results to return individually they are put into the matrix M9.
Perturbations in planetary orbits by the other planets are not taken into account by this routine, which can cause errors of a few minutes of arc in the calculated positions, especially for Jupiter and Saturn. However the accuracy is considered sufficient for the intended purpose of the program.
On entry Y=year, M=month, D=date in month, T=universal time (decimal hours), L0(1)=observer's longitude, L0(2)=observer's latitude.
Returns M9 as a matrix of 10 rows and 14 columns. Row 1 is for the Sun, row 2 is for the Moon. Rows 3 to 10 are for the planets Mercury to Pluto, in order outwards. The meaning of each column is:
Column  1  right ascension (decimal hours)
Column  2  declination (decimal degrees)
Column  3  current altitude (decimal degrees)
Column  4  current azimuth (decimal degrees, North=0)
Column  5  rise time (decimal hours, local time)
Column  6  set time
Column  7  azimuth at rise
Column  8  azimuth at setting
Column  9  magnitude (i.e. visual brightness)
Column 10  angular size (seconds of arc)
Column 11  phase (the fraction of the side facing the earth illuminated, range 0-1)
Column 12  distance from Earth (in AU)
Column 13  time the planet culminates, i.e is highest in the sky (decimal hours, local time)
Column 14  maximum altitude, at culmination
Column 15  Heliocentric longitude

Program .AstTrAn


FNROOT(DEG®RAD(X)-K*SIN(X)- DEG®RAD(N),X,N)ÞX:¿
2*ATAN(Ö((1+K)/(1-K))* TAN(X/2))ÞN:¿
 		Calculates the 'true anomaly' from the 'mean anomaly' and eccentricity.
To a first approximation planets could be considered to move in circular orbits around the Sun. The angle which the line from the Sun to such a planet would make at a given date with the direction of the vernal equinox is called its mean anomaly. In reality planets move in elliptical orbits with a known eccentricity, which is a measure of how much it deviates from a circle. The true anomaly is the angle between the real planet and the vernal equinox, as seen from the Sun. Unfortunately the true anomaly cannot be calculated exactly but is found from the mean anomaly by a numeric approximation.
The simplest approximation is n=M+(360/p)*e*sin(M)
where M is the mean anomaly, e is the eccentricity and n is the true anomaly. This is accurate to a few minutes of arc if the eccentricity is small, but is not really accurate enough for planets with a high eccentricity such as Mercury.
A much more accurate value for the true anomaly can be found by first solving Kepler's equation for the eccentric anomaly E, given by:
E-e*sin(E)=M  where angles are measured in radians. This equation has to be solved numerically.
The true anomaly is then calculated as:
n=2*arctan(Ö((1+e)/(1-e))*tan(E/2))
Routine .AstTrAn uses the calculator's built-in root finder to solve the equation for eccentric anomaly and converts it to the true anomaly.
On entry N=mean anomaly, K=eccentricity.
Returns N=true anomaly

Program .AstMoon


RUN ".AstJDN": J-2447891.5ÞI:¿
360/365.242191*IÞF:¿
F9(F+279.403303-282.768422)ÞV:¿
F9(F+360/p*0.016713*SIN(V)+279.403303)ÞF:¿
F9(13.1763966*I+318.351648)ÞC:¿
F9(C-0.1114041*I-36.34041)ÞB:¿
F9(318.510107-0.052953*I)ÞN:¿
1.2739*SIN(2*C-2*F-B)ÞQ:¿
0.1858*SIN(V)ÞK:¿
B+Q-K-0.37*SIN(V)ÞB: ¿
6.2886*SIN(B)ÞZ:¿
C+Q+Z-K+0.214*SIN(2*B)ÞC:¿
(1-0.0549²)/(1+0.0549*COS(B+Z))ÞZ: ¿
FÞP:¿
C+0.6583*SIN(2*C-2*F)ÞF:¿
(1-COS(F-P))/2*-SIGN(F-P-180*(SIGN(F-P)-1)-180)ÞP:¿
N-0.16*SIN(V)ÞN:¿
SIN(F-N)*COS(5.145396)ÞQ:¿
COS(F-N)ÞX:¿
F9(ARG((X,Q))+N)ÞC:¿
ASIN(SIN(F-N)*SIN(5.145396)ÞB:¿
RUN ".AstEc2E":¿
0.9507/ZÞN:¿
RUN ".AstPara":¿
1865/ZÞN: ¿
LOG(Z²/P²)/LOG(2.51)-12.73ÞS: ¿
 		Calculates the position (right ascension and declination) of the Moon on the current date, plus other lunar data as detailed below.
One difference for calculations relating to the Moon rather than the Sun or planets is that the Moon is relatively close to the Earth and therefore the observer's precise location makes a difference to the apparent position of the Moon against the background stars. This parallax effect is very small for the planets and the directly calculated coordinates, which strictly are only correct for an observer at the centre of the Earth, are accurate enough. For the Moon though parallax can make a difference of about one degree of arc in the Moon's apparent position and this is taken into account by calling the routine .AstPara.
On entry Y=year, M=month, D=date in month, T=universal time (decimal hours).
Returns R=right ascension, E=declination, Z=distance from Earth (in multiples of the Moon's average distance, not AU as for the planets), N=angular diameter (seconds of arc), S=approximate magnitude, P=phase (fraction illuminated, -1 to +1)
Note: The figure for the visual magnitude may be somewhat inaccurate because the Moon appears much brighter when close to full than a day or two either side. There are two reasons for this:
  1. The rough surface of the Moon means that away from full small areas of the surface are in shadow from nearby higher areas, making the average brightness less. At full the shadows disappear as seen from Earth.
  2. The fine dust covering the Moon's surface preferentially scatters light back the way it has come so that the surface is considerably brighter when the light source is directly behind the observer than when it is slightly to one side.

If the Moon is waxing the phase is given as positive while if it is waning the phase is shown as negative.

Program .AstPara


RUN ".AstU2LS":¿
(L-R)*15ÞL:¿
ATAN(0.996647*TAN(L0(2)))ÞX:¿
0.996647*SIN(X)+L0(3)/6378140* SIN(L0(2))ÞQ:¿
COS(X)+L0(3)/6378140* COS(L0(2))ÞX:¿
1/SIN(N)ÞN:¿
ATAN(X*SIN(L)/(N*COS(E)-X* COS(L)))ÞI:¿
F8(R-I/15)ÞR:¿
ATAN(COS(L+I)*(N*SIN(E)-Q)/ (N*COS(E)*COS(L)-X))ÞE:¿
 		Corrects right ascension and declination, which are correct as seen from the centre of the Earth, for geocentric parallax, i.e. that the observer is actually on the surface of the Earth. The apparent direction of a nearby astronomical object depends on the observer's latitude and longitude, very slightly on the observer's altitude above sea level, and also on the fact that the Earth is not a perfect sphere but is flattened along a line between the two poles. This routine corrects for parallax.
On entry N=horizontal parallax (decimal degrees), R=geocentric right ascension (decimal hours), E=geocentric declination (decimal degrees), L0(1)=observer's longitude, L0(2)=observer's latitude, L0(3)=observer's altitude (metres).
Returns R=corrected right ascension, E=corrected declination
Horizontal parallax is the parallax angle when the object is on the observer's horizon, when the angle is greatest. It is approximately the angle subtended by the radius of the Earth as seen from the object.

Program .AstMoRS


TÞO: ¿
RUN ".AstMoon":¿
RUN ".AstRSet":¿
IF ABS(H)<90 THEN HÞT:¿
RUN ".AstMoon":¿
RUN ".AstRSet":¿
HÞG:AÞW:¿
SÞT:¿
RUN ".AstMoon":¿
RUN ".AstRSet":¿
F0(S+N/108E3)ÞS: WÞA:¿
F0(G-N/108E3)ÞH: OÞT:¿
END:¿
HÞM9(2,5):SÞM9(2,6): AÞM9(2,7):ZÞM9(2,8):¿
 		Calculates the times and azimuths of moonrise and moonset on the current day.
As explained in the section for .AstSRS this is complicated by the fact that we need to know the precise position of the Moon at the time of moonrise before we can work out the time of moonrise. As before this is solved by first calculating the approximate times of rising and setting and then recalculating the Moon's positions at these times to obtain more accurate times. A correction is also made for the Moon's finite diameter.
On entry Y=year, M=month, D=date in month.
Returns H=time of moonrise (decimal hours, local time), S=time of moonset, A=azimuth at rising (decimal degrees), Z=azimuth at setting.
These four values are also entered into matrix M9 row 2, columns 5 to 8.
Remember that moonset may be before moonrise on a particular day. If the Moon is circumpolar from the observer's location then H>+90 while if it never rises then H<-90.

Program .AstDisp


0ÞXmin:3.9ÞXmax:0ÞYmax:7ÞYmin:ERASE:¿
1-CÞQ:¿
CASE¿
IF W==0 THEN DISPXY Q;0;1;"Sunday":END¿
IF W==1 THEN DISPXY Q;0;1;"Monday":END¿
IF W==2 THEN DISPXY Q;0;1;"Tuesday":END¿
IF W==3 THEN DISPXY Q;0;1;"Wednesday":END¿
IF W==4 THEN DISPXY Q;0;1;"Thursday":END¿
IF W==5 THEN DISPXY Q;0;1;"Friday":END¿
IF W==6 THEN DISPXY Q;0;1;"Saturday":END ¿
END:¿
DISPXY 2.2-C;0;1;Y:¿
IF M£9 THEN DISPXY 2.8-C;0;1;"0"M:
ELSE DISPXY 2.8-C;0;1;M:END:¿
IF D£9 THEN DISPXY 3.2-C;0;1;"0"D:
ELSE DISPXY 3.2-C;0;1;D:END:¿
2ÞFormat:DISPXY 3.6-C;0;1;®HMS(F0(T))" Local time": ¿
IF ABS(K)<90 THEN DISPXY 6-C;0;1;"Twilight start "®HMS(K):
DISPXY 8.5-C;0;1;"Twilight end "®HMS(N):END:1ÞFormat: ¿
IF K>90 THEN DISPXY 6-C;0;1;"Twilight all night":END:¿
IF K<-90 THEN DISPXY 6-C;0;1;"No twilight":END:¿
DISPXY 11-C;0;1;"Long "INT(L0(1))"°"ABS(INT(FRAC(L0(1))*60)"'":
DISPXY 13-C;0;1;"Lat "INT(L0(2))"°"ABS(INT(FRAC(L0(2))*60)"'":¿
DISPXY 15-C;0;1;"GMT+ "L0(4)" hours":¿
IF R==1 THEN DISPXY 0;2;1;"Sun":
DISPXY 0;3;1;"Moon":DISPXY 0;4;1;"Mercury":
DISPXY 0;5;1;"Venus":DISPXY 0;6;1;"Mars":END:¿
IF R==6 THEN DISPXY 0;2;1;"Ceres":
DISPXY 0;3;1;"Jupiter":DISPXY 0;4;1;"Saturn":
DISPXY 0;5;1;"Uranus":DISPXY 0;6;1;"Neptune":END:¿
IF R==11 THEN DISPXY 0;2;1;"Pluto":DISPXY 0;3;1;"Eris":END: ¿
DISPXY 2-C;0.9;1;"Ri Asc":DISPXY 3-C;0.9;1;"Declin":¿
DISPXY 4-C;0.9;1;"Alti":DISPXY 5-C;0.9;1;"Azim":
DISPXY 6-C;0.9;1;"Rises":DISPXY 7-C;0.9;1;"Sets":
DISPXY 8-C;0.9;1;"Rise Az":DISPXY 9-C;0.9;1;"Set Az":
DISPXY 10-C;0.9;1;"Mag":DISPXY 11-C;0.9;1;"Ang Dia":
DISPXY 12-C;0.9;1;"Phase":DISPXY 13-C;0.9;1;"Dist AU":
DISPXY 14-C;0.9;1;"Culmin":
DISPXY 15-C;0.9;1;"Max Alt Helio L": ¿
LINE 0;1.7;4;1.7:LINE 0.9;1.7;0.9;7: ¿
FOR P=0 TO 1+3*(R<11);¿
IF C<3 THEN M9(P+R,1)ÞQ:
DISPXY 2-C;P+2;1;INT(Q)"h"INT(FRAC(Q)*60)"m":
M9(P+R,2)ÞQ:
IF Q<0 THEN DISPXY 3-C;P+2;1;"-":ABS(Q)ÞQ:
ELSE DISPXY 3-C;P+2;1;"+":END:
DISPXY 3.15-C;P+2;1;INT(Q)"°"INT(FRAC(Q)*60)"'":END:¿
IF C£3 THEN M9(P+R,3)ÞQ:
IF Q<0 THEN DISPXY 4-C;P+2;1;"-":ABS(Q)ÞQ:
ELSE DISPXY 4-C;P+2;1;"+":END:
DISPXY 4.14-C;P+2;1;INT(Q)"°"INT(FRAC(Q)*60)"'":
M9(P+R,4)ÞQ:
DISPXY 5-C;P+2;1;INT(Q)"°"INT(FRAC(Q)*60)"'":END:¿
IF ABS(M9(P+R,5))<90 THEN IF C³3 AND C£5
THEN 2ÞFormat:DISPXY 6-C;P+2;1;®HMS(M9(P+R,5)):
DISPXY 7-C;P+2;1;®HMS(M9(P+R,6)):1ÞFormat:END:¿
IF C³5 AND C£7 THEN M9(P+R,7)ÞQ:
DISPXY 8-C;P+2;1;INT(Q)"°"INT(FRAC(Q)*60)"'":
M9(P+R,8)ÞQ:
DISPXY 9-C;P+2;1;INT(Q)"°"INT(FRAC(Q)*60)"'":END:¿
ELSE IF C³3 AND C£5 THEN IF M9(P+R,5)>90 THEN
DISPXY 6-C;P+2;1;"Circumpolar":
ELSE DISPXY 6-C;P+2;1;"Never rises":END:END: ¿
END:¿
IF C³7 AND C£9 THEN
DISPXY 10-C;P+2;1;ROUND(M9(P+R,9),1):
DISPXY 11-C;P+2;1;ROUND(M9(P+R,10),1)"¨":END:¿
IF C³9 AND C£11 THEN
DISPXY 12-C;P+2;1;ROUND(M9(P+R,11)*100,1)"%":
DISPXY 13-C;P+2;1;ROUND(M9(P+R,12),2):END:¿
IF C³11 THEN 2ÞFormat:
DISPXY 14-C;P+2;1;®HMS(M9(P+R,13)):
1ÞFormat:M9(P+R,14)ÞQ:DISPXY 15-C;P+2;1;INT(Q)"°"INT(FRAC(ABS(Q))*60)"'":
M9(P+R,15)ÞQ:
DISPXY 16-C;P+2;1;INT(Q)"°"INT(FRAC(ABS(Q))*60)"'":
END: ¿
END: 		This long routine displays all the calculated data for the Sun, Moon and planets on the calculator's screen in tabular format. It also shows the date and time for which the calculations apply, the times of start and end of astronomical twilight, and the observer's longitude and latitude on the top line of the display.
The data which appears in the table has been stored in matrix M9. The limited size of the LCD display means that only five rows (for five planets) and three columns of data can be shown at once so parameters are passed to the routine to define which row and column appears in the top left corner.
When the display has been printed the routine returns to the calling program, which is responsible for ensuring there is a pause to read the display and for scrolling around the table.
The display format used for values in the table is to convert the decimal values to degrees and minutes or hours and minutes, as follows:
  • Times are shown as HH.MM, eg. 14.23. These are all local times.
  • Right ascensions are in the form HHhMMm, e.g. 19h45m. Leading zeros are omitted.
  • Angles (declination, altitude, azimuth) appear as DDD°MM', e.g. -65°24'.
On entry M9() holds the data to display, Y=year, M=month, D=date in month, T=local time, W=day of week, K=local time of morning twilight, N=local time of evening twilight, R=row of M9 to begin display at, C=column of M9 to begin display at.
Does not return anything useful.
Note that for a proper display C must be an odd number between 1 and 13, R must be 1, 6 or 11.

Program .AstCulm


FOR P=1 TO 10;¿
M9(P,1)ÞL:¿
RUN ".AstLS2U":¿
F0(U)ÞM9(P,13):¿
ASIN(SIN(M9(P,2))*SIN(L0(2))+ COS(M9(P,2))*COS(L0(2))) ÞM9(P,14):¿
END: 		Calculates the time of culmination and maximum altitude reached by the Sun, Moon and planets on the current day.
An astronomical object culminates, or transits, when it crosses the observer's central meridian, i.e. when it is due south in the northern hemisphere, or due north in the southern hemisphere. This is also the time when it reaches its greatest altitude above the horizon so is often the best time to view the object.
Transit occurs when the object's right ascension is equal to the local sidereal time, so by setting LST to RA and converting to universal time the transit time is easily found. The altitude at culmination is also calculated, though without any correction for refraction. If the particular object never rises then the altitude shown will be negative. Times for the Sun and Moon may be a few minutes out because no allowance is made for their movement in right ascension between the current time and transit time.
On entry M9(P,1)=right ascension, M9(P,2)=declination, Y=year, M=month, D=date in month. Returns M9(P,13)=time of culmination (local time, decimal hours), M9(P,14)=altitude at culmination (degrees).

Program .AstEast


INT(Y/100)ÞK:¿
(19*(Y MOD 19)+K-INT(K/4)-
INT((K-INT((K+8)/25)+1)/3)+15) MOD 30ÞQ:¿
Y MOD 100ÞK:¿
(32+2*(INT(Y/100) MOD 4)+2*INT(K/4)-Q-
(K MOD 4)) MOD 7ÞN:¿
Q+N-7*INT(((Y MOD 19)+11*Q+22*N)/451)+114ÞK:¿
INT(K/31)ÞN:¿
(K MOD 31)+1ÞK:¿
IF N==3 THEN MSGBOX "Easter Sunday: "K" March "Y:¿
ELSE MSGBOX "Easter Sunday: "K" April "Y:¿
END:¿
 		Calculates the date of Easter Sunday in the currently set year, i.e. the year which is used for planetary calculations.
Easter, commemorating the resurrection, is the most important date in the year for the Christian religion. Unfortunately none of the Gospels record exactly what the date was, other than that it was after the Jewish feast of the Passover. The Council of Nicaea, called by emperor Constantine in 325 AD, decreed that Easter Sunday should be the first Sunday after the first full moon after the vernal equinox, so long as this did not coincide with the start of Passover, in which case Easter was to be the following Sunday. The vernal equinox is the date when the Sun appears to pass from south to north of the celestial equator, and is always within a day of March 21st.
This definition of Easter was problematic because it depends on the interaction between the solar ('tropical') year of 365.24219 days and the synodic month (the time for one cycle of lunar phases) of 29.53059 days. Of course neither of these were known to high accuracy in 325 AD which made predicting the date of Easter in advance, without actually observing the full moon, tricky.
The cycle of dates of Easter does eventually repeat, but only after 532 years.
This program uses a method devised in 1876 to calculate the date of Easter for any year in the Gregorian calendar, i.e. from 1583 onwards. It will be wrong for years before 1583.
On entry Y=year.
Returns N=month of Easter (3=March, 4=April), K=date within month of Easter Sunday.
Shows result as a message box.

Program .AstJul


RUN ".AstJDN":¿
MSGBOX "Date "Y"-"M"-"D" Local T. "ROUND(®HMS(F0(T)),2)" Julian Day N° "J:¿ 		Displays the Julian Day Number of the currently set date and time. Note that the correct offset from Greenwich Mean Time needs to stored in L0(4) to convert local civil time to Universal Time, since Julian Day numbers are based on UT.
The routine simply calls the Julian Day Number calculating program listed earlier.
On entry Y=year, M=month, D=date in month, T=Universal Time as a decimal number of hours, i.e. in (hours+minutes/60) form.
Returns nothing.
Shows result as a message box.