To: Bob Vance

From: Ron Maddalena

Subject: Basic Astronomical Formulae in H316 Computer

The following describes the H316 code used for calculating such quantities as
aberration, nutation, etc. The view I have taken is NOT to compare the results
of such calculations with the J2000 standard but to the standard from which the
code was originally taken (i.e., Newcomb et al formulations). Thus, the results
of the H316 routines do not produce the most current suggested values but, for
the 140-ft purposes, produces accurate enough results. Since the code uses
approximations to the Newcomb way of doing things, I have checked the code for
errors which will be generated by these approximations as we head into the 21st
century. I have checked whether constants have the correct values as well as
determined how often we should update these constants. I have also studied the
code for algorithm errors, such as overflows and pathological problems but I
have probably missed a few which will only rear their ugly heads when the code
falls flat on its face.

In the following, I have given references as to where one may look for
discussions of the various subtopics, locations in the H316 code where the
calculations are performed, the algorithm used by the H316, some notes for the
programmer, and the values and locations of constants used by the algorithms.
(Note that the locations I give for variables and code are from what the code
was like on this date and may not reflect the code you are going to modify.
Also, watch out for typos in the constants I give)

Those items which I think are important for us to consider now, in the future,
or on a yearly basis are highlighted with asterisks

References:

Meeus, pp 61-62

Explanatory Supplement, p 30

Location in code: (03000 - 03226)

Algorithm: Uses an approximation to Newcomb's formulation

T0 = (initial epoch - 1900) / 100 (usually 0.5)

T = (final epoch - initial epoch) / 100

zeta = (2304".250 + 1".396*T0)*T + 0".302*T^2 + ignored terms

z = zeta + 0".791*T^2 + ignored terms

theta = (2004".682 - 0".853*T0)*T - 0".426*T^2 + ignored terms

Notes:

1. Errors caused by ignoring terms ammounts to 0".03 in alpha and delta

2. Coefficients need never be changed

3. Algorithm overflows in 4 centuries

4. Initial epoch cannot be < 1900; final - initial cannot be > 400 or < 0

Constants:

KZETA1 (03212) = 2304.250 = 0.001 777 971 turns

KZETA2 (03214) = 1.396 = 0.000 001 077 turns

KZETA3 (03216) = 0.302 = 0.000 002 330 turns

KZETT1 (03220) = 0.791 = 0.000 000 2330 turns

KTHET1 (03222) = 2004.682 = 0.001 546 823 turns

KTHET2 (03224) = 0.853 = 0.000 000 658 turns

KTHET3 (03226) = 0.426 = 0.000 000 329 turns

K1900 (03200) = 1900

K100 (03210) = 100

References: None

Location in code: (03227 - 03310)

Algorithm: Calculates number of days which have passed since a specified day with a specified Julian date (the specified day now happens to be 0h UT Jan 1, 1983)

jcjnd = Julian date of 0h UT Jan 1, 1983

ut-adj = 5h (difference in hours between EST and UT)

JD = 365*(year - 1983) + jcjnd + int[ (year-1983+2)/4 ] + day-of-year + ut-adj

Notes:

1. Note the way the number of leap days is calculated which depends upon the year chosen to start counting. JD's are stored as their real value minus 2400 000 (i.e., JD for 12h Jan 0 1990 is stored as 47892.0 as opposed to its actual value of 1447892.0.

2. Will NOT take care of the problem that the year 2000 is NOT a leap year. So, code may produce a wrong JD on March 1, 2000 and on Jan 1, 2001.

3. To change date from which JD is calculated to, for example, 2001:

a. change all J1983 references to J2001

b. update value of J2001 and JCJND

c. change line in code (near Loc 03257) which reads ADD =2 to ADD =4 (this is the number of years to the next leap year + 1)

Constants:

J1983 (03311) = 1983 (year of JCJND)

JCJND (03306) = 454 335.5 = Julian date of Jan 1, 0h UT 1983

J365 (03304) = 365 (days in one year)

UTADJ (03310) = 5 hours = 0.208333333 turns = difference in Long. between Greenwich and Green Bank.

References:

Meeus, p. 77

Smart, p. 238

Explanatory Supplement, p. 38

Location in code: (see Constants below)

Algorithm: H316 stores the obliquity (e) as a series of constants which
contain the values of cos(e), sin(e), tan(e), 20".4958*cos(e) {where
20".4958 = aberration constant}

Notes:

1. Every year e should change by 0".5 each year, so we need to update these constants each year. e should be for the middle of the year in question although its value for the beginning of the year is also OK.

2. e can be derived with sufficient accuracy (0".05) with the following:

T = (JD at 0h UT on July 1 of year in question - 2415020.) / 36525e = 23.452 294 - 0.013 0125*T - 0.000 001 64*T^2 + 0.000 000 503*T^3

4. e can also be gotten from the Astronomical Ephemeris, Section B, page B18. For 1990.5, e = 23.440527

Constants:

KCEPS (03412) = cos(e) = 0.917 473 48 in 1990

KSEPS (03410) = sin(e) = 0.397 796 95 in 1990

ATEPS (03544) = tan(e) = 0.433 578 69 in 1990

ACC1 (03540) = 20."4958*cos(e) = 0.000 091 166 075 918 radians in 1990

(All of these constants need to be changed every year; 20."4958 = aberration constant and never needs changing -- see below)

References:

Meeus, p. 68

Explanatory Supplement, pp. 47-48

Lang, pp. 500-501

Smart, pp. 179-193

Location in code: (03502 - 03544)

Algorithm: Brakes up the calculation of aberration into three parts and returns the values associated with each part. Uses aberration constant (20".4958), cos(e), and tan(e), where e is the obliquity.

term1 = 20".4958*sin(long. of Sun)

term2 = 20".4958*cos(e)*cos(long. of Sun)

term3 = term2*tan(e) = 20".4958*sin(e)*cos(long. of Sun)

Notes:

1. e needs to be changed every year... (See Obliquity of the Ecliptic)

2. 20".4958 need not be to so many digits... 20".496 would be enough. The constant of aberration changes too slowly to worry about.

Constants:

ACD1 (03542) = 20".4958 = 0.000 009 936 644 radians

ATEPS, ACC1 (see obliquity of the ecliptic above)

References:

Meeus, pp 65-66

Smart, pp 231-235

Explanatory Supplement, pp. 41-45

Location in code: (03312 - 03412)

Algorithm: Uses an adequate approximation to the long-period terms in Nutation

T = (JD - 2415020)

omega = long. of moon's ascending node = 259.1833 deg - T * 1934.1420 deg / 36525 + ignored terms

Nut. in Long. = -17".2327*sin(omega) + ignored terms

Nut. in Obliquity = 9".2100*cos(omega) + ignored terms

Notes:

1. Errors introduced by ignoring higher terms amounts to only 0".09 in Omega and 0".136 and 0".45 in Long. and Obliq. for 1990. These errors are sufficiently small that we need never to update this algorithm.

Constants:

KLCSA1 (03406) = -17".2327 = -0.000 083 546 487 turns

KOCCA1 (03404) = 9."2100 = 0.000 044 651 34 turns

KJD1900 (03376) = 15020.0 (Julian date of 1900 Jan 0 - 2400 000)

KOMEGA1 (03400) = 259.1833 = 0.719 953 5417 turns

KOMEGA2 (03402) = 1934.1420 deg / 36525 = 0.000 147 094 23 turns

References:

Meeus, p. 75

Explanatory Supplement, p. 98

Location in code: (See constants below)

Algorithm: H316 uses stored values of the eccentricity. It is stored as ecc/pi
and 0.625*ecc

Notes:

1. Since ecc. changes by .2% each year, we must change it yearly.

2. Ecc. can be gotten from the Astronomical Ephemeris, Section C, page C1 (using d = 0.5 days) or, with sufficient accuracy (0.3%) from the following:

T = (JD at 12h UT 0 Jan - 2415020.0) / 36525

ecc = 0.016 751 04 - 0.000 0418*T - 0.000 000 126*T^2

For 1990, ecc = 0.016 712 54

Constants:

L10PIE (03500) = ecc / pi = 0.005 320 014 514 for 1990

L625E (03476) = 0.625*ecc = 0.010 445 824 072 for 1990

References:

Meeus, pp. 75-77

Explanatory Supplement, pp 98-103

Location: (03414 - 03500)

Algorithm: Uses an approximation to the algorithm in Meeus and in the Expl.Suppl.

T = number of days since a particular date = (JD - JD at 12h UT 0 Jan 1983) = (JD - JD(1983) )

L = Sun's Mean Long. = L1 + L2*T + ignored terms

M = Sun's Mean Anomaly = M1 + M2*T + ignored terms

C = Sun's Equation of Center = L10PIE*sin(M) + L625E*sin(2*M)*L10PIE + ignored terms

L10PIE = eccentricity / pi

L652E = 0.625*eccentricity

L0 = Sun's Geometric Longitude = L + C

Notes:

1. The above equations for L and M work well throughout a year but, over a few years, significant errors (few arcsec) will accumulate. L2 and M2 change insignificantly with time so we need never change them over the next few centuries but L1 and M1 do change from year to year. In order to make the errors in L and M less than 1 arcsec over the course of a year, we must change the "particular date" in the calculation of T as well as the L1 and M1 constants every year. Also, we need to change L10PIE and L625E as pointed out in the eccentricity section above.

2. L10PIE and L625E are from the eccentricity calculation above (see eccentricity above)

3. L1, L2, M1, M2 can be gotten from section C, page C1 of the Astronomical Almanac (use d = 0.5). For 1990, they give L1 = 279.896 269 deg, L2 = 0.985 647 36 deg, M1 = 357.129 887 deg, M2 = 0.985 600 28 deg.

4. L1 and M1 can also be calculated with good accuracy (< 1") from:

T = (JD at 12h UT 0 Jan - 2415020.0) / 36525

L1 = 279.696 68 + 36 000.768 92*T + 0.000 3025*T^2

M1 = 358.475 83 + 35 999.049 75*T - 0.000 150*T^2 - 0.000 0033*T^3

The values of L2 and M2 are:

M2 = 35 999.049 75 deg / 36525

L2 = 36 000.768 92 deg / 36525

5. Ignored terms are on the order of 2 arcsec, which is barely allowable.

6. All references to LE1983 must be changed every year to LEyear (e.g. LE1990). These references now occur at locations 03422 and 03474.

Constants:

LMLC1 (03470) = L1 = 0.777 489 636 turns for 1990

LMLC2 (03472) = L2 = 0.002 737 909 26 turns (need never be changed)

LMAC1 (03464) = M1 = 0.992 027 464 turns for 1990

LMAC2 (03466) = M2 = 0.002 737 778 52 turns (need never be changed)

L10PIE (see eccentricity above)

L625E (see eccentricity above)

LE1983 (should now be LE1990) (03474) = -(JD at 12h UT on Jan 0 - 2400 000) = -47892 for 1990

References:

Lang, pp 535 - 536, 494-497

Methods of Experimental Physics, pp 308-314

Location in code: (22000 - 22356)

Algorithm: Uses an exact formulation

Vearth = Ve*sin(hour angle)*cos(dec. at date)*cos(geocentric lattitude)

Ve = equatorial rotation of the earth

Notes:

1. The latitude is stored as Cos(lat.); HOWEVER, it is in error! What is stored is the geodetic latitude and NOT the geocentric. The geocentric latitude of the 140-ft telescope is 38 deg 15' 01".28 so the cos(lat.) should be 0.785 313 089. The error in latitude generates a 0.25% error in Vearth (which amount to 0.0012 km/sec). As will be shown, this is trivial in comparison to other velocity errors (see note 2). (See below the section entitled "Other Problems".)

The geocentric latitude can be calculated from the following:f = 1 / 298.25 = Earth's flattening factor

tan(geocentric lat.) = (1-f)^2 * tan(geodetic lat.)

2. The velocity shift caused by the motion of the Earth-Moon barycenter is ignored, although this amounts to a 0.1 km/sec error in velocities! Also ignored are planetary perturbations which can amount to 0.01 km/sec.

Thus, unless we go through a lot of trouble adding these two additional corrections to velocity corrections, we need to only calculate vel. to within 0.1 km/sec. We easily achieve that with this and the algorithms which follow.

Constants:

VEARTH (22320) = Ve = 0.465 km/sec

CLAT (22262) = cos(lat) = 0.783 280 722 (actually should be the value [0.785 313 089] )

References:

Lang, p. 535

Methods of Experimental Physics, pp 308-314

Location in code: (22000 - 22356)

Algorithm: Uses a rearranged version of eq. 5-131 in Lang.

VSun = VSUN*cos(Lat.)*[sin(Long. Sun)*cos(Long)-cos(Long. Sun)*sin(Long)] +

VSUN*eccentricity*cos(Lat.)*[sin(Long. Per.)*cos(Long)-

cos(Long. Per.)*sin(Long)]

Lat, Long = Ecliptic coordinates of positon of source

Long. Sun = Longitude of the Sun (calculated using above algorithm)

obliquity = 23.452 294 - 0.013 0125*T + ignored terms (used in the calculation of Lat and Long from alpha and delta)

eccentricity = 0.016 751 04 - 0.000 0418*T + ignored terms

Long. Per. = Longitude of the Sun at Perigee = 281.220 833 + 1.719 175*T + ignored terms

T = (JD - 2415020) / 36525

VSUN = orbital velocity of the Sun

Notes:

1. Ignored terms amounts to a 1".3 error in Long. of Per., 3" in obliquity, and a 0.005% error in eccentricity. These amount to 0.00019, 0.0005, and 0.0015 km/sec errors in heliocentric velocities. Note that the calculation for the Sun's latitude also introduces an error of 2" or 0.00029 km/sec.

2. Although eccentricity was entered elswhere as constants, the heleocentric velocity correction calculates it anew. Replacing this calculation of ecc. with the value for the year (and updating that value every year) will reduce the error in VSun by a factor of ten.3. Like eccentricity, the obliquity is calculated here while a value is entered as a constant elswhere in the program. And, again, replacing the calculation of obliq. with a constant (updated every year) would further reduce the error in VSun.

4. The total error in heliocentric velocities is 0.004 km/sec which is small in comparison to the error introduced by ignoring Barycenter motion (0.1 km/sec).

5. The value of V, as given by eq. 5-133 in Lang is wrong! The P value he used is wrong! The actual value of V should be 29.76 km/sec.

6. While the value in Eq. 5-133 is incorrect, the equation is correct so V changes with time as the eccentricity changes. But the change in V only amounts to 0.0006 km/sec in 50 years so we can ignore it.

Constants:

JULCEN (22342) = 36525

JD1900 (22340) = 15020 = Julin date for 0 Jan 1900 - 2400 000.0

OBLC2 (22336) = 0.013 0125 deg = -0.000 036 146 turns

OBLC1 (22334) = 23.452 294 deg = 0.065 145 2611 turns

ECCC1 (22324) = 0.016 751 04

ECCC2 (22326) = -0.000 0418

PGEEC1 (22330) = 281.220 833 deg = 0.781 168 9814 turns

PGEEC2 (22332) = 1.719 175 deg = 0.004 775 4861 turns

VSUN (22316) = -29.79 km/sec

References:

Lang, p. 535, 525-527

Methods of Experimental Physics, pp 308-314

Location in code: (22000 - 22356)

Algorithm: Uses an exact formulation

Vlsr = VLSR*(VXS'*cos(alpha)*cos(delta) + VYS'*sin(alpha)*cos(delta) +

VZS'*sin(delta)

alpha, delta = Apparent coordinates

VXS', VYS', VZS' = apparent rectangular coordinates of Sun's proper motion

VLSR = Velocity of the Sun with respect to LSR

Notes:

1. VXS', VYS', VZS' are calculated from the values of VXS, VYS, and VZS which are the rectangular coordinates of the solar motion for 1950 and from the Julian date of the observations (i.e., they are precessed and then used)

VXS = cos(a0)*cos(d0); VYS = sin(a0)*cos(d0); VZS = sin(d0) where a0, d0 = 1950 coordinates of the solar motion. (a0 = 18h 01m 55.1s d0 = 30 deg 00' 04" which are 18h and 30 deg for 1900 but precessed to 1950).

2. Although Lang gives VLSR as 19.5 km/sec, it is almost universally assigned a value of 20 km/sec.

3. No additional errors in velocities are introduced by this calculation so the error quoted in Note 4 of Heliocentric velocity (0.004 km/sec) also applies here. Again, this is swamped by the 0.1 km/sec error possibly introduced by ignoring the Barycenter motion.

Constants:

VXS (22302) = 0.007 245 511

VYS (22304) = -0.865 984 917

VZS (22306) = 0.500 017 405

VLSR (22352) = -20.0 km/sec

I have noted that the latitude of the telescope is in error by 2 arcsec in
the calculation of refraction and in the matrix for converting to and from
horizontal coordinate systems. The stored value is: 38 deg 26' 17" while
the actual geodetic latitude of the 140-ft is 38 deg 26' 15".409. This
error is probably OK in these cases. (I got the 140-ft positions from the Dept.
of the Air Force survey of 1970.)

Meeus, Jean, "Astronomical Formulae for Calculators," (1980)

Smart, W. M., "Textbook on Spherical Astronomy," (Cambridge Univ, 1977).

Lang, Kenneth R., " Astrophysical Formulae," (Springer Verlag, 1980).

"Explanatory Supplement to the Astronomical Ephemeris ... ," (1961).

Gordon, M. A., in "Methods of Experimental Physics, Vol. 12, Astrophysics, Part C: Radio Observations" (Academic Press, 1976)