Overview
-The first 3 functions ( KG , OB , KPL ) are M-Code routines.
-The other ones are focal programs.
-The origin of time is 2000/01/01 0h TT ( 12 hours before J2000.0
)
-Flag F04 must be cleared for Gregorian dates
-Flag F04 must be set for Julian dates.
-The longitudes are measured positively Eastward from the meridian
of Greenwich
| XROM | Function | Desciption |
| 06,00
06,01 06,02 06,03 06,04 06,05 06,06 06,07 06,08 06,09 06,10 06,11 06,12 06,13 06,14 06,15 06,16 06,17 06,18 06,19 06,20 06,21 06,22 06,23 06,24 06,25 06,26 06,27 06,28 06,29 06,30 06,31 06,32 06,33 06,34 06,35 06,36 06,37 06,38 06,39 06,40 06,41 06,42 06,43 06,44 06,45 06,46 06,47 06;48 06,49 06,50 06,51 06,52 06,53 06,54 06,55 06,56 06,57 06,58 06,59 06,60 06,61 06,62 06,63 |
-ASTRO2010
KG OB KPL "J2000" "DT" "EE" "R-S" "S-R" "MST" -EPHEMERIS "EPH" "SUN" "MO" "ME" "VE" "MA" "JU" "SA" "UR" "NE" "PL" "N" "K" "O" "L" -SATELLITES "SATJSU" -COMETS "PHEM" "GAUSS" -TRANSCOOR "EQ-ECL" "ECL-EQ" "EQ-AZ" "AZ-EQ" "PREC" "PREQ" "IOO" -REFRACTION "H-H0" "H0-H" "REFR" -MISCELLAN "BOSS" "BSTAR" "COSMO" "DELTAT" "DT-KY" "EASTER" "ECLIPSE" "FIX" "GNB+1" "GNBP" "IFMD" "KY-DT" "LST" "MNODE" "NUT" "OBL" "PHASE" "RTS" "STYLM" "ZODIAC" |
Section Header
Gaussian Gravitational Constant Obliquity of the Ecliptic on 2000/01/01 0h Kepler's Equation Date >>> Number of days since 2000/01/01 0h Number of days since 2000/01/01 0h >>> Date Equatorial <> Ecliptic ( subroutine ) Rectangular >>> Spherical Spherical >>> Rectangular Mean Sidereal Time at Greenwich Section Header Initialization The Sun The Moon Mercury Venus Mars Jupiter Saturn Uranus Neptune Pluto "L" "N" "O" "K" are 4 subroutines that are used to compute the position of celestial bodies. "N" may be useful for itself, however. Section Header Satellites of Jupiter(4), Saturn(7) & Uranus(5) Section Header Parabolic, Hyperbolic, Elliptic Motions Preliminary Orbit Determination by Gauss' Method Section Header Equatorial > Ecliptic Ecliptic > Equatorial Equatorial > Azimuthal Azimuthal > Equatorial Precession ( ecliptic coordinates ) Precession ( equatorial coordinates ) Orbital elements from one equinox to another one Section Header True Altitude > Apparent Altitude ( mean conditions Apparent Altitude > True Altitude of temp & press ) A more general refraction program Section Header Barycenter of the Solar System Binary Stars Redshift > Cosmological Distances Calculation of Delta T = TT-UT for a given year Julian/Gregorian Calendar > Hindu Calendar Date of Easter Solar & Lunar Eclipses Positional fix from the heights of 2 stars Gravitational (n+1)-body Problem Gravitational n-body Problem Illuminated Fraction of the Moon's Disk Hindu Calendar > Julian/Gregorian Calendar Local Solar Time Passage of the Moon through the Nodes Nutation ( simplified formulae ) Obliquity of the Ecliptic Phases of the Moon Rising, Transit, Setting Sidereal & Tropical years and Lunar Months Date & Time for a given longitude of the Sun |
Gaussian Gravitational Constant
KG simply returns k = 0.01720209895 in
X-register.
Obliquity of the Ecliptic on 2000/01/01 0h
Likewise, OB returns 23.4392796
in X-register,
which is the mean obliquity of the ecliptic ( expressed
in degrees ) on 2000/01/01 0h TT
This M-Code routine uses Newton's method to solve Kepler's equation:
M = E - e sin E
if 0 <= e < 1 ( elliptic
case )
M = e sinh E - E if
e > 1 ( hyperbolic case )
E = eccentric anomaly,
where M = mean anomaly
e = eccenticity
| STACK | INPUTS | OUTPUTS |
| Y | e | e |
| X | M | E |
| L | / | M |
Examples:
e = 1.2 ENTER^
M = 0.4 XEQ "KPL" returns
E = 0.987853767
e = 0.9 ENTER^
M = 7° XEQ "KPL"
returns E = 40.46693453°
Notes:
-When e < 1 , the angles are expressed in degrees, whatever the angular
mode is.
-When e < 1 , first take M modulo 360 to get a better accuracy.
-In some pathological cases - like M = 7° , e = 0.999 -
more than 40 iterations are needed before the solution ( 52°27026153
) is found.
Date >>> Number of days since 2000/01/01 0h
-"J2000" takes a date expressed in YMD format ( Gregorian Calendar
if CF 04 , Julian Calendar if SF 04 )
and calculate the number of days since 2000/01/01 at 0h. ( i-e
the Julian Day number minus 2451544.5 )
-Z & T-registers are lost.
-Y register is saved in Y & Z registers, and the day of the week
is returned in T-register ( Sunday = 0 , Monday = 1 , Tuesday = 2 , ....
, Saturday = 6 ).
Data Registers: /
Flag: F04
SF04 for the Julian Calendar
CF04 for the Gregorian Calendar
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| T | / | dow |
| Z | / | Y |
| Y | Y | Y |
| X | YYYY.MNDDdd | N |
Examples:
CF 04 2134.0404 ( Gregorian ) XEQ "J2000"
>>> N = 49036 RDN
RDN gives 0 = Sunday
SF 04 1234.0428 ( Julian
) XEQ "J2000" >>>
N = -279651 RDN RDN gives 5 =
Friday
Notes:
-The days may have decimals
-You can choose the Julian Calendar after 1582/10/15
and you can choose the Gregorian Calendar before 1582/10/04
-This program works for all dates: positive, 0 or negative.
-The "B.C." years are counted astronomically:
year 0 = 1 B.C.
year -1 = 2 B.C. ...etc...
Number of days since 2000/01/01 0h >>> Date
"DT" performs the reciprocal conversion.
Data Registers: /
Flag: F04
SF04 for the Julian Calendar
CF04 for the Gregorian Calendar
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| T | / | Y |
| Z | / | Y |
| Y | Y | Y |
| X | N | YYYY.MNDDdd |
Examples:
CF 04 N = 49036 XEQ "DT"
>>> 2134.0404
SF 04 N = -279651 XEQ "DT" >>>
1234.0428
Note:
-"J2000" & "DT" are separated by a RTN, so you can easily use them
to convert a date in the Julian Calendar into a date in the Gregorian Calendar.
-For instance: SF 04 1234.0428
XEQ "J2000" >>> -279651 CF 04 R/S >>>
1234.0505
-Thus 1234 April 4th Julian = 1234 May 5th Gregorian.
Equatorial <> Ecliptic ( subroutine )
-Many transformations of coordinates use the same type of formulae which appear in the equatorial-ecliptic conversion, namely:
sin b = cos e sin d - sin e cos d sin a
cos b cos l = cos d cos a
cos b sin l = sin e sin d +
cos e cos d sin a
-So, "EE" is extensively used in the "TRANSCOOR" section.
Data Registers: /
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Z | e | e |
| Y | d | b |
| X | a | l |
Example: if a = 116°328942 , d = 28°026183 and e = 23°4392911 ( in DEG mode )
23.4392911 ENTER^
28.026183 ENTER^
116.328942 XEQ "EE" >>>> l
= 113°215630 RDN b =
6°684170
Note:
-Like "R-S" and "S-R" , "EE" works in all angular modes.
x = r cos b cos l
-We use the formulae: y = r cos b sin l
with the built-in R-P function.
z = r sin b
where x , y , z = rectangular coordinates,
r = ( x2 + y2 + z2 )1/2
, l = longitude , b = latitude
Data Registers: /
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| T | T | T |
| Z | z | b |
| Y | y | l |
| X | x | r |
| L | / | (x2+y2)1/2 |
Example: x = 3 ; y = 4 ; z = -7 find the spherical coordinates ( in DEG mode )
-7 ENTER^
4 ENTER^
3 XEQ "R-S" r
= 8.602325267
RDN l
= 53.13010235°
RDN b = -54.46232221°
"R-S" works in all angular modes.
-The same formulas are used with P-R
Data Registers: /
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| T | T | T |
| Z | b | z |
| Y | l | y |
| X | r | x |
| L | / | (x2+y2)1/2 |
Example: r = 10 ; l = 124° ; b = 37° find x ; y ; z ( in DEG mode )
37 ENTER^
124 ENTER^
10 XEQ "S-R" x = -4.465913097
RDN y = 6.620988446
RDN z = 6.018150232
-"S-R" works in all angular modes
Mean Sidereal Time at Greenwich
"MST" computes the mean sidereal time at Greenwich
• T = the number of millenia since 2000/01/01 0h is stored
into R00
Data Register: R00
= Number of millenia since 2000/01/01 0h
Flag: /
Subroutine: "J2000"
| STACK | INPUTS | OUTPUTS |
| Y | YYYY.MNDD | / |
| X | hh.mnss (UT) | HH.MNSS (MST) |
Example: Find the mean sidereal time at Greenwich on 2010/07/16 at 7h41m UT
2010.0716 ENTER^
7.41
XEQ "MST" >>>> MST = 3h17m10s
and R00 = T = 0.01053886417
"EPH" stores your longitude in R16 measured positively
eastwards
from the meridian of Greenwich.
your latitude in R17 positive = North , negative = South
the value of DeltaT = TT-UT in R18 expressed in seconds,
( in 2010 , TT-UT is about 66 seconds )
"EPH" also calculates MST in R19 and stores the number of millenia since 2000/01/01 0h TT in R00
• Suppose we are at the US Naval Observatory at Washington ( D.C. )
Longitude = 77°03'56" W = -77°03'56"
Latitude = 38°55'17"
N = +38°55'17"
• It is 7h41m UT , on July 16th 2010
• And we take TT-UT = 66 seconds
XEQ "EPH" the HP41 displays "LONG^LAT"
-77.0356 ENTER^
38.5517 R/S
next PROMPT "TT-UT=SEC?"
66 R/S next PROMPT "DATE^TIME?"
2010.0716 ENTER^
7.41
R/S the HP-41
returns the date in Y & the time in X
>>> R16-R17-R18-R19 are initialized, R00 = 0.01053886626
( DeltaT has been added to R00 after executing "MST" ) and the HP-41 is
set in DEG mode.
We suppose these values remain unchanged
in all the examples below
-We are now at the beginning of "SUN"
-All the programs "SUN" "MO" "ME" "VE" "MA" "JU" "SA" "UR" "NE" "PL"
return
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
R03 = geocentric longitude in decimal degrees = l
R04 = geocentric latitude ------------------
= b
R05 = distance to the Earth in Astronomical Units = r ( except
for the Moon: R05 = the Moon's parallax in sexagesimal degrees )
R06 = Right Ascension in hh.mnss = R.A
R07 = Declination in
° ' '' = Decl
R08 = heliocentric longitude in decimal degrees = L ( except
for the Sun and the Moon )
R09 = heliocentric latitude ------------------- =
B ---------------------------------
R10 = radius vector in Astronomical Units
= R ---------------------------------
R11 = Azimuth in ° ' '' reckoned clockwise from
the South = Az
R12 = True Altitude in ° ' '' = h
R13 = Apparent Altitude in ° ' '' ( after correction
for refraction ) = h0
R14 = Elongation from the Sun = Elg
R01 = X & R02 = Y are the rectangular coordinates of the Sun:
"SUN" may be executed first to get the exact geocentric or topocentric
coordinates of the planets.
The effects of aberration and nutation are neglected: we obtain
mean
coordinates.
The correction for parallax is neglected too, except for the
Moon's altitude ( it may reach 1 degree )
These ephemerides have an accuracy of about 0.01° in
the heliocentric longitudes over the interval 1000-3000.
( except for Pluto: its coordinates are obtained with this precision
over 1880-2110 only )
The Sun's coordinates are more accurate. However, errors in the
geocentric coordinates are increased when the distance between the Earth
and the planet
is small and in this case, for Venus and Mars, the error in
the obtained geocentric longitude can reach 1 arcminute, perhaps a little
more.
All these coordinates are referred to the mean ecliptic and
equinox of the date.
Notes:
-If R00 is correctly initialized but R16-R17-R18-R19 are not, the topocentric
coordinates will be wrong,
but the right ascension, declination, heliocentric & geocentric
longitudes & latitudes and the distances from the Sun or the Earth
will be exact !
-If CF 01, all the coordinates described above are computed
-If SF 01, the stack registers X & Y contain the Sun's radius
vector & its longitude, without disturbing R01, R02, ....
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
CF 01 XEQ "SUN" yields Az = 216°27'26"
= R11
RDN h = -20°37'58" = R12
RDN h0 = -20°37'58" = R13
( when h0 < 0 , h0 is meaningless )
RDN Elg = 0 ... logical !
R03 = l = 113°692
R06 = R.A. = 7h42m14s
R04 = b = 0
R07 = Decl = 21°21'37"
R05 = r = 1.016434 AU
and the rectangular coordinates of the Sun: R01 = X = -0.408421 & R02 = Y = 0.930760
• SF 01 XEQ "SUN" returns X = 1.016434
& Y = 113°692
-If you've just executed "SUN" with CF 01, simply press R/S instead
of XEQ "MO"
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
XEQ "MO" >>>>
Az = 145°13'58" = R11
RDN h = -49°32'56" = R12
RDN h0 = -49°32'56" = R13
( when h0 < 0 , h0 is meaningless )
RDN Elg = 62°383
R03 = l = 175°957
R06 = R.A. =11h37m06s
R04 = b = -5°073
R07 = Decl = -3°02'57"
R05 = parallax = 0°59'48"
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
XEQ "ME" >>>>
Az = 198°51'53" = R11
RDN h = -30°14'17" = R12
RDN h0 = -30°14'17" = R13
( when h0 < 0 , h0 is meaningless )
RDN Elg = 18°528
R03 = l = 132°157
R06 = R.A. = 9h00m19s
R08 = L = 186°609
R04 = b = 1°556
R07 = Decl = 18°38'33"
R09 = B = 4°687
R05 = r = 1.19459 AU
R10 = R = 0.39700 AU
-If you've just XEQ "ME" , simply press R/S instead of
XEQ "VE"
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
XEQ "VE" >>>>
Az = 171°34'59" = R11
RDN h = -40°29'28" = R12
RDN h0 = -40°29'28" = R13
( when h0 < 0 , h0 is meaningless )
RDN Elg = 42°846
R03 = l = 156°525
R06 = R.A. = 10h34m53s
R08 = L = 229°255
R04 = b = 1°184
R07 = Decl = 10°12'59"
R09 = B = 1°570
R05 = r = 0.96043 AU
R10 = R = 0.72394 AU
-If you've just XEQ "VE" , simply press R/S instead of
XEQ "MA"
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
XEQ "MA" >>>>
Az = 151°20'12" = R11
RDN h = -43°12'11" = R12
RDN h0 = -43°12'11" = R13
( when h0 < 0 , h0 is meaningless )
RDN Elg = 58°114
R03 = l = 171°803
R06 = R.A. = 11h30m57s
R08 = L = 204°193
R04 = b = 0°676
R07 = Decl = 3°52'18"
R09 = B = 0°796
R05 = r = 1.8976 AU
R10 = R = 1.6113 AU
-If you've just XEQ "MA" , simply press R/S instead of
XEQ "JU"
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
XEQ "JU"
>>>> Az = 315°49'15" = R11
RDN h = 41°41'30" = R12
RDN h0 = 41°42'35" = R13
... at last a planet above the horizon !
RDN Elg = 110°370
R03 = l = 3°316
R06 = R.A. = 0h14m20s
R08 = L = -7°751
R04 = b = -1°361
R07 = Decl = 0°04'11"
R09 = B = -1°238
R05 = r = 4.519 AU
R10 = R = 4.965 AU
-If you've just XEQ "JU" , simply press R/S instead of
XEQ "SA"
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
XEQ "SA"
>>>> Az = 140°54'06" = R11
RDN h = -41°17'24" = R12
RDN h0 = -41°17'24" = R13
( when h0 < 0 , h0 is meaningless )
RDN Elg = 65°925
R03 = l = 179°597
R06 = R.A. = 12h02m08s
R08 = L = 185°186
R04 = b = 2°268
R07 = Decl = 2°14'28"
R09 = B = 2°356
R05 = r = 9.904 AU
R10 = R = 9.535 AU
-If you've just XEQ "SA" , simply press R/S instead of
XEQ "UR"
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
XEQ "UR"
>>>> Az = 319°25'07" = R11
RDN h = 42°40'41"
= R12
RDN h0 = 42°41'43"
= R13
RDN Elg = 113°153
R03 = l = 0°537
R06 = R.A. = 0h03m11s
R08 = L = -2°129
R04 = b = -0°765
R07 = Decl = -0°29'18"
R09 = B = -0°749
R05 = r = 19.672 AU
R10 = R = 20.093 AU
-If you've just XEQ "UR" , simply press R/S instead of
XEQ "NE"
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
XEQ "NE"
>>>> Az = 2°07'56" = R11
RDN h = 38°30'15" = R12
RDN h0 = 38°31'27" = R13
RDN Elg = 145°508
R03 = l = -31°819
R06 = R.A. = 22h02m04s
R08 = L = -32°918
R04 = b = -0°473
R07 = Decl = -12°32'59"
R09 = B = -0°460
R05 = r = 29.176 AU
R10 = R = 30.019 AU
-If you've just XEQ "NE" , simply press R/S instead of
XEQ "PL"
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
XEQ "PL"
>>>> Az = 55°29'06" = R11
RDN h = 10°53'39" = R12
RDN h0 = 10°58'25" = R13
RDN Elg = 159°289
R03 = l = -86°413
R06 = R.A. = 18h15m03s
R08 = L = 274°218
R04 = b = 5°081
R07 = Decl = -18°18'35"
R09 = B = 4°929
R05 = r = 30.912 AU
R10 = R = 31.865 AU
-This subroutine takes e in Y-register and M
in X-register ( M = mean anomaly , e = eccentricity < 1 here )
-And it returns v in Y-register and r
/ q in X-register where
v = true anomaly
r = radius vector
q = perihelion distance
-This is used, for instance, to study the motion of a comet, an asteroid....
on an elliptic orbit.
| STACK | INPUTS | OUTPUTS |
| Y | e < 1 | v |
| X | M | r/q |
Example:
e = 0.9 ENTER^
M = 7° XEQ "N"
>>>> r / q = 3.152974273
X<>Y v = 116°2027755
Note:
-The HP-41 must be in DEG mode.
Satellites of Jupiter, Saturn and Uranus
-This program allows you to get the configuration of
4 satellites of Jupiter ( Io - Europa - Ganymede - Callisto
)
7 satellites of Saturn ( Mimas - Enceladus - Tethys - Dione
- Rhea -Titan - Hyperion )
5 satellites of Uranus ( Ariel - Umbriel -Titania - Oberon -
Miranda )
-The x-axis coincides with the equator of the planet.
-The center of the Planet is the origin and x , y are measured in equatorial
radii.
y ( North )
|
|
|
( East ) ------------Jup/Sat/Ura------------------ x (
West )
|
|
( South )
Data Registers: R00 thru R20 are used for temporay data storage and when the program stops:
Mimas - Enceladus - Tethys - Dione - Rhea - Titan - Hyperion
R01 = x1 R03 = x2
R05 = x3 R07 = x4
R09 = x5 R11 = x6
R13 = x7
R02 = y1 R04 = y2
R06 = y3 R08 = y4
R10 = y5 R12 = y6
R14 = y7
or Io - Europa - Ganymede - Callisto
or Ariel - Umbriel - Titania - Oberon - Miranda
>>>> R19 = - sin DE where DE is the planetocentric declination of the Earth.
Flags: F01 F02
F03 F04 F05 F06 F07
-Flag nn is set when the distance Earth-Satellite n is shorter
than the distance Earth-Planet ( 0 < nn < 8 )
Subroutines: "MST" "J2000"
"SUN" "JU" "SA" "UR"
| STACK | INPUTS | OUTPUTS |
| Z | YYYY.MNDD | / |
| Y | HH.MNSS(TT) | y1 |
| X | 1 or 2 or 3 | x1 |
X = 1 for the Satellites
of Jupiter
X = 2 ------------------
Saturn
X = 3 ------------------
Uranus
Example:
On 2009/07/16 at 12h41m TT
• Jovian Satellites
2009.0716 ENTER^
12.41
ENTER^
1 XEQ "SATJSU"
>>>> x1 = -1.075
X<>Y y1 = -0.053
and in registers R01 thru R08:
Io - Europa - Ganymede - Callisto
x1 = -1.075 x2
= 9.261 x3
= 13.120 x4
= -12.053 DE
= 0°52
y1 = -0.053 y2
= 0.019 y3
= -0.059 y4
= 0.223
-Flags F01 and F03 are set: Io & Ganymede are closer to the
Earth than Jupiter.
• Saturnian Satellites
2009.0716 ENTER^
12.41
ENTER^
2 R/S
>>>> x1 = -2.426
X<>Y y1 = -0.058
and in registers R01 thru R14:
Mimas - Enceladus - Tethys - Dione - Rhea - Titan - Hyperion
x1 = -2.426 x2
= -0.945 x3
= 2.034 x4 = -5.336
x5 = -8.749 x6
= -11.495 x7 = -10.289
DE
= -2°63
y1 = -0.058 y2
= -0.176 y3
= -0.112 y4 = 0.150
y5 = -0.039 y6
= 0.626 y7
= 0.764
-Flags F04 F05 F06 F07 are set,
whence Dione, Rhea, Titan and Hyperion are closer to the Earth than
Saturn.
• Uranian Satellites
2009.0716 ENTER^
12.41
ENTER^
3 R/S
>>>> x1 = 6.052
X<>Y y1 = -0.667
and in registers R01 thru R10:
Ariel - Umbriel - Titania - Oberon - Miranda
x1 = 6.052
x2 = -7.880 x3
= 5.296 x4
= 22.755 x5 = 4.557
DE = 8°76
y1 = -0.667 y2
= 1.039 y3
= -2.483 y4
= 0.265 y5 = 0.449
-Flags F01 & F03 are set, so Ariel & Titania are
closer to the Earth than Uranus.
Notes:
-The accuracy is of order of a few hundredths of the planet's radius.
-Hyperion's coordinates are less accurate than those of the other Saturnian
Satellites.
-For Hyperion, the series converge slowly and several terms should
be added to get more accurate results.
-The complete series may be downloaded from the server
ftp://ftp.imcce.fr
-For the Uranian Satellites, the error in the y-values may reach 0.1
Parabolic, Hyperbolic, Elliptic Motions
Parabolic Motion
-Let s = Tan v/2 where v = true anomaly,
s is the unique solution of the cubic equation:
s3 + 3s - W = 0
where W = (3k/sqrt(2)) (t-T) q -3/2
and k = 0.01720209895 = Gaussian gravitational constant
-Then, the radius vector is obtained by r = q(1+s2)
Hyperbolic Motion
-For hyperbolic orbits, Kepler's equation is: k(t-T)
((e-1)/q)3/2 = e Sinh E - E
-Newton's method is used in the M-Code routine KPL.
-Then, the true anomaly v is obtained by Tan v/2 = ((e+1)/(e-1))1/2
Tanh E/2
and the radius vector by r = q(1+e)/(1+e cos v)
Elliptic Motion
-We solve Kepler's equation M = E - e sin E ( in radians
) where M = k(t-T) ((1-e)/q)3/2
= mean anomaly and E = excentric anomaly
-This equation is solved by Newton's method and the M-code routine
KPL.
-Then, the true anomaly v is obtained by Tan v/2 = ((e+1)/(1-e))1/2
Tan E/2
and the radius vector is r = q(1+e)/(1+e cos v)
Data Registers: R00 = Number of millenia since 2000/01/01 0h
R01 & R02 = rectangular ecliptic coordinates of the Sun
R03 = geocentric longitude in decimal degrees = l
R04 = geocentric latitude ------------------
= b
R05 = distance to the Earth in Astronomical Units = r
R06 = Right Ascension in hh.mnss = R.A
R07 = Declination in
° ' '' = Decl
R08 = heliocentric longitude in decimal degrees = L
R09 = heliocentric latitude -------------------
= B
R10 = radius vector in Astronomical Units
= R
R11 = Azimuth in ° ' '' = Az reckoned
clockwise from the South
R12 = True Altitude in ° ' '' = h
R13 = Apparent Altitude in ° ' '' = h0 ( after
correction for refraction )
R14 = Elongation from the Sun = Elg
R16 = Longitude ( positive East ) in ° ' "
R17 = Latitude in ° ' "
R18 = Delta T in seconds
R19 = Mean sidereal time at Greenwich
R15 & R20 are unused
R21 = T = time of passage in perihelion - expressed in
days since 2000/01/01 0h TT ( not 12h )
R22 = q = perihelion distance in AU
R23 = e = eccentricity
R24 = i = inclination in degree
referred to
R25 = omega = argument of the perihelion in degree
the ecliptic
R26 = OMEGA = longitude of the ascending node in degree
of the date
Flag: /
Subroutines: "MST" "SUN" "N"
......
| STACK | INPUTS | OUTPUTS |
| T | / | Elg = elongation from the Sun ( deg ) |
| Z | / | h0 = apparent altitude ( ° ' '' ) |
| Y | / | h = true altitude ( ° ' '' ) |
| X | / | Az = Azimuth ( ° ' '' ) |
Example1: Calculate the position of comet Kohler for 1977 September 29 0h UT using the following elements
T = 1977 November 10.5659 TT
i2000 = 48°7131
q = 0.990662
omega2000 = 163°4788
e = 1
OMEGA2000 = 182°1660
i = 48°7160
-First, we reduce the elements to 1977/11/29 it gives
omega = 163°4798
use for instance "IOO"
OMEGA = 181°8582
-We have T = -8086.4341 days from 2000/01/01 0h
XEQ "PHEM" >>>> the HP-41 displays "T^Q^E?"
-8086.4341 ENTER^
0.990662
ENTER^
1
R/S the HP-41 displays
"I^O^OO?"
48.7160
ENTER^
163.4798 ENTER^
181.8582
R/S the HP-41 displays
"LONG^LAT?" and, if we were at the USNO:
-77.0356
ENTER^
38.5517
R/S next
PROMPT "TT-UT=SEC?"
49 R/S next PROMPT "DATE^TIME?"
------------------------------------------------------------------------------------------------------------------------
1977.0929 ENTER^
0
R/S Az
= 78°07'22" = R11
RDN
h = 46°36'14" = R12
RDN
h0 = 46°37'08" = R13
RDN Elg
= 62°505 = R14
R03 = l = -121°944
R06 = R.A. = 16h19m10s
R08 = L = 302°797
R04 = b = 40°960
R07 = Decl = 20°16'22"
R09 = B = 44°329
R05 = r = 1.306 AU
R10 = R = 1.225 AU
>>> To calculate its position for another date, place the date & time ( UT ) in Y & X-registers and R/S
For instance 1977.1212 ENTER^
20 R/S
yields Az = 334°24'25" ... etc ...
Example2:
Calculate the position of comet C/2007 T1 for 2008 January 01 6h
UT using the following elements
T = 2007 December 12.49731 TT
i2000 = 117°649041
q = 0.969480
omega2000 = 233°671201
e = 1.000785
OMEGA2000 = 111°418623
i = 117°64857
-Reduce the elements to 2008/01/01.25
omega = 233°67226
use for instance "IOO"
OMEGA = 111°53086
-We have T = 2902.49731 days from 2000/01/01 0h
XEQ "PHEM" >>>> the HP-41 displays "T^Q^E?"
2902.49731 ENTER^
0.969480
ENTER^
1.000785
R/S the HP-41 displays
"I^O^OO?"
117.64857 ENTER^
233.67226 ENTER^
111.53086
R/S the HP-41 displays
"LONG^LAT?" and, if we were at the USNO:
-77.0356
ENTER^
38.5517
R/S next
PROMPT "TT-UT=SEC?"
65 R/S next PROMPT "DATE^TIME?"
------------------------------------------------------------------------------------------------------------------------
2008.0101 ENTER^
6
R/S Az
= 320°13'55" = R11
RDN
h = -59°25'36" = R12
RDN
h0 = -59°25'36" = R13
RDN Elg
= 39°158 = R14
R03 = l = -99°228
R06 = R.A. = 17h02m53s
R08 = L = 219°595
R04 = b = -34°699
R07 = Decl = -57°40'52"
R09 = B = -61°144
R05 = r = 1.582 AU
R10 = R = 1.029 AU
>>> To calculate its position for another date, place the date & time ( UT ) in Y & X-registers and R/S
For instance 2008.0131 ENTER^
23 R/S
yields Az = 357°55'18" ... etc ...
Example3: Calculate the position
of comet C/2007 K6 for 2007 December 01 0h UT using the
following elements
T = 2007 July 01.47533 TT
i2000 = 105°063204
q = 3.432968
omega2000 = 337°140230
e = 0.984585
OMEGA2000 = 298°075386
i = 105°06377
-We reduce the elements to 2007/12/01 it gives omega
= 337°13934 use for instance "IOO"
OMEGA = 298°18570
T = 2007 July 01.47533 = 2738.47533 days from 2000/01/01 0h
XEQ "PHEM" >>>> the HP-41 displays "T^Q^E?"
2738.47533 ENTER^
3.432968
ENTER^
0.984585
R/S the HP-41 displays
"I^O^OO?"
105.06377 ENTER^
337.13934 ENTER^
298.18570
R/S the HP-41 displays
"LONG^LAT?" and, if we were at the USNO:
-77.0356
ENTER^
38.5517
R/S next
PROMPT "TT-UT=SEC?"
65 R/S next PROMPT "DATE^TIME?"
------------------------------------------------------------------------------------------------------------------------
2007.1201 ENTER^
0
R/S Az
= 62°27'44" = R11
RDN
h = 8°14'54" = R12
RDN
h0 = 8°21'05" = R13
RDN Elg =
38°523 = R14
R03 = l = -73°633
R06 = R.A. = 19h07m27s
R08 = L = 295°894
R04 = b = 7°068
R07 = Decl = -15°25'02"
R09 = B = 8°451
R05 = r = 4.426 AU
R10 = R = 3.706 AU
>>> To calculate its position for another date, place the date & time ( UT ) in Y & X-registers and R/S
For instance 2007.1207 ENTER^
23.4523 R/S
yields Az = 64°46'26" ... etc ...
Notes:
-When the HP-41 displays "T^Q^E?" , simply press R/S without any digit
entry if the 6 elements are already stored in R21 thru R26
-This is the case if you've just used "GAUSS" to find the orbital elements
of a new comet.
-If you have to calculate several positions, i , omega , OMEGA may
be regarded as constant for several days or weeks
( though "IOO" is relatively easy to use... )
-The program uses the ecliptic and equator of the date and that could
be criticized because this is not an inertial frame.
-However, the resulting errors are usually negligible - provided the
intervals of time remain of the order of a few days or weeks or even months.
Preliminary Orbit Determination by Gauss' Method
-Gauss' method computes approximate values of T , q , e , i ,
omega , OMEGA from 3 observations of a comet.
-This program also uses the Herrick-Gibbs formulae which are 4th-order
accurate in time.
-The method fails if the inclination i = 0 or is very small.
-At an instant t1 the right ascension of the
comet = RA1 its declination = d1 and
the rectangular ecliptic coordinates of the Sun are X1
& Y1
-At an instant t2 the right ascension of the
comet = RA2 its declination = d2 and
the rectangular ecliptic coordinates of the Sun are X2
& Y2
-At an instant t3 the right ascension of the
comet = RA3 its declination = d3 and
the rectangular ecliptic coordinates of the Sun are X3
& Y3
-We must have t1 < t2 < t3
>>> The results will be more accurate if t3 - t2 is equal ( or nearly equal ) to t2 - t1 ( this value is usually of the order of a few days )
---------------------------------------------------------FORMULAE-------------------------------------------------------------------------
-Let
li
= cos RAi cos di
mi
= sin obl sin di + cos obl cos di sin RAi
where obl is the obliquity of the ecliptic, i =
1 , 2 , 3
ni
= cos obl sin di - sin obl cos di sin RAi
-Let Di = distance Earth-Comet at the instant ti , ui ( li , mi , ni ) , -Ri = position vector of the Earth , ri = position vector of the Comet
ri = Di ui - Ri ( || ui || = 1 )
-If we neglect the planetary perturbations, the 3 position-vectors ri lie in the same plane, so there are coefficients c1 , c3 such that
r2 = c1 r1 + c3 r3 and likewise for the velocity at the second instant: V2 = -d1 r1 + d2 r2 + d3 r3
-One can prove that, approximately: ci = Ai + ( 1 + B2 / r23 ) Bi / r23 ( i = 1 , 3 ) with B2 = [ k2(t3 - t1)2/12 ] ( 1 + A1 A3 )
and di = Gi + Hi / ri3 ( i = 1 , 2 , 3 )
where A1 =
(t3 - t2)/(t3 - t1) ,
B1 = k2A1 [ (t3 - t1)2
- (t3 - t2)2 ]/6
with k = 0.01720209895 = Gaussian gravitational
constant
A3 = (t2 - t1)/(t3 - t1)
, B3 = k2A3 [ (t3
- t1)2 - (t2 - t1)2
]/6
and
G1 = (t3 - t2)/[(t3 - t1)(t2
- t1)] G3 = (t2
- t1)/[(t3 - t1)(t3 - t2)]
G2 = G1 - G3
H1 = k2(t3 - t2)/12
H3 = k2(t2 - t1)/12
H2 = H1 - H3
-Let A( A1 X1 - X2
+ A3 X3 , A1 Y1 - Y2
+ A3 Y3 , A1 Z1 - Z2
+ A3 Z3 )
and B( B1 X1 + B3
X3 , B1 Y1 + B3 Y3
, B1 Z1 + B3 Z3 )
( all the Zi's ~ 0 if we use the rectangular ecliptic coordinates
of the Sun )
P1 =
A.(
u2
x u3 )/D
Q1 = B.( u2 x u3
)/D
P2 =
A.(
u1
x u3 )/D
Q2 = B.( u1 x u3
)/D where
D = u1.( u2 x u3
) ( .
= dot product , x = cross product )
P3 =
A.(
u1
x u2 )/D
Q3 = B.( u1 x u2
)/D
-The distance Earth-Comet at the central instant is found by solving the equation
D2 = P2 + ( 1 + B2 / r23 ) Q2 / r23 with r2 = D22 - 2 D2 ( l2 X2 + m2 Y2 ) + X22 + Y22
-Then, D1 = [ P1 + ( 1 + B2
/ r23 ) Q1 / r23
] / [ A1 + ( 1 + B2 / r23 )
B1 / r23 ]
and D3 = [ P3 + (
1 + B2 / r23 ) Q3 / r23
] / [ A3 + ( 1 + B2 / r23 )
B3 / r23 ]
and the rectangular ecliptic coordinates of the comet at the 3 instants are:
xi
= li Di - Xi
yi
= mi Di - Yi
i = 1 , 2 , 3 so we know
the 3 vectors ri's and the velocity V2
= -d1 r1 + d2r2
+ d3 r3 may be computed
zi
= ni Di
-We calculate the cross-product h = r2 x V2 and then:
-The parameter p = h2/k2
-The inclination i = Acos hz/h = Atan2 (hx2+hy2)1/2/hz
you can also use C = r1 x r3
instead of h to get i and OMEGA
-The longitude of the ascending node OMEGA = Atan2 hx/(-hy)
( vi + omega ) is obtained by
ri Sin ( vi + omega ) =
zi / Sin i
ri Cos ( vi + omega ) =
xi Cos OMEGA + yi Sin OMEGA
whence (v3 - v1)/2
( v + omega ) may also be computed by the formula:
v + omega = sign(z) Acos [ ( x Cos OMEGA + y Sin OMEGA ) / r ] with perhaps less accurate results in some cases.
-The relations: 2e Sin (v3 + v1)/2
= ( p/r1 - p/r3 ) / Sin (v3 - v1)/2
give e and (v3 + v1)/2 whence v1
and omega
2e Cos (v3 + v1)/2 = ( p/r1 + p/r3
- 2 ) / Cos (v3 - v1)/2
-And the time of passage in perihelion is computed via Kepler's equation.
------------------------------------------------------------------------------------------------------------------------------------------------------
-Registers R01 thru R15 are to be initialized.
-However, if you set F00, the HP-41 displays several PROMPTs and calculate
the rectangular coordinates of the Sun X & Y
Data Registers: R00 & R16 to R42: temp ( Registers R01 thru R15 are to be initialized )
• R01 = t1
• R06 = t2
• R11 = t3
expressed in days since 2000/01/01 0h TT
• R02 = RA1
• R07 = RA2
• R12 = RA3
in hh.mnss
• R03 = d1
• R08 = d2
• R13 = d3
in °. ' ''
• R04 = X1
• R09 = X2
• R14 = X3
in Astronomical Units
• R05 = Y1
• R10 = Y2
• R15 = Y3
---------------------
>>>> When the program stops,
R21 = T ( in days from 2000/01/01 0h )
R22 = q ( in Astronomical Units )
R23 = e
R24 = i ( in degrees )
R25 = omega ( in degrees )
R26 = OMEGA ( in degrees )
R27 = n = mean motion ( in degrees/day ) - virtual for a hyperbolic
orbit.
Flag: F00
Subroutines: "SUN" "MST"
| STACK | INPUT | OUTPUT |
| T | / | i |
| Z | / | e |
| Y | / | q |
| X | D | T |
where D is an estimation of the distance Comet-Earth ( in AU ) at the 2nd instant ( take for instance D = 1 )
Example: You observe a comet on 2007/07/01 0h , 2007/07/05 0h , 2007/07/09 0h TT the positions were:
Dates 2007/07/01
0h
2007/07/05 0h
2007/07/09 0h
RA
14h27m26
14h16m34s
14h06m38s
Decl
-39°30'56"
-38°44'07"
-37°52'59"
• SF 00 1 XEQ "GAUSS" >>> "DATE^TIME1?"
2007.0701 ENTER^
0
R/S
>>> "RA^DECL1?"
14.2726
ENTER^
-39.3056
R/S
>>> "DATE^TIME2?"
2007.0705 ENTER^
0
R/S
>>> "RA^DECL2?"
14.1634
ENTER^
-38.4407
R/S
>>> "DATE^TIME3?"
2007.0709 ENTER^
0
R/S
>>> "RA^DECL3?"
14.0638
ENTER^
-37.5259
R/S
>>> T = 2817.84871 = R21 ( in
days from 2000/01/01 0h TT ) Activate the
turbo mode...
RDN q = 0.70439 = R22
RDN e = 0.76153 = R23
RDN i = 9°70768 = R24
and R25 = omega = -1°22487
R26 = OMEGA
= 3°49104 Remark that
T = 2007/09/18.84871 ( use "DT" to get this result
)
R27 = n = 0.19415
deg/day
Notes:
-"GAUSS" solves an equation by the secant method.
-This equation may have several solutions ( 2 or sometimes 3 )
-If the solution that is returned is not satisfactory, XEQ "GAUSS"
again with CF 00 and another guess in X-register.
-Sometimes, the algorithm converges to a negative - meaningless - distance
D = Comet-Earth.
-In this case, it stops and displays "D=?"
-Place a different initial guess ( say 3 ) and press R/S
-Gauss' method is very sensitive to observational errors, so take the average of successive results to increase the accuracy.
-"GAUSS" works for elliptic, hyperbolic and parabolic orbits, though
the case e = 1 exactly will be very rare...
-The orbital elements are automatically stored in the proper registers
to be used by "PHEM"
-If you execute "PHEM" with the orbital elements found above ( take
TT-UT = 0 ),
you'll find - almost - the same right ascensions & declinations
-The largest errors occur at the 2nd instant 2007/07/05 0h:
"PHEM" returns for this date: RA = 14h16m34s029
& Decl = -38°44'07"173 in registers R06
& R07
instead of
RA = 14h16m34s & Decl
= -38°44'07"
Data Registers: R00 = number
of millenia since 2000/01/01 0h
R11 & R12: temp
Flags: /
Subroutines: "J2000" "OBL"
"EE"
| STACK | INPUTS | OUTPUTS |
| Z | YYYY.MNDDdd | em |
| Y | decl ( ° . ' " ) | b ( ° . ' " ) |
| X | RA (hh.mnss) | l ( ° . ' " ) |
where l = ecliptic longitude , b
= ecliptic latitude , em = mean obliquity of the ecliptic
Example: On 2134/04/04 at 0h if Right-Ascension = 12h34m56s and Declination = 25°12'49" ( mean coordinates )
2134.0404 ENTER^
25.1249 ENTER^
12.3456 XEQ
"EQ-ECL" >>>> l = 177°13'44"69
RDN b = 26°27'18"71
Data Registers: R00 = number
of millenia since 2000/01/01 0h
R11 & R12: temp
Flags: /
Subroutines: "J2000" "OBL"
"EE"
| STACK | INPUTS | OUTPUTS |
| Z | YYYY.MNDDdd | -em |
| Y | b ( ° . ' " ) | Decl ( ° . ' " ) |
| X | l ( ° . ' " ) | RA (hh.mnss) |
where l = ecliptic longitude , b = ecliptic latitude , em = mean obliquity of the ecliptic
Example: On 2134/04/04 at 0h if l = 177°13'44"69 and b = 26°27'18"71 ( mean coordinates )
2134.0404 ENTER^
26.271871 ENTER^
177.134469 XEQ "ECL-EQ"
>>>> RA = 12h34m56s00
RDN decl = 25°12'49"00
Data Registers: R00 = number
of millenia since 2000/01/01 0h (
Registers R16 & R17 are to be initialized before executing "EQ-AZ"
)
R11 thru R15: temp
• R16 = Longitude ( positive East ) in ° ' "
• R17 = Latitude in ° ' "
Flags: /
Subroutines: "J2000"
"OBL" "MST" "EE"
| STACK | INPUTS | OUTPUTS |
| T | YYYY.MNDD | / |
| Z | HH.MNSS (UT) | h0 ( ° . ' " ) |
| Y | decl ( ° . ' " ) | h ( ° . ' " ) |
| X | RA (hh.mnss) | Az ( ° . ' " ) |
Azimuths are measured clockwise ( westwards ) from South
h = Altitude
h0 = Apparent Altitude ( corrected for refraction
)
Example: On 2005/12/12 at
20h51m29s (UT) we have R.A. = 7h41m16s & decl
= 60°21'37" ( mean coordinates ),
Compute the azimuthal coordinates at the US Naval Observatory.
• Longitude = 77°03'56" W =
-77°03'56" STO 16
• Latitude =
38°55'17" N = +38°55'17" STO 17
2005.1212 ENTER^
20.5129 ENTER^
60.2137 ENTER^
7.4116
XEQ "EQ-AZ" >>>> Az = 190°56'31"
RDN h = 10°55'58"
RDN h0 = 11°00'47"
Data Registers: R00 = number
of millenia since 2000/01/01 0h (
Registers R16 & R17 are to be initialized before executing "AZ-EQ"
)
R11 thru R15: temp
• R16 = Longitude ( positive East ) in ° ' "
• R17 = Latitude in ° ' "
Flags: /
Subroutines: "J2000"
"OBL" "MST" "EE"
| STACK | INPUTS | OUTPUTS |
| T | YYYY.MNDD | / |
| Z | HH.MNSS (UT) | / |
| Y | altitude ( ° . ' " ) | decl ( ° . ' " ) |
| X | azimuth ( ° . ' " ) | RA (hh.mnss) |
-Azimuths are measured clockwise ( westwards ) from South
Example: On 2005/12/12 at
20h51m29s (UT) we have Azimuth = 190°56'31"
& Altitude = 10°55'58" ( mean coordinates
),
Compute the equatorial coordinates at the US Naval Observatory.
• Longitude = 77°03'56" W =
-77°03'56" STO 16
• Latitude =
38°55'17" N = +38°55'17" STO 17
2005.1212 ENTER^
20.5129
ENTER^
10.5558
ENTER^
190.5631 XEQ
"AZ-EQ" >>>> R.A. = 7h41m16s
RDN decl = 60°21'37"
Precession ( ecliptic coordinates )
Data Registers: R00 & R11
thru R14: temp
Flag: F00
Subroutines: "J2000" "EE"
| STACK | INPUTS | OUTPUTS |
| T | YYYY.MNDDdd1 | / |
| Z | YYYY.MNDDdd2 | / |
| Y | b1 ( ° . ' " ) | b2 ( ° . ' " ) |
| X | l1 ( ° . ' " ) | l2 ( ° . ' " ) |
where l = ecliptic longitudes , b = ecliptic latitudes
Example: The mean ecliptic coordinates of
Sirius on 1600/04/04 0h are l1 = 98°30'58"32
; b1 = -39°39'17"79
Calculate the ecliptic coordinates on 2134/12/12 0h
1600.0404 ENTER^
2134.1212 ENTER^
-39.391779 ENTER^
98.305832 XEQ "PREC" >>>>
l2
= 105°57'44"15
RDN b2 = -39°35'19"15
-The polynomial coefficients of the precession angles ( expressed in
degrees ) are rounded to 6 decimals, with T expressed in millenia from
J2000.
-So, the accuracy is of the order of 0.01 arcsecond over the
time span [ 1000 ; 3000 ]
Precession ( equatorial coordinates )
Data Registers: R00 & R11
thru R14: temp
Flag: F00
Subroutines: "J2000" "EE"
"PREC"
| STACK | INPUTS | OUTPUTS |
| T | YYYY.MNDDdd1 | / |
| Z | YYYY.MNDDdd2 | / |
| Y | decl1 ( ° . ' " ) | decl2 ( ° . ' " ) |
| X | RA1 (hh.mnss) | RA2 (hh.mnss) |
Example: The mean coordinates of Sirius on
1600/04/04 0h are AR = 6h27m17s88 ; Decl = -16°21'56"34
Calculate the equatorial coordinates on 2134/12/12 0h
1600.0404 ENTER^
2134.1212 ENTER^
-16.215634 ENTER^
6.271788 XEQ "PREQ"
>>>> RA = 6h51m10s79
RDN Decl = -16°52'22"05
Orbital elements from one equinox to another one
Data Registers: • R00 = Date2 = YYYY.MNDDdd2 ( Register R00 is to be initialized before executing "IOO" )
R01 thru R07: temp
Flag: F00 is cleared at the end
Subroutines: "J2000"
"EE"
| STACK | INPUTS | OUTPUTS |
| T | YYYY.MNDDdd1 | i2 |
| Z | OMEGA1 | OMEGA2 |
| Y | omega1 | omega2 |
| X | i1 | i2 |
Example: We have i1
= 12°789 omega1 = 49°345 OMEGA1
= 166°234 on 1600/01/01
Compute i2 , omega2 , OMEGA2
on 2900/12/12
2900.1212 STO 00
1600.0101 ENTER^
166.234 ENTER^
49.345
ENTER^
12.789
XEQ "IOOM" >>>> i2
= 12.619940°
RDN omega2 =
49.370109°
RDN OMEGA2 = 184.401887°
True Altitude > Apparent Altitude ( mean conditions
of temp & press )
-The light coming from a star is curved by the atmosphere,
so that its apparent altitude differs from its true altitude.
-Knowing the refraction R allows to convert the apparent altitude
h0 and the true altitude h of a given star:
h = h0 - R
-"H-H0" & "H0-H" use data from the Pulkovo Refraction Tables.
Temperature:
+15°Celsius
Pressure:
1013.25 mbar
Light wave-length:
0.590 µm
Partial pressure of water vapor:
0 ( dry air )
Latitude:
45°
Observer's altitude:
0 ( i-e at sea-level )
Formula: R ~ (1°/62.6)/Tan( h + 5.459/( h + 19.272/( h + 6.942 ) ) ) where h is expressed in degrees
-Errors are smaller than 0"8 over the whole range [ -0°32'58"
, 90° ]
Data Registers: /
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Z | / | h |
| Y | / | R |
| X | h | h0 |
-All angles expressed in ° ' " ( but Z-output in degrees and decimals )
Example: With h = 1°30'00"
1.30 XEQ "H-H0" >>>> h0 =
1°48'38"6
RDN R = 0°18'38"6
Apparent Altitude > True Altitude ( mean conditions
of temp & press )
Formula: R ~ (1°/62.83)/Tan( h0 + 4.208/( h0 + 14.978/( h0 + 5.906 ) ) ) where h0 is expressed in degrees
-Errors are smaller than 0"34 over the whole range [ 0°
, 90° ]
Data Registers: /
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Z | / | h0 |
| Y | / | R |
| X | h0 | h |
-All angles expressed in ° ' " ( but Z-output in degrees and decimals )
Examples: With h0 = 1°30'00" ; h0 = 27° ; h0 = 0
1.30 XEQ "H0-H" >>>> h =
1°09'42"6 RDN R = 0°20'17"4
27 XEQ "H0-H" >>>>
h = 26°58'08"3 RDN R = 0°01'51"7
0 XEQ "H0-H"
>>>> h = -0°32'58"0 RDN
R = 0°32'58"0
A more general refraction program
"REFR" is a generalization of "H0-H"
Formulae: The standard refraction is computed by R0 ~ (57"298)/Tan( h0 + 4.208/( h0 + 14.978/( h0 + 5.906 ) ) )
R = R0 ( P/1013.25 ) [ 286.68
/ ( t + 271.68 ) ] [ 1 - f / 6579 - ( f/426 )2 ] [ 1
+ { (15-t) + (15-t)2 / 196 } exp(-h0/7) / (1+h0)
/ 157 ]
x [ 1 + (P-1013.25) exp(-0.4h0)/12181 ] [ 1 -
( f + f2/16 ) exp(-0.18h0) / (1+h0)
/ 5198 ]
where
t = Temperature
P = Pressure
f = Humidity
Data Registers: R00 = h0
( Registers R01 thru R03 are to be initialized before executing
"REFR" )
• R01 = t ( in °C )
( -10 <= t <= +30 )
• R02 = P ( in mbar )
( 700 <= P <= 1100 )
• R03 = f ( in mbar )
( 0 <= f <= 20 )
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Y | / | R |
| X | h0 | h |
( All angles are expressed in ° ' " )
Example: t = 0°C ; P = 900 mbar ; f = 12 mbar
0 STO 01 900 STO 02 12 STO 03
-If h0 = 0°
0 XEQ "REFR" >>>> h = -0°33'32"
X<>Y R = 0°33'32"
-If h0 = 10°23'45"
10.2345 R/S
>>>> h = 10°19'02"8 X<>Y
R = 0°04'42"2
-If h0 = 49°12'34"
49.1234 R/S
>>>> h = 49°11'47"9 X<>Y
R = 0°00'46"1
Barycenter of the Solar System
-This program calculates the rectangular heliocentric ecliptic coordinates
X , Y of the center of mass of the Solar system ( we make the approximation
Z = 0 )
-It takes a date YYYY.MNDD and returns the distance G between
the Solar system barycenter B and the center of the Sun S.
-Only the 4 giant planets are taken into account: the other planets
have a negligible effect.
-The positions of these planets are not computed very accurately: the
terms in e2 are neglected ( e = eccentricity ),
and only one term in the perturbations is taken into account.
Data Registers: R00 = T in millennia
since 2000/01/01 0h TT
R01 = X
R02 = Y R03 to R11: temp
Flags: /
Subroutine: "J2000"
| STACK | INPUTS | OUTPUTS |
| Y | / | µ |
| X | YYYY.MNDD | G |
where G = distance Barycenter-Center of the
Sun ( unit = solar radius )
and µ
= ( Sx , SB )
Example:
2008.0824 XEQ "BOSS"
>>>> G = 1.0761
X<>Y µ = -74°0782
Since G > 1 , the Solar system barycenter is currently outside
the Sun's globe.
G will be smaller than one in April 2010
We have also X = R01 = 0.2952
Y = R02 = -1.0348
Other Checkings: on 1983/03/15
G = 2.098 ( exact result - almost the maximum possible
value )
on 1990/04/23 G = 0.069 ( exact
value = 0.066 )
on 2130/03/10 G = 0.016 ( exact
value = 0.021 )
Note:
Errors seem smaller than 0.01 over the time span [ 1000
, 3000 ]
µ may be inaccurate if G is very small
-This program calculates the polar coordinates ( r , µ ) of a
binary star ( S1 ; S2 )
and the eccentricity e' of the apparent orbit, as seen from the
Earth.
North
|
S2 |
\ |
\ µ |
r \ |
r = the angular distance = S1S2
\ |
µ = the position angle = ( S1N , S2N
)
\|
S1
Data Registers: R07-R08: temp
• R00 = P = period of revolution ( expressed in years )
• R01 = T = time of perihelion passage ( given as a year and
decimals )
• R02 = a = semimajor axis ( in arcseconds )
• R03 = e = eccentricity
• R04 = i = inclination of the orbit
• R05 = omega = longitude of the periastron.
• R06 = OMEGA = position angle of the ascending node
Flag: F00
if CF 00 you have to initialize R00 thru R06
if SF 00 the HP-41 displays several PROMPTs
Subroutine: "N"
| STACK | INPUTS | OUTPUTS |
| Z | / | e' |
| Y | / | µ |
| X | YYYY.yy | r |
Example: The elements of Gamma Virginis are:
P = 168.68 years
T = 2005.13
i = 148.0°
a = 3.697"
omega = 256.5°
e = 0.885
OMEGA = 36.9°
-Calculate r and µ for the epoch 2010.25 ( ~ 2010/04/01 )
SF 00
2010.25 XEQ "BSTAR" >>>>
" P^T^A^E?"
168.68 ENTER^
2005.13 ENTER^
3.697 ENTER^
0.885
R/S
>>>> " I^O^OO?"
148.0 ENTER^
256.5 ENTER^
36.9
R/S
>>>> r = 1.544"
RDN µ = 19.66°
RDN e' = 0.844
-Likewise 2015 R/S
>>>> r = 2.352"
RDN µ = 4.92°
( RDN e' = 0.844 )
Redshift > Cosmological Distances
-Knowing the cosmological redshift z of a galaxy, and
the values of the current parameters (Omega)mat , (Omega)lambda
, (Omega)rad
"COSMO" computes several "distances", the recessional velocity,
the scale factor and the age of the Universe.
D
= light-travel time distance ( in light-years ) = light-travel time
( in years )
D0 = comoving
distance ( in light-years )
DL = luminosity-distance
( in light-years )
VR = recessional
velocity ( speed of light = 1 )
t0
= age of the Universe ( in years )
k
= curvature ( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic
)
Formulae:
t0 = (c/H0)
§01 y. [ (Omega)lambda.y4
+ ( 1-(Omega)tot ).y2 + (Omega)mat.y +
(Omega)rad ] -1/2 dy
D = (c/H0) §y(em)1
y. [ (Omega)lambda.y4 + ( 1-(Omega)tot
).y2 + (Omega)mat.y + (Omega)rad
] -1/2 dy
y(em) = y at the instant of emission
D0 = (c/H0) §y(em)1
[ (Omega)lambda.y4 + ( 1-(Omega)tot
).y2 + (Omega)mat.y + (Omega)rad
] -1/2 dy
z + 1 = R0/R = 1/yem
where (Omega)tot = (Omega)mat + (Omega)lambda + (Omega)rad
F(x) = Sinh(x) if k = -1
DL = R0 ( z + 1 ) F(D0/R0)
where F(x) = x
if k = 0
F(x) = Sin(x) if k = +1
VR = H0.D0
-These integrals are calculated by Gauss-Legendre 3-point formula (
the change of variable y = u2 facilitates the numerical
integration )
-The number of subintervals is doubled until 2 consecutive rounded
values are equal.
-Thus, the accuracy is controlled by the display settings
( avoid FIX 9 and SCI 9 , it could lead to infinite loops ! )
Data Registers: R00 thru R17 ( Registers R11-R12-R13-R14 are to be initialized before executing "COSMO" )
• R11 = (Omega)mat
density parameter ( matter )
• R12 = (Omega)lambda
cosmological constant
• R13 = (Omega)rad
radiation parameter
• R14 = 1/H0
= inverse of Hubble's constant expressed in years ( for instance
13.7 109 )
-when the program stops: R00 = k
( +1 = Spherical , 0 = Euclidean , -1 = Hyperbolic )
R01 = D , R02 = D0 , R03 = DL , R04 = VR
, R05 = t0 , R06 = R0 ( current scale factor )
Flags: F00 set flag F00
if you don't want to (re-)calculate the age of the Universe t0
F10
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| T | / | VR |
| Z | / | DL |
| Y | / | D0 |
| X | z | D |
| L | / | t0 |
Example1: With (Omega)mat = 0.044 , (Omega)lambda = 0.519 , (Omega)rad = 0.000049 & 1/H0 = 13.7 109
0.044 STO 11
0.000049 STO 13
0.519 STO 12
13.7 E9 STO 14
and, for instance, FIX 4
• z = 7 XEQ "COSMO" >>>>
D = 1.4046 E10 l-y
= R01
RDN D0 = 3.4012 E10 l-y
= R02
RDN DL = 4.1177 E11 l-y
= R03
RDN VR = 2.4826
= R04 ( unit c = 1 )
LASTX t0 = 1.5507 E10 years
= R05
and R00 = k = -1 ( hyperbolic space )
R06 =
R0 = 2.0725 E10 light-years
-For another redshift - in the same Universe - SF 00, place z in X-register
and R/S
Example2: With
(Omega)mat = 0.271 , (Omega)lambda
= 0.732 , (Omega)rad = 0.000049
& 1/H0 = 13.7 109
( these parameters are suggested by WMAP = NASA's Wilkinson Microwave Anisotropy
Probe )
0.271 STO 11
0.000049 STO 13
0.732 STO 12
13.7 E9 STO 14
and, for instance, FIX 7
• z = 7 XEQ "COSMO" >>>>
D = 1.2823287 E10 l-y
= R01
RDN D0 = 2.8607721 E10
l-y = R02
RDN DL = 2.2835499 E11
l-y = R03
RDN VR = 2.0881542
= R04 ( unit c = 1 )
LASTX t0 = 1.3596706 E10
years = R05
and R00 = k = +1 ( spherical space )
R06 =
R0 = 2.4810862 E11 light-years
Notes:
-Actually, WMAP suggests that the cosmological parameters lead to a
flat space ( k = 0 ), but this is still uncertain...
-In most cases, (Omega)rad has a negligible effect,
except in the very early age of the Universe.
- t0 cannot be computed if there is no Big-Bang.
-Set F00 in this case.
-The value 1/H0 = 13771721429 years would produce - almost
- exactly: Hubble's constant = H0 = 71 km/s/Mpc
-These numbers are linked by the relation: ( 1/H0
) ( years ) = 977.79222143 109 / H0 ( km/s/Mpc
)
-We have supposed that our Universe is isotropic & homogeneous,
which is only an approximation.
-See reference [7] for a more general point of view.
Calculation of DeltaT = TT-UT for a given year
-This routine computes the difference between Terrestrial Time ( or
Ephemeris Time ) and Universal Time for a given year y.
-It uses several polynomial approximations, all expressed in seconds:
y <
2
>>>> DeltaT ~ 32 ( y - 1820 )2
10 -4 - 20
2 < y < 1600
>>>> DeltaT ~ 1587 - 556 t + 66.7
t2 + 0.33 t3 - 0.523 t4
- t5/291 + t6/533 with t
= y/100 - 10
1600 < y < 1800 >>>>
DeltaT ~ 8 - 16 t + 96 t2 - 75 t3 - 38
t4 + 38 t5
with t = y/100 - 17
1800 < y < 1998 >>>>
2 polynomials are used. See reference [1] page 80.
1998 < y < 2150 >>>>
DeltaT ~ 64 + 35 t - 112 t2 + 444 t3
- 874 t4 + 46.4 t5
with t = y/100 - 20
y > 2150
>>>> DeltaT ~ 32 ( y - 1820 )2
10 -4 - 20
| STACK | INPUTS | OUTPUTS |
| X | year | Delta T (sec) |
Examples:
-2000 XEQ "DELTAT" >>>> 46675.68
s
400 XEQ "DELTAT" >>>>
6702.33 s
1880 XEQ "DELTAT" >>>> -5.42
s
1984 XEQ "DELTAT" >>>>
53.82 s
2010 XEQ "DELTAT" >>>>
66.80 s
2150 XEQ "DELTAT" >>>> 328.48
s
-The accuracy is of the order of 1 second between 1800 & 2010
-For other polynomial approximations, see also: http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html
Julian/Gregorian Calendar > Hindu Calendar
-"DT-KY" & "KY-DT" use the Old Hindu Calendars
that began on -3101 January 23rd ( Gregorian
)
= -3101 February 18th ( Julian )
= 0.0101 ( Kali Yuga )
( -3101 = 3102 BC )
-It is based on a mean sidereal year = 1577917500 / 4320000 =
210389 / 576 = 365 + 149/576 ~ 365.25868 days
-However, if you set flag F01, this value is replaced by a more accurate
approximation, i-e 172401/472 ~ 365.2563559 days.
( in 2000 , one sidereal year = 365.256363004 days )
-There are 12 months which may have 30 or 31 days.
Formulae: Given a date
YYYY.MNDD in the Gregorian or Julian calendar, we first compute
the number of days N since 2000/01/01 , then:
Let s = N + 1863079.25 , the date yyyy.mndd in the Hindu Calendar is obtained as follows:
y = FLOOR ( 576 s / 210389 )
m = 1 + [ FLOOR ( 6912 s / 210389 ) ] MOD
12
d = 1 + FLOOR [ s mod ( 210389 / 6912
) ]
Data Registers: R00-R01-R02: temp
Flag: F01-F04 CF 01 =
traditional old Hindu Calendar CF 04 = Gregorian
Calendar
SF 01 = "new-old" Hindu Calendar
SF 04 = Julian Calendar
Subroutine:
"J2000"
| STACK | INPUT | OUTPUT |
| X | YYYY.MNDDGr/Ju | YYYY.MNDDKy |
Examples: With CF 04 ( Gregorian Calendar )
• CF 01 Traditional Hindu Calendar
-3101.0123 XEQ "DT-KY" >>>>
0.0101
1979.0716 XEQ "DT-KY" >>>>
5080.0331
2010.0404 XEQ "DT-KY" >>>>
5110.1219
• SF 01 "New-Old" Hindu Calendar
-3101.0123 XEQ "DT-KY" >>>>
0.0101
1979.0716 XEQ "DT-KY" >>>>
5080.0412
2010.0404 XEQ "DT-KY" >>>>
5111.0101
Notes:
-In the Hindu calendar, the day begins at 6 a.m.
-So we also assume that a date written yyyy.mndd in the Julian
or Gregorian calendar represents a time after 6 a.m.
-The reverse conversion is performed by "KY-DT"
-This program takes a year in X-register and returns the date of Easter
for this year.
Data Registers: /
Flag: F04
CF 04 = Gregorian calendar
SF 04 = Julian calendar
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| X | YYYY | YYYY.MNDD |
Examples:
CF 04 2010 XEQ "EASTER"
>>>> 2010.0404
SF 04 1243
R/S
>>>> 1243.0412
-This program calculates the dates of Lunar & Solar eclipses with
an accuracy of about 5 minutes between 1000 & 3000 and probably a little
more...
Data Registers: R00 = 4.k
R05 = 0 or 1 R10 = 2.F
R15 = 7.38264721
R01 = m
R06: temp
R11 = 2.M'
R16: temp
R02 = m'
R07 = E
R12 = M'+M-180°
R03 = gamma R08 = M' -
180° R13 = M'-M-180°
R04 = phase
R09 = M
R14 = 2M
For a Solar eclipse, gamma =
the least "distance" from the center of the Earth to the axis of the Moon's
shadow ( unit = Earth's radius )
For a Lunar eclipse, gamma =
the least "distance" from the center of the Moon to the axis of the Earth's
shadow ( unit = Earth's radius )
Flag: F02 is set to indicate a Hybrid
( solar ) eclipse.
Subroutines: "J2000" & "DT"
| STACK | INPUT | OUTPUTS |
| T | / | m' |
| Z | / | m |
| Y | / | HH.MNSS |
| X | yyyy.mndd | YYYY.MNDD |
where yyyy.mndd is an approximate date
YYYY.MNDD is
the date of the eclipse
HH.MNSS
is the time of the maximum eclipse
For a Lunar Eclipse: m = magnitude
of the umbral eclipse
if m < 0 there is no umbral eclipse
m , m' > 1 for a total eclipse,
m' = magnitude of the penumbral eclipse
m , m' < 1 for a partial eclipse
For a Solar Eclipse: m
= magnitude of the eclipse
m > 1 for a total or hybrid eclipse
m < 1 for a partial eclipse
m' = 0 for a partial eclipse
m' = 1 for an annular eclipse
m' = m for a total or hybrid eclipse
F02
will be set for a hybrid eclipse
Examples:
3000.1016 XEQ "ECLIPSE" >>>> " SOLAR" is displayed. Press the backarrow key
3000.1019
Flag F02 is set, so we have a Hybrid ( Solar ) eclipse
RDN 16.0807
RDN m = 1.0027
R03 = gamma = -0.2167
RDN m' = 1
R04 = phase = 0 ( New Moon )
-The exact instant is 3000/1019 at 16h10m16s and the exact m =
1.0049
>>>> To search for the next eclipse, simply press
R/S
-In the last example, it yields: " LUNAR" press the backarrow key
3000/11/04 RDN 5h44m53s RDN m = -0.4165 RDN m' = 0.5983 R03 = gamma = -1.2355 R04 = 2 ( Full Moon )
-Since m < 0 , there is no eclipse in the umbra:
this Lunar eclipse is a penumbral one.
Positional Fix from the heights of 2 Stars
-"FIX" calculates your longitude & latitude from the altitudes of
2 stars ( or planets ... )
-There are usually 2 solutions.
-Several PROMPTs allow you to store the data in the proper registers
R01 thru R10.
-We have to solve the non-linear system:
sin b sin
d1 + cos b cos d1 cos ( ST1 - a1
+ L ) = sin h1
( ST = mean Sidereal Time at Greenwich , hi
= altitude )
sin b
sin d2 + cos b cos d2 cos ( ST2 - a2
+ L ) = sin h2
( L = longitude , b = latitude , ai = right-ascension
, di = declination )
may be solved without any iterative method:
-Let x = sin b After eliminating L from these 2 equations, we find that x is one solution of the quadratic equation A x2 + B x + C = 0
A = tan2 d1 + tan2 d2 + sin2
( ST2 - a2 - ST1 + a1 ) - 2
tan d1 tan d2 cos ( ST2 - a2
- ST1 + a1 )
with B
= -2 sin h1 sin d1 / cos2 d1
- 2 sin h2 sin d2 / cos2 d2
+ 2 ( sin h1 sin d2 + sin h2 sin d1
) cos ( ST2 - a2 - ST1 + a1
) / ( cos d1 cos d2 )
C = sin2 h1 / cos2 d1 + sin2
h2 / cos2 d2 - sin2 ( ST2
- a2 - ST1 + a1 ) - 2 sin h1
sin h2 cos ( ST2 - a2 - ST1
+ a1 ) / ( cos d1 cos d2 )
-After solving this quadratic equation, we find the longitude L by using the R-P conversion and:
2 cos [ L + ( ST2 - a2 + ST1 - a1
)/2 ] cos [ ( ST2 - a2 - ST1 + a1
)/2 ] = E + F
2 sin [ L + ( ST2 - a2 + ST1 - a1
)/2 ] sin [ ( ST2 - a2 - ST1 + a1
)/2 ] = E - F
with E = ( sin
h1 - x sin d1 ) / [ cos d1 ( 1 - x2
)1/2 ]
and F
= ( sin h2 - x sin d2 ) / [ cos d2 ( 1
- x2 )1/2 ]
Data Registers: R00 & R11 thru R19: temp
• R01 = Date1
• R06 = Date2
YYYY.MNDD
• R02 = Time1
• R07 = Time2
HH.MNSS ( UT )
• R03 = R.A1
• R08 = R.A2
hh.mnss
• R04 = Decl1
• R09 = Decl2
° . ' ''
• R05 = h1
• R10 = h2
° . ' ''
Flags: /
Subroutines: "MST" "H0-H"
| STACK | INPUTS | OUTPUTS |
| T | / | b2 ( ° ' " ) |
| Z | / | L2 ( ° ' " ) |
| Y | / | b1 ( ° ' " ) |
| X | / | L1 ( ° ' " ) |
where L = longitude ( positive East ) , b = latitude
Example: on 2006/04/04,
• 6h00m (UT) alpha Lyr
- RA = 18h37m09s , Decl = 38°47'21" , Altitude = 29°41'59"
( apparent )
• 6h01m (UT) zeta UMaj -
RA = 13h24m10s , Decl = 54°53'34" , Altitude = 60°29'56"
( apparent )
XEQ "FIX" >>>> "DATE^TIME1?"
2006.0404 ENTER^
6
R/S "RA^DECL^H1?"
18.3709
ENTER^
38.4721
ENTER^
29.4159
R/S "DATE^TIME2?"
2006.0404 ENTER^
6.01
R/S "RA^DECL^H2?"
13.2410
ENTER^
54.5354
ENTER^
60.2956
R/S It
yields:
L1 = -136°36'29''
RDN b1 =
77°27'26''
RDN L2 = -74°40'43'
RDN b2 =
25°50'00''
-The 2 solutions are ( L1 , b1 ) & ( L2 , b2 )
Notes:
-"FIX" takes the refraction into account and corrects the apparent altitudes by calling "H0-H"
-If you know that you are on the Atlantic Ocean ( near the Bahamas !
), the correct solution is ( L2 , b2 )
-Otherwise - or to increase accuracy - use another star and take the
average of the results.
-The program doesn't work if a declination = +/-90°
but in this case, the system is very easy to solve:
-If, for instance, d1 = 90° then b = h1
and after substituting this result into the 2nd equation sin
b sin d2 + cos b cos d2 cos ( ST2 - a2
+ L ) = sin h2
ACOS easily gives the 2 L-values.
-It doesn't really matter to take this case into account because no
star has a declination exactly equal to 90° or -90°
-However, roundoff errors may be important - at least in the longitude
- if the declination is about +/-89°59'
-But it still works well for the Pole Star.
We have to solve a system of 3n second order differential equations
of the type: d2y/dt2 = f (y)
The time t and the first derivative y' = dy/dt do not appear
explicitly in this trajectory problem.
So, we can use a special 4th order Runge-Kutta method, namely:
y(t+h) = y(t) + h ( y'(t) + k1/6 + 2k2/3
)
y'(t+h) = y'(t) + k1/6 + 2k2/3 + k3/6
where k1 = h.f (y) ; k2 = h.f(y+h.y'/2+h.k1/8)
; k3 = h.f(y+h.y'+h.k2/2)
Let M1( x1, y1,z1)
; ....... ; Mn( xn, yn,zn)
be n celestial bodies with initial velocities V1(
x'1, y'1,z'1) ; ....... ; Vn(
x'n, y'n,z'n) at an instant t
and m1 , ...... , mn
the n masses of these bodies.
We want to know their future positions and velocities at an instant
t + N.h ( h = step size ; N = number of steps )
d2xi /dt2 = SUMj#i
Gmj ( xj - xi ) / [ ( xi -
xj )2+ ( yi - yj )2
+ ( zi - zj )2 ]3/2
So, we have to solve the system:
d2yi /dt2 = SUMj#i
Gmj ( yj - yi ) / [ ( xi -
xj )2 + ( yi - yj )2
+ ( zi - zj )2 ]3/2
i , j = 1 , ..... , n
d2zi /dt2 = SUMj#i
Gmj ( zj - zi ) / [ ( xi -
xj )2 + ( yi - yj )2
+ ( zi - zj )2 ]3/2
where G is the gravitational constant.
"GNBP" can solve the 15-body problem.
Data Registers:
R00 thru R26+19n are used
R00 = n = number of bodies
R27+n = x1
R27+4n = x'1
R24 = G = constant of gravitation
R28+n = y1
R28+4n = y'1
R25 = h = step size
R29+n = z1
R29+4n = z'1
R26 = N = number of steps
................
...................
R27 = m1
R28 = m2
R24+4n = xn
R24+7n = x'n
...............
R25+4n = yn
R25+7n = y'n
R26+n = mn
R26+4n = zn
R26+7n = z'n
Flag: F10
Subroutines: /
>>> Thus, the n masses, the 3n position coordinates and the 3n velocity coordinates are stored into contiguous registers, starting at R27
N.B.
We use G = k2 where k = 0.01720209895
is Gaussian gravitational constant.
The units are Astronomical Unit , Solar mass and day.
| STACK | INPUT | OUTPUT |
| X | / | bbb.eee |
where bbb.eee = control number of the solution.
Example: Let's imagine a system of three stars
M1 ( 2 ; 0 ;0 ) M2
( 0 ; 4 ; 0 ) M3 ( 0 ; 0 ; 1 )
unit = 1 AU
with initial velocities V1
( 0 ; 0.03 ; 0 ) V2 ( 0 ; 0 ; 0.01 )
V3 ( -0.02 ; 0 ; 0 )
unit = 1 AU per day
and masses m1 = 2
; m2 = 1 ; m3 = 3
unit = 1 Solar mass
-Calculate the positions and velocities 10 days later with stepsize
= h = 10 , number of steps = 1
XEQ "GNBP" >>>> " M1=?"
2 R/S
" M2=?"
1 R/S
" M3=?"
3 R/S
" M4=?"
since there are only 3 stars, press R/S without any digit entry
R/S " X1^Y1^Z1"
2 ENTER^
0 ENTER^
0 R/S
"VX1^VY1^VZ1"
0 ENTER^
0.03 ENTER^
0 R/S
" X2^Y2^Z2"
0 ENTER^
4 ENTER^
0 R/S
"VX2^VY2^VZ2"
0 ENTER^
0 ENTER^
0.01 R/S
" X3^Y3^Z3"
0 ENTER^
0 ENTER^
1 R/S
"VX3^VY3^VZ3"
-0.02 ENTER^
0 ENTER^
0 R/S
" H^Nb?"
Stepsize ENTER Number of steps
10 ENTER^
1 R/S
>>>> 30.047 = R01
( Turbo mode quite useful ! )
>>> The new positions and velocities have replaced the previous
ones in registers R30 thru R47 and we get:
| Register | input | h = 10 ; N=1 | h = 5 ; N=2 | |
| m | R27=m1
R28=m2 R29=m3 |
2
1 3 |
undisturbed | undisturbed |
|
p
|
R30=x1
R31=y1 R32=z1 -------- R33=x2 R34=y2 R36=z2 -------- R36=x3 R37=y3 R38=z3 |
2
0 0 ----- 0 4 0 ----- 0 0 1 |
1.992077590
0.300333861 0.003673761 -------------- 0.000661665 3.996080594 0.100603408 -------------- -0.194938948 0.001083895 0.997349690 |
1.992077585
0.300333570 0.003673682 -------------- 0.000661669 3.996080575 0.100603412 -------------- -0.194938946 0.001084095 0.997349741 |
| v
e l . |
R39=x'1
R40=y'1 R41=z'1 -------- R42=x'2 R43=y'2 R44=z'2 -------- R45=x'3 R46=y'3 R47=z'3 |
0
0.03 0 ----- 0 0 0.01 ----- -0.02 0 0 |
-0.001550090
0.030038159 0.000706688 -------------- 0.000132598 -0.000790384 0.010117548 -------------- -0.019010806 0.000238022 -0.000510308 |
-0.001550083
0.030038158 0.000706684 -------------- 0.000132598 -0.000790385 0.010117549 -------------- -0.019010811 0.000238023 -0.000510306 |
-To estimate the accuracy of the results, we can perform the calculation
with a half step size ( h = 5 days ) and N = 2 instead of 1:
We obtain the results in the last column: errors are smaller
than 0.000001 AU
- To continue, simply press R/S the HP-41 displays " H^Nb?"
Press R/S without any digit entry if you keep the same stepsize & number of steps.
Notes:
-If you get NONEXISTENT during the initialization, increase the SIZE
( at least SIZE 27+19n = 084 in the example above )
-If you apply "GNBP" or "GNB+1" to the whole Solar system, h
= 1day gives an error of about 7 10-6 AU in the position
of Mercury after 88 days.
-The relativitic effects are not taken into account.
Gravitational (n+1)-body Problem
When computing planetary trajectories for instance, it's often
better to know directly the heliocentric coordinates by choosing the Sun
as the origin.
In this case however, the reference frame is no more inertial
, the equations become more complex and therefore, the program is longer.
On the other hand, it executes faster and requires less data
registers ( the 16-body problem can be solved ). The equations are now:
d2xi /dt2 = - Gm0 xi/ri3
- SUMj Gmj xj/rj3
+ SUMj#i Gmj ( xj -
xi ) / rij3
d2yi /dt2 = - Gm0 yi/ri3
- SUMj Gmj yj/rj3
+ SUMj#i Gmj ( yj -
yi ) / rij3
i , j = 1 , ..... , n
d2zi /dt2 = - Gm0 zi/ri3
- SUMj Gmj zj/rj3
+ SUMj#i Gmj ( zj -
zi ) / rij3
where ri = ( xi2 + yi2
+ zi2 ) 1/2 and
rij = [ ( xi - xj )2 + ( yi
- yj )2 + ( zi - zj )2
] 1/2
Data Registers: R00 thru R30+19n are used
R00 = n = Nb of bodies - 1
R31+n = x1
R31+4n = x'1
R27 = G = constant of gravitation
R32+n = y1
R32+4n = y'1
R28 = h = step size
R33+n = z1
R33+4n = z'1
R29 = N = number of steps
................
...................
R30 = m0 = mass of the origin
R31 = m1
R28+4n = xn
R28+7n = x'n
...............
R29+4n = yn
R29+7n = y'n
R30+4n = zn
R30+7n = z'n
R30+n = mn
Flag: F10
Subroutines: /
>>> Thus, the n masses, the 3n position coordinates and the 3n velocity
coordinates are to be stored into contiguous registers, starting at R30
| STACK | INPUT | OUTPUT |
| X | / | bbb.eee |
where bbb.eee = control number of the solution.
Example: let's take the same 3 stars and let's choose the third star as the origin. The coordinates become:
M1 ( 2 ; 0 ; -1 ) M2 ( 0 ; 4 ; -1 ) and V1 ( 0.02 ; 0.03 ; 0 ) V2 ( 0.02 ; 0 ; 0.01 )
-Calculate the positions and velocities 10 days later with stepsize = h = 10 , number of steps = 1
XEQ "GNB+1" >>>> " M0=?"
3 R/S
" M1=?"
2 R/S
" M2=?"
1 R/S
" M3=?"
since there are only 3 stars, press R/S without any digit entry
R/S " X1^Y1^Z1"
2 ENTER^
0 ENTER^
-1 R/S
"VX1^VY1^VZ1"
0.02 ENTER^
0.03 ENTER^
0 R/S
" X2^Y2^Z2"
0 ENTER^
4 ENTER^
-1 R/S
"VX2^VY2^VZ2"
0.02 ENTER^
0 ENTER^
0.01 R/S
" H^Nb?"
Stepsize ENTER Number of steps
10 ENTER^
1 R/S
>>>> 33.044 = R01
( Turbo mode quite useful again ! )
>>> the new postions and velocities have replaced the previous
ones in registers R33 thru R44 ( last column )
| Register | input | h = 10 ; N=1 | |
| m | R30=m0
R31=m1 R32=m2 |
3
2 1 |
undisturbed |
| p
o s. |
R33=x1
R34=y1 R35=z1 -------- R36=x2 R37=y2 R38=z2 |
2
0 -1 ----- 0 4 -1 |
2.187016538
0.299249966 -0.993675929 -------------- 0.195600614 3.994996700 -0.896746283 |
| v
e l. |
R39=x'1
R40=y'1 R41=z'1 -------- R42=x'2 R43=y'2 R44=z'2 |
0.02
0.03 0 ----- 0.02 0 0.01 |
0.017460717
0.029800137 0.001216996 -------------- 0.019143404 -0.001028406 0.010627856 |
Notes:
-If you get NONEXISTENT during the initialization, increase the SIZE
( at least SIZE 31+19n = 069 in the example above )
-If you apply "GNBP" or "GNB+1" to the whole Solar system, h
= 1day gives an error of about 7 10-6 AU in the position
of Mercury after 88 days.
-The relativitic effects are not taken into account.
Execution time:
-Here are the results of a few tests ( for 1 step ) and a real
HP-41
| n | GNBP | GNB+1 |
| 3 | 1mn52 | 1mn28 |
| 5 | 5mn04 | 4mn27 |
| 10 | 19mn50 | 18mn51 |
-Since 10 minutes become 1 second with V41, we can integrate the motion
of the 9 major planets in a reasonable time !
Illuminated Fraction of the Moon's Disk
-Given a date and a time expressed in TT, this short program returns
the illuminated fraction of the Moon at that instant.
-The results - though not extremely accurate - are valid over the time
span [1000,3000] at least.
Data Registers: R00 = T in millennia since
2000/01/01 0h
R01 = D , R02 = M' then M' - 2D
Flags: /
Subroutine: "J2000"
| STACK | INPUTS | OUTPUTS |
| Y | YYYY.MNDD | / |
| X | HH.MNSS | f |
where f = illuminated fraction of the Moon's disk ( 0 < f < 1 )
Examples: With CF 04 - Gregorian Calendar.
2010.0716 ENTER^
7.41
XEQ "IFMD" >>>> f = 0.2688 = 27%
3000.0404 ENTER^
0
R/S
>>>> f = 0.4942 = 49%
1001.0714 ENTER^
8.3916
R/S
>>>> f = 1.0000 = 100% ( almost
)
Hindu Calendar > Julian/Gregorian Calendar
-See "DT-KY" for the definitions.
Data Registers: R00-R01-R02: temp
Flag: F01-F04 CF 01 =
traditional old Hindu Calendar CF 04 = Gregorian
Calendar
SF 01 = "new-old" Hindu Calendar
SF 04 = Julian Calendar
Subroutine:
"DT"
| STACK | INPUT | OUTPUT |
| X | YYYY.MNDDKy | YYYY.MNDDGr/Ju |
Examples: With CF 04 ( Gregorian Calendar )
• CF 01 Traditional Hindu Calendar
0.0101
XEQ "KY-DT" >>>> -3101.0123
5080.0331 XEQ "KY-DT" >>>>
1979.0716
5110.1219 XEQ "KY-DT" >>>>
2010.0404
• SF 01 "New-Old" Hindu Calendar
0.0101
XEQ "KY-DT" >>>> -3101.0123
5080.0412 XEQ "KY-DT" >>>>
1979.0716
5111.0101 XEQ "KY-DT" >>>>
2010.0404
-"DT-KY" & "KY-DT" are just separated by a RTN
-This program calculates the local solar time defined by: LST
= 12 h. + the hour angle of the Sun.
-12h are added because it's more convenient to regard the solar day
as commencing at midnight rather than at noon.
-The equation of time ( i-e the difference between apparent and mean
time ) is obtained by the following formula:
E = y.sin 2L - 2e.sin M + 4ey.sin M cos 2L - (y2/2).sin 4L - (5e2/4).sin 2M ( in radians )
where
y = tan2(obl/2)
obl = obliquity of the ecliptic
L = Sun's mean longitude
e = eccentricity of the Earth's orbit
M = Sun's mean anomaly
Data Registers: R00 thru R05
Flag: /
Subroutine: "J2000"
| STACK | INPUTS | OUTPUTS |
| Z | YYYY.MNDD | / |
| Y | HH.MNSS (UT) | / |
| X | Longitude ( ° . ' " ) | LST = HH.MNSS |
Example: In the U.S. Naval Observatory at Washington ( D.C.) ( Longitude = 77°03'56"W ) on 2010/07/16 at 16h24m41s UT
2010.0716 ENTER^
16.2441
ENTER^
-77.0356 XEQ
"LST" >>>> LST = 11h10mn21s
Notes:
-Longitudes are measured positively westward from the meridian of Greenwich
-The accuracy is of the order of 1 or 2 seconds between A.D. 1000 &
3000
Passage of the Moon through the Nodes
-This program takes an approximate date and returns the time when the
Moon passes at the ascending or the descending node of its orbit ( Moon's
latitude = 0 ).
-The accuracy is about 5 minutes over the interval [1900,2100]
Data Registers: R00 thru R04: temp
( R04 is an even integer for the ascending node or an odd integer for the
descending node )
Flags: /
Subroutines: "J2000" & "DT"
| STACK | INPUTS | OUTPUTS |
| Z | / | HH.MNSS |
| Y | / | YYYY.MNDD |
| X | yyyy.mndd | +/-1 |
X-output = +1 for a passage at the ascending
node
X-output = -1 for a passage at the descending
node
Examples:
2008.0825 XEQ "MNODE" >>>>
1
RDN 2008.0912
RDN 18.2530
-The Moon passes through the ascending node on 2008/09/12 at 18h25m TT
2008.0915
R/S
>>>> -1
RDN 2008.0925
RDN 16.2137
-The Moon passes through the descending node on 2008/09/25 at 16h22m
TT
Nutation ( simplified formulae )
-The complete series involve more than 1000 terms, and the following
program uses the first few terms only:
-The nutation in longitude Dpsi is obtained with an accuracy of
0"5
and the nutation in obliquity De -----------------------------
0"1
Data Register:
• R00 = T = number of millenia since 2000/01/01 0h
( this register is to be initialized before executing "NUT" )
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Y | / | De ( deg ) |
| X | / | Dpsi ( deg ) |
where Dpsi is the nutation in longitude
and De is the nutation
in obliquity.
Example: With T = 1 ( 3000/01/08 0h )
1 STO 00 XEQ "NUT" >>>> Dpsi
= 0.00353°
X<>Y De = -0.00198°
"OBL" uses the series given in reference [5] to compute the mean
obliquity of the ecliptic em
Data Register:
• R00 = T = number of millenia since 2000/01/01 0h
( this register is to be initialized before executing "OBL" )
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Y | / | T |
| X | / | em ( °. ' " ) |
where em is the mean obliquity
Example: With T = 1 ( 3000/01/08 0h )
1 STO 00 XEQ "OBL" >>>> em = 23°18'35"02
-The accuracy is of the order of 0"01 over the time-span [1000,3000]
-The phases of the Moon are usually defined by the difference between
the longitudes of the Moon and the Sun ( 0° 90° 180° 270°
)
-In the following routine, we use this definition if flag F00 is set,
but if flag F00 is clear,
we use another definition i-e when the illuminated fraction
of the Moon's disk equals 0 , 1/2 , 1 , 1/2
-Thus - unlike the standard definition - the instant of a Solar eclipse
is identical to the instant of the New Moon
( likewise for a Lunar eclipse and a Full Moon ) which
is perhaps more logical...
-This program is actually combined with "ECLIPSE"
-Accuracy is about 5 minutes over the time-span 1000-3000.
Data Registers: R00 thru R16: temp
R04 = phase
Flag: F00
SF 00 = Standard definition of phases
CF 00 = Non-Standard definition of phases
Subroutines: "J2000" & "DT"
| STACK | INPUTS | OUTPUTS |
| Z | / | HH.MNSS |
| Y | / | YYYY.MNDD |
| X | yyyy.mndd | Phase |
where yyyy.mndd is an approximate date
Phase = 0 for a New Moon
Phase = 1 for a First Quarter
Phase = 2 for a Full Moon
Phase = 3 for a Last Quarter
YYYY.MNDD & HH.MNSS are the date & time ( in Terrestrial Time ) of the Phase
Example: CF 00 SF 00
2008.0801 XEQ "PHASE" >>>>
0
0
RDN 2008.0801
2008.0801
RDN 10.2426
10.1618
So, there is a New Moon on 2008/08/01 at 10h24m26s TT 10h16m18s TT with the standard definition
-Likewise, you find a First Quarter on 2008/08/08 at
20h02m45s TT
20h21m20s TT ---------------------------
Full Moon on 2008/08/16 at 21h12m40s
TT 21h19m46s
TT ---------------------------
Last Quarter on 2008/08/24 at 00h08m41s TT
23h52m48s TT on 2008/08/23 ----------
Note:
-LBL "PHASE" is the 1st line of the program
-LBL "ECLIPSE" is line 132
-"RTS" calculates the times of rising, transit and setting of a celestial
body on a day D
and the geometric altitude at the moment of transit.
-The coordinates ( right ascension and declination ) must be obtained
independently on days D-1, D , D+1 ( 0h UT )
Data Registers: R00 and R07 thru R15 are used for temporary data storage. ( Registers R01 to R06 are to be initialized before executing "RTS" )
( on day D-1 ) ( on day D ) ( on day D+1 ) ( at 0h )
• R01 = RA(-1) • R03 = RA(0) •
R05 = RA(+1) RA =
right ascension ( hh.mnss )
R11 = Longitude
• R02 = decl(-1) • R04 = decl(0) • R06 =
decl(+1) decl = declination
( ° . ' " )
R12 = Latitude
Flags: Set Flag F00
for the Sun
Set Flag F02 for the Moon ( Clear these
2 flags for the other celestial bodies )
Subroutine: "MST"
| STACK | INPUTS | OUTPUTS |
| T | / | h = altitude* ( ° . ' " ) |
| Z | Date ( YYYY.MNDD ) | setting(UT) (hh.mnss) |
| Y | Longitude ( ° . ' " ) | transit(UT) (hh.mnss) |
| X | latitude ( ° . ' " ) | rising (UT) (hh.mnss) |
* where h = the geometric altitude of the center of the body at the time of transit.
Example1: Find the times of rising, transit and setting of the Sun at Seattle ( 122°19'51"W ; 47°36'23"N ) on 2005/01/27
-First, we find the right ascension and declination of the Sun on 2005/01/26
, 2005/01/27 , 2005/01/28 at 0h UT ( cf "Astronomical Ephemeris for
the HP-41" )
and we store these values into R01 thru R06:
20.3404 STO 01
20.3813 STO 03 20.4221
STO 05
-18.4417 STO 02
-18.2900 STO 04 -18.1324 STO
06
then, SF 00 CF 02
2005.0127
ENTER^
-122.1951 ENTER^
47.3623 XEQ "RTS" >>>> Rising
= 15h41m39s (UT)
RDN Transit = 20h22m12s (UT)
RDN Setting = 25h03m17s (UT) ( i-e
1h03m17s on 2005/01/28 )
RDN altitude of the Sun
at the time of transit = 24°07'48"
-For the Moon, CF 00 SF 02
Example2: Same questions for Sirius. The right ascension and declination are 6h45m23s and -16°43'18" for this date.
6.4523
STO 01 STO 03 STO 05
the coordinates may be regarded as constant.
-16.4318 STO 02
STO 04 STO 06
then, CF 00 CF 02
2005.0127 ENTER^ ( or
2005.0127 RCL 11 RCL 12 XEQ
"RTS" )
-122.1951 ENTER^
47.3623 R/S
>>>> Rising = 1h42m05s (UT)
RDN Transit = 6h28m09s (UT)
RDN Setting = 11h14m14s (UT)
RDN
altitude of Sirius at the time of transit = 25°40'19"
Notes:
-If the body is circumpolar, there will be no rising and setting: the
star will remain above ( or below ) the horizon.
-In this case however, "RTS" yields ( approximately ) the same values
in X, Y, Z registers
-For example, with the pole star ( RA = 2h37m39s ; Decl = 89°17'10"
) "RTS" gives X = Y = Z = 2h21m06s
but since T = altitude of the star at the instant of transit
= 48°19'13"
we can deduce that X and Z are meaningless results.
Sidereal & Tropical years and Lunar Months
-"STYLM" computes the mean sidereal year, tropical year and lunar
( synodic ) month for a given instant.
-The following formulas are used:
1 sidereal year = 365.2563630 + 0.0000115
T + 0.0000003 T2
1 tropical year = 365.2421904 - 0.0006153
T - 0.0000061 T2 + 0.0002652 T3 + 0.0000254
T4 - 0.0000338 T5
1 lunar month = 29.53058886 +
0.0002496 T - 0.0000364 T2 + 0.0000235
T3
where T is measured in unit of 10000 julian years from J2000.0
Data Registers: R00 = sidereal year
R01 = tropical year
Flags: /
Subroutines: /
| STACK | INPUTS | OUTPUTS |
| Z | / | lunar month |
| Y | / | tropical year |
| X | year | sidereal year |
All the results are expressed in days
Examples:
2000 XEQ "STYLM" >>>> 365.2563630
days
RDN 365.2421904 days
RDN 29.53058886 days
5000 R/S
>>>> 365.2563665
RDN 365.2420125
RDN 29.53059608
12000 R/S
>>>> 365.2563748
RDN 365.2418258
RDN 29.53061253
-3000 R/S
>>>> 365.2563573
RDN 365.2424660
RDN 29.53057518
-8000 R/S
>>>> 365.2563518
RDN 365.2425936
RDN 29.53055791
Notes:
-The input may be fractional.
-The outputs are less accurate before -8000 and after 12000 especially
for the lunar months,
because the tidal acceleration of the Moon is not known with
a great precision.
Date & Time for a given longitude of the Sun
-"ZODIAC" calculates the instant of a year when the longitude of the
Sun equals a given value.
-The nutation & aberration are taken into account.
Data Register: R00 & R01: temp
Flag: F01
Subroutines: "SUN" "DT"
| STACK | INPUTS | OUTPUTS |
| Y | YYYY | hh.mnss (TT) |
| X | L | YYYY.MNDD |
where YYYY = year and L = longitude is expressed in decimal degrees
Examples: For the same year 2010 , it yields
2010 ENTER^
0 XEQ
"ZODIAC" >>>> 2010.0320
X<>Y 17.3202
which is the instant of the spring equinox
2010 ENTER^
90
R/S >>>>
2010.0621
X<>Y 11.2910
which is the instant of the summer solstice
2010 ENTER^
180 R/S
>>>> 2010.0923
X<>Y 3.1257
which is the instant of the autumn equinox
2010 ENTER^
270 R/S
>>>> 2010.1221
X<>Y 23.4055
which is the instant of the winter solstice
-Of course, L is not necessarily a multiple of 90°, for instance:
2010 ENTER^
234 R/S
>>>> 2010.1116
X<>Y 11.3131
-However, for an instant before the spring equinox, subtract
360° to L.
-Otherwise, the result would concern the following year !
References:
[1] Jean Meeus - "Astronomical Algorithms" - Willmann-Bell
- ISBN 0-943396-61-1
[2] Robin M. Green - "Spherical Astronomy" - Cambridge University
Press - ISBN 0-521-31779-7
[3] John D. Anderson - JPL Technical Report n° 32-497
[4] Dan Boulet - "Methods of Orbit Determination for the Micro
Computer" - Willmann-Bell - ISBN 0-943396-34-4
[5] United States Naval Observatory, Circular n° 179 http://aa.usno.navy.mil/publications/docs/circular_179.html
[6] Nachum Dershowitz & Edward M. Reingold - "Calendrical
Calculations" - Cambridge University Press - ISBN 978-0-521-70238-6
[7] Kantowski, Kao, Thomas - "Distance-Redshift Relations in
Inhomogeneous Friedmann-Lemaitre-Robertson-Walker Cosmology" -
The Astrophysical Journal, 545: 549-560