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)