Wellcome to the 6th Port of Call on our journey to Slovenian Nautical Almanac. Things we've covered so far:
- Explained what and why, we're doing it
- Described the JPL DE 405 planetary Ephemerides files and calculation of J2000 Rectangular Coordinates based on Sloar System Barycenter
- Calculated Geocentric Equatorial Astrometric coordinates from XYZ values, calculated in 2
- Learned how to get proper Ephemerides time from the basic Universal Coordinated Time, i.e. our measure of time on Earth
Now we need to obtain apparent positions of the Sun, Moon ands the planets.
As neither the ecliptic pole, nor celestial pole is fixed in space in relation to distant stars, and we have calculated the positions based on ecliptic and celestial pole at J2000, we need to apply three corrections to Astrometric RA and Declination:
1. Precession
The biggest effect on celestial pole movement is caused by the combined effect of gravitational pull by the Sun and the Moon, and to lesser degree, planets. This causes Celestial pole to "sweep" a cone around the pole of Ecliptic, with a half-angle of about 23,5°. Complete base circle takes some 26000 years to complete. And since the position of Celestial pole is constantly changing over time, the position of Vernal Equinox moves as well.
To correct this I followed an approach by Jan Meeus in his Astronomical Algorithms book, published by Willman-Bell. It was this book, which I was able to obtain second hand on a local book sell, that made me think I can do all this in the first place. It is really a fantasic book on the subject. Not to mention the fact that it devotes a lot of pages to basic arithetics and math in general as well.
1.1. Coarse correction
The following formulae may be used, if no great accuracy is needed and if the dates between observation and standard epoch are not very far apart:
ΔRA = ((m + n * sin(RA) * tan(Dec) * y
ΔDec = (n * cos(RA)) * y
where m and n are not really constants, but van be calculated by the following equations:
m = 46,1244" + 0,0279" * T
n = 20,04315" + 0,00855" * T
where T is the number of centuries between J2000 (JD = 2451545) and The date in question. y is the number of years between J2000.
The code, which solves this is shown bellow:
Public Class Precession Public Function CoarseCorrection(ByVal Jultime As Decimal, ByVal EC As Equatorial) As Equatorial 'Centuries since J2000 Dim T As Decimal = (Jultime - 2451545) / 36525 'Years since J2000 Dim y As Decimal = T * 100 Dim m As Decimal = (46.1244 + 0.0279 * T) * Constants.dtr / 3600 Dim n As Decimal = (20.04315 - 0.00855 * T) * Constants.dtr / 3600 'Results Dim Corr As New Equatorial Corr.RightAscension = (m + n * Math.Sin(EC.RightAscension) * Math.Tan(EC.Declination)) * y Corr.Declination = (n * Math.Cos(EC.RightAscension)) * y Return Corr End Function End Class
2.2. Fine Correction
When a great accuracy is required, the rgirous procedure is used, again based on Meuus book:
Declare T as number of Julian Centuries between observation Date and J2000:
T = (JDE - 2451545) / 36525
Then we define three aux variables: ζ, z, and θ (In degrees decimal):
ζ = (2306.2181 * T + 0.30188 * T ^ 2 + 0.017998 * T ^ 3) / 3600
z = (2306.2181 * T + 1.09468 * T ^ 2 + 0.018203 * T ^ 3) / 3600
θ = (2004.3109 * T + 0.42665 * T ^ 2 + 0.41833 * T ^ 3) / 3600
Then we define three aux variables, A, B, and C:
A = Cos(Dec) * Sin(RA + ζ)
B = Cos(θ) * Cos(Dec) * Cos(RA + ζ) - Sin(θ) * Sin(Dec)
C = Sin(θ) * Cos(Dec) * Cos(RA + ζ) - Cos(θ) * Sin(Dec)
We finally calculate corrected Right Ascension and Declination by:
RA(c) = Atan2(A, B) + Z
Dec(c) = Asin(C)
If celestial body is close to either North or South Celestial Pole (Say, that it's Declination is above N 85°, or below S85°), then we use the following formula for Declination:
Dec(c) = Acos(Sqrt(A ^ 2 + B ^ 2))
To obtain corrections, i.e., ΔRA, and ΔDec:
ΔRA = RA(c) - RA
ΔDec = Dec(c) - Dec
The code for a rigorous method, is shown bellow:
''' <summary> ''' Calculates Fine precession correction in Right Ascension and Declination ''' as given in Meeus: Astronomical Algoritms, chapter 13, page ''' </summary> ''' <param name="Jultime">Julian day</param> ''' <param name="EC">Input Equatorial Coordinates</param> ''' <returns>EquatorialCoordinates Object containg the corrections for RA and Declination (In Radians)</returns> Public Function FineCorrection(ByVal Jultime As Decimal, ByVal EC As Equatorial) As Equatorial Dim Corr As New Equatorial Dim T As Decimal = (Jultime - 2451545) / 36525 Dim Zeta As Decimal = Constants.dtr * (2306.2181 * T + 0.30188 * T ^ 2 + 0.017998 * T ^ 3) / 3600 Dim Z As Decimal = Constants.dtr * (2306.2181 * T + 1.09468 * T ^ 2 + 0.018203 * T ^ 3) / 3600 Dim Theta As Decimal = Constants.dtr * (2004.3109 * T + 0.42665 * T ^ 2 + 0.41833 * T ^ 3) / 3600 Dim A As Decimal = Math.Cos(EC.Declination) * Math.Sin(EC.RightAscension + Zeta) Dim B As Decimal = Math.Cos(Theta) * Math.Cos(EC.Declination) * Math.Cos(EC.RightAscension + Zeta) - Math.Sin(Theta) * Math.Sin(EC.Declination) Dim C As Decimal = Math.Sin(Theta) * Math.Cos(EC.Declination) * Math.Cos(EC.RightAscension + Zeta) + Math.Cos(Theta) * Math.Sin(EC.Declination) Dim Corr_EC As New Equatorial Dim RA As Decimal RA = Math.Atan2(A, B) + Z RA = MathFunctions.TruncateRad(RA) Corr.RightAscension = RA - EC.RightAscension Dim Dec As Decimal = Math.Asin(C) Corr.Declination = Dec - EC.Declination Return Corr End Function
2. Nutation
Nutation can be seen as a first order correction to Precession, described above. It is caused by gravitational pull of the moon on the earth. The Earth's crust is not uniform, and massive amounts of water in the oceans is not distributed evenly around the globe. This causes slight shift in position of Earth's axis of rotation, with a period of approximately 18,6 years.
I had found a beautifull picture to describe both phenomena on Wikipedia: By User Herbye (German Wikipedia). Designed by Dr. H. Sulzer - Original, CC BY-SA 3.0
To find the required corrections for RA and Dec, we follow the IAU 1980 theory. First we again define T as number of Julian Centuries between Date and J2000:
T = (JDE - 2451545) / 36525
Then we calculate Mean obliquity of the Ecliptic, that is mean angle between Earth's equator and the plane of Ecliptic:
ε0 = (84381.448 - 46.815 * T - 0.00059 * T ^ 2 + 0.001813 * T ^ 3) / 3600
To get True Obliquity of Ecliptic ε one needs to apply the correction called Δε (DeltaEps). This quantity, together with correction for Nutaton in Longitude, called Δψ (DeltaPsi) can be calculated as follows:
First compute main Delaunay arguments (5 in total):
- Mean anomaly of the Moon L:
L = (485868.249036 + 1717915923.2178 * T + 31.8792 * T^2 + 0.051635 * T^3 - 0.00024470 * T^4) / 3600.0 - Mean anomaly of the Sun Lp:
Lp = (1287104.79305 + 129596581.0481 * T - 0.5532 * T^2 + 0.000136 * T^3 - 0.00001149 * T^4) / 3600.0 - Mean argument of the latitude of the Moon F
F = (335779.526232 + 1739527262.8478 * T - 12.7512 * T^2 - 0.001037 * T^3 + 0.00000417 * T^4) / 3600.0 - Mean elongation of the Moon from the Sun D:
D = (1072260.70369 + 1602961601.2090 * T - 6.3706 * T^2 + 0.006593 * T^3 - 0.00003169 * T^4) / 3600.0 - Mean longitude of the ascending node of the Moon Om:
Om = (450160.398036 - 6962890.5431 * T + 7.4722 * T^2 + 0.007702 * T^3 - 0.00005939 * T^4) / 3600.0
To calculate Δε and Δψ we use IAU2000B Data, which I found on Jay Tanner website, which also provides the computation method.
There are 77 sets of data. Each consists of:
- 5 Delaunay argument multipliers (M1 to M5)
- 3 Longitude coefficients (Lon1 to Lon3), measured in 10000000 of Arc seconds
- 3 Obliquity coefficients (Obl1 to Obl), measured in 10000000 of Arc seconds
For each set of data we first compute the full Argument (Arg):
Arg = m1 * L + m2 * Lp + m3 * F + m4 * D + m5 * Om
Then we compute each row Tems for Δε and Δψ as:
dPsiTerm = (Lon1 + Lon2 * T) * Sin(Arg) + Lon3 * Cos(Arg)
dEpsTerm = (Obl1 + Obl2 * T) * Cos(Arg) + Obl3 * Sin(Arg)
Summing up these terms and dividing the sum by 36000000000 will give us Δε and Δψ in decimal degrees.
Δψ = SUM(dPsiTerm) / 36000000000
Δε = SUM(dEpsTerm) / 36000000000
And the True obliquity of the ecliptic is then:
ε = ε0 + Δε
To correct RA and Declination for nutation, we follow the following code:
ΔRA = (Cos(ε) + Sin(ε) * Sin(RA) * Tan(Dec)) * Δψ - (Cos(RA) * Tan(Dec)) * Δε
ΔDec = (Sin(ε) * Cos(RA)) * Δψ + Sin(RA) * Δε
To put this into code, we first declare a class Called Nutation:
Public Class Nutation Public Property MeanObliquity As Decimal Public Property TrueObliquity As Decimal Public Property dEps As Decimal Public Property dPsi As Decimal End Class
IAU200BNutationItem as a POCO, to hold the IAU 2000B Nutation terms as described above:Public Class IAU2000BNutationItem Public Property m1 As Integer Public Property m2 As Integer Public Property m3 As Integer Public Property m4 As Integer Public Property m5 As Integer Public Property Lon1 As Decimal Public Property Lon2 As Decimal Public Property Lon3 As Decimal Public Property Obl1 As Decimal Public Property Obl2 As Decimal Public Property Obl3 As Decimal ''' <summary> ''' Constructor with parameters ''' </summary> Public Sub New(ByVal m1 As Integer, ByVal m2 As Integer, ByVal m3 As Integer, ByVal m4 As Integer, ByVal m5 As Integer, ByVal Lon1 As Decimal, ByVal Lon2 As Decimal, ByVal Lon3 As Decimal, ByVal Obl1 As Decimal, ByVal Obl2 As Decimal, ByVal Obl3 As Decimal) With Me .m1 = m1 .m2 = m2 .m3 = m3 .m4 = m4 .m5 = m5 .Lon1 = Lon1 .Lon2 = Lon2 .Lon3 = Lon3 .Obl1 = Obl1 .Obl2 = Obl2 .Obl3 = Obl3 End With End Sub End Class
IAU2000NutationItem objects, we create a class, called
IAU2000BNutationList, which inherits from
List(of T):
''' <summary> ''' Represent the list of IAU2000BNutationItem Items, needed to compute the Nutation parameters ''' </summary> Public Class IAU2000BNutationList Inherits List(Of IAU2000BNutationItem) ''' <summary> ''' Default Constructor ''' </summary> Public Sub New() 'If any data exitst inside, remove it: Me.Clear() 'Add the data to the list: Me.Add(New IAU2000BNutationItem(0, 0, 0, 0, 1, -172064161, -174666, 33386, 92052331, 9086, 15377)) Me.Add(New IAU2000BNutationItem(0, 0, 2, -2, 2, -13170906, -1675, 13696, 5730336, -3015, -4587)) Me.Add(New IAU2000BNutationItem(0, 0, 2, 0, 2, -2276413, -234, 2796, 978459, -485, 1374)) Me.Add(New IAU2000BNutationItem(0, 0, 0, 0, 2, 2074554, 207, -698, -897492, 470, -291)) Me.Add(New IAU2000BNutationItem(0, 1, 0, 0, 0, 1475877, -3633, 11817, 73871, -184, -1924)) Me.Add(New IAU2000BNutationItem(0, 1, 2, -2, 2, -516821, 1226, -524, 224386, -677, -174)) Me.Add(New IAU2000BNutationItem(1, 0, 0, 0, 0, 711159, 73, -872, -6750, 0, 358)) Me.Add(New IAU2000BNutationItem(0, 0, 2, 0, 1, -387298, -367, 380, 200728, 18, 318)) Me.Add(New IAU2000BNutationItem(1, 0, 2, 0, 2, -301461, -36, 816, 129025, -63, 367)) Me.Add(New IAU2000BNutationItem(0, -1, 2, -2, 2, 215829, -494, 111, -95929, 299, 132)) Me.Add(New IAU2000BNutationItem(0, 0, 2, -2, 1, 128227, 137, 181, -68982, -9, 39)) Me.Add(New IAU2000BNutationItem(-1, 0, 2, 0, 2, 123457, 11, 19, -53311, 32, -4)) Me.Add(New IAU2000BNutationItem(-1, 0, 0, 2, 0, 156994, 10, -168, -1235, 0, 82)) Me.Add(New IAU2000BNutationItem(1, 0, 0, 0, 1, 63110, 63, 27, -33228, 0, -9)) Me.Add(New IAU2000BNutationItem(-1, 0, 0, 0, 1, -57976, -63, -189, 31429, 0, -75)) Me.Add(New IAU2000BNutationItem(-1, 0, 2, 2, 2, -59641, -11, 149, 25543, -11, 66)) Me.Add(New IAU2000BNutationItem(1, 0, 2, 0, 1, -51613, -42, 129, 26366, 0, 78)) Me.Add(New IAU2000BNutationItem(-2, 0, 2, 0, 1, 45893, 50, 31, -24236, -10, 20)) Me.Add(New IAU2000BNutationItem(0, 0, 0, 2, 0, 63384, 11, -150, -1220, 0, 29)) Me.Add(New IAU2000BNutationItem(0, 0, 2, 2, 2, -38571, -1, 158, 16452, -11, 68)) Me.Add(New IAU2000BNutationItem(0, -2, 2, -2, 2, 32481, 0, 0, -13870, 0, 0)) Me.Add(New IAU2000BNutationItem(-2, 0, 0, 2, 0, -47722, 0, -18, 477, 0, -25)) Me.Add(New IAU2000BNutationItem(2, 0, 2, 0, 2, -31046, -1, 131, 13238, -11, 59)) Me.Add(New IAU2000BNutationItem(1, 0, 2, -2, 2, 28593, 0, -1, -12338, 10, -3)) Me.Add(New IAU2000BNutationItem(-1, 0, 2, 0, 1, 20441, 21, 10, -10758, 0, -3)) Me.Add(New IAU2000BNutationItem(2, 0, 0, 0, 0, 29243, 0, -74, -609, 0, 13)) Me.Add(New IAU2000BNutationItem(0, 0, 2, 0, 0, 25887, 0, -66, -550, 0, 11)) Me.Add(New IAU2000BNutationItem(0, 1, 0, 0, 1, -14053, -25, 79, 8551, -2, -45)) Me.Add(New IAU2000BNutationItem(-1, 0, 0, 2, 1, 15164, 10, 11, -8001, 0, -1)) Me.Add(New IAU2000BNutationItem(0, 2, 2, -2, 2, -15794, 72, -16, 6850, -42, -5)) Me.Add(New IAU2000BNutationItem(0, 0, -2, 2, 0, 21783, 0, 13, -167, 0, 13)) Me.Add(New IAU2000BNutationItem(1, 0, 0, -2, 1, -12873, -10, -37, 6953, 0, -14)) Me.Add(New IAU2000BNutationItem(0, -1, 0, 0, 1, -12654, 11, 63, 6415, 0, 26)) Me.Add(New IAU2000BNutationItem(-1, 0, 2, 2, 1, -10204, 0, 25, 5222, 0, 15)) Me.Add(New IAU2000BNutationItem(0, 2, 0, 0, 0, 16707, -85, -10, 168, -1, 10)) Me.Add(New IAU2000BNutationItem(1, 0, 2, 2, 2, -7691, 0, 44, 3268, 0, 19)) Me.Add(New IAU2000BNutationItem(-2, 0, 2, 0, 0, -11024, 0, -14, 104, 0, 2)) Me.Add(New IAU2000BNutationItem(0, 1, 2, 0, 2, 7566, -21, -11, -3250, 0, -5)) Me.Add(New IAU2000BNutationItem(0, 0, 2, 2, 1, -6637, -11, 25, 3353, 0, 14)) Me.Add(New IAU2000BNutationItem(0, -1, 2, 0, 2, -7141, 21, 8, 3070, 0, 4)) Me.Add(New IAU2000BNutationItem(0, 0, 0, 2, 1, -6302, -11, 2, 3272, 0, 4)) Me.Add(New IAU2000BNutationItem(1, 0, 2, -2, 1, 5800, 10, 2, -3045, 0, -1)) Me.Add(New IAU2000BNutationItem(2, 0, 2, -2, 2, 6443, 0, -7, -2768, 0, -4)) Me.Add(New IAU2000BNutationItem(-2, 0, 0, 2, 1, -5774, -11, -15, 3041, 0, -5)) Me.Add(New IAU2000BNutationItem(2, 0, 2, 0, 1, -5350, 0, 21, 2695, 0, 12)) Me.Add(New IAU2000BNutationItem(0, -1, 2, -2, 1, -4752, -11, -3, 2719, 0, -3)) Me.Add(New IAU2000BNutationItem(0, 0, 0, -2, 1, -4940, -11, -21, 2720, 0, -9)) Me.Add(New IAU2000BNutationItem(-1, -1, 0, 2, 0, 7350, 0, -8, -51, 0, 4)) Me.Add(New IAU2000BNutationItem(2, 0, 0, -2, 1, 4065, 0, 6, -2206, 0, 1)) Me.Add(New IAU2000BNutationItem(1, 0, 0, 2, 0, 6579, 0, -24, -199, 0, 2)) Me.Add(New IAU2000BNutationItem(0, 1, 2, -2, 1, 3579, 0, 5, -1900, 0, 1)) Me.Add(New IAU2000BNutationItem(1, -1, 0, 0, 0, 4725, 0, -6, -41, 0, 3)) Me.Add(New IAU2000BNutationItem(-2, 0, 2, 0, 2, -3075, 0, -2, 1313, 0, -1)) Me.Add(New IAU2000BNutationItem(3, 0, 2, 0, 2, -2904, 0, 15, 1233, 0, 7)) Me.Add(New IAU2000BNutationItem(0, -1, 0, 2, 0, 4348, 0, -10, -81, 0, 2)) Me.Add(New IAU2000BNutationItem(1, -1, 2, 0, 2, -2878, 0, 8, 1232, 0, 4)) Me.Add(New IAU2000BNutationItem(0, 0, 0, 1, 0, -4230, 0, 5, -20, 0, -2)) Me.Add(New IAU2000BNutationItem(-1, -1, 2, 2, 2, -2819, 0, 7, 1207, 0, 3)) Me.Add(New IAU2000BNutationItem(-1, 0, 2, 0, 0, -4056, 0, 5, 40, 0, -2)) Me.Add(New IAU2000BNutationItem(0, -1, 2, 2, 2, -2647, 0, 11, 1129, 0, 5)) Me.Add(New IAU2000BNutationItem(-2, 0, 0, 0, 1, -2294, 0, -10, 1266, 0, -4)) Me.Add(New IAU2000BNutationItem(1, 1, 2, 0, 2, 2481, 0, -7, -1062, 0, -3)) Me.Add(New IAU2000BNutationItem(2, 0, 0, 0, 1, 2179, 0, -2, -1129, 0, -2)) Me.Add(New IAU2000BNutationItem(-1, 1, 0, 1, 0, 3276, 0, 1, -9, 0, 0)) Me.Add(New IAU2000BNutationItem(1, 1, 0, 0, 0, -3389, 0, 5, 35, 0, -2)) Me.Add(New IAU2000BNutationItem(1, 0, 2, 0, 0, 3339, 0, -13, -107, 0, 1)) Me.Add(New IAU2000BNutationItem(-1, 0, 2, -2, 1, -1987, 0, -6, 1073, 0, -2)) Me.Add(New IAU2000BNutationItem(1, 0, 0, 0, 2, -1981, 0, 0, 854, 0, 0)) Me.Add(New IAU2000BNutationItem(-1, 0, 0, 1, 0, 4026, 0, -353, -553, 0, -139)) Me.Add(New IAU2000BNutationItem(0, 0, 2, 1, 2, 1660, 0, -5, -710, 0, -2)) Me.Add(New IAU2000BNutationItem(-1, 0, 2, 4, 2, -1521, 0, 9, 647, 0, 4)) Me.Add(New IAU2000BNutationItem(-1, 1, 0, 1, 1, 1314, 0, 0, -700, 0, 0)) Me.Add(New IAU2000BNutationItem(0, -2, 2, -2, 1, -1283, 0, 0, 672, 0, 0)) Me.Add(New IAU2000BNutationItem(1, 0, 2, 2, 1, -1331, 0, 8, 663, 0, 4)) Me.Add(New IAU2000BNutationItem(-2, 0, 2, 2, 2, 1383, 0, -2, -594, 0, -2)) Me.Add(New IAU2000BNutationItem(-1, 0, 0, 0, 2, 1405, 0, 4, -610, 0, 2)) Me.Add(New IAU2000BNutationItem(1, 1, 2, -2, 2, 1290, 0, 0, -556, 0, 0)) End Sub End Class
IAU2000BNutationCalculation class, with two methods, GetNutation, which returns Nutation object for a given JDE, and NutationCorrection, which returns the corrections for Right Ascension and Declination for a given JDE and Equatorial Coordinates:''' <summary> ''' Class to represent Nutation parameter calculation and Equatorial coordinates corrections (in RA and Declination) ''' </summary> Public Class IAU2000BNutationCalculation Private IAUData As New IAU2000BNutationList ''' <summary> ''' Gets the Nutation object for a given JDE ''' </summary> ''' <param name="JDE">Julian Ephemeris date</param> ''' <returns>Nutation object, containing ε0, ε, Δψ, and Δε</returns> Public Function GetNutation(ByVal JDE As Decimal) As Nutation 'Output Nutation object Dim Result As New Nutation 'Julian Centuries since J2000: Dim T As Decimal = (JDE - 2451545) / 36525 'Compute Main Delaunay Arguments: 'Mean anomaly of the Moon in radians Dim L As Decimal = Constants.dtr * (485868.249036 + 1717915923.2178 * T + 31.8792 * T ^ 2 + 0.051635 * T ^ 3 - 0.0002447 * T ^ 4) / 3600.0 'Mean anomaly of the Sun in radians Dim Lp As Decimal = Constants.dtr * (1287104.79305 + 129596581.0481 * T - 0.5532 * T ^ 2 + 0.000136 * T ^ 3 - 0.00001149 * T ^ 4) / 3600.0 'Mean argument of the latitude of the Moon in radians Dim F As Decimal = Constants.dtr * (335779.526232 + 1739527262.8478 * T - 12.7512 * T ^ 2 - 0.001037 * T ^ 3 + 0.00000417 * T ^ 4) / 3600.0 'Mean elongation of the Moon from the Sun in radians Dim D As Decimal = Constants.dtr * (1072260.70369 + 1602961601.209 * T - 6.3706 * T ^ 2 + 0.006593 * T ^ 3 - 0.00003169 * T ^ 4) / 3600.0 'Mean longitude of the ascending node of the Moon in radians Dim Om As Decimal = Constants.dtr * (450160.398036 - 6962890.5431 * T + 7.4722 * T ^ 2 + 0.007702 * T ^ 3 - 0.00005939 * T ^ 4) / 3600.0 'Now loop through the IAUData and compute DeltaEps and DeltaPsi values 'Main Argument Dim Arg As Decimal Dim item As IAU2000BNutationItem 'Reset Δψ and Δε Values Result.dEps = 0 Result.dPsi = 0 For i As Integer = 0 To IAUData.Count - 1 'Get the IAU2000BNutationItem item = IAUData(i) 'Calculate the Main Argument Arg = item.m1 * L + item.m2 * Lp + item.m3 * F + item.m4 * D + item.m5 * Om 'Add individual data term to Δψ and Δε Result.dPsi += (item.Lon1 + item.Lon2 * T) * Math.Sin(Arg) + item.Lon3 * Math.Cos(Arg) Result.dEps += (item.Obl1 + item.Obl2 * T) * Math.Cos(Arg) + item.Obl3 * Math.Sin(Arg) Next 'Convert Δψ and Δε to radians from 10000000 of ArcSeconds Result.dEps = Constants.dtr * Result.dEps / 36000000000 Result.dPsi = Constants.dtr * Result.dPsi / 36000000000 'Mean Obliquity of the ecliptic in radians: Result.MeanObliquity = Constants.dtr * (84381.448 - 46.815 * T - 0.00059 * T ^ 2 + 0.001813 * T ^ 3) / 3600 'True obliquity of the Ecliptic in radians: Result.TrueObliquity = Result.MeanObliquity + Result.dEps Return Result End Function ''' <summary> ''' Calculates corrections in Right ascension and Declination for Nutation. ''' </summary> ''' <param name="EC">Equatorial coordinates to be corrected</param> ''' <param name="JDE">Julian Ephemerides Day</param> ''' <returns>Equatorial coordinar</returns> Public Function NutationCorrection(ByVal EC As Equatorial, ByVal JDE As Decimal) As Equatorial Dim corr As New Equatorial 'Get Nutation parameters Dim nut As Nutation = GetNutation(JDE) Dim Eps As Decimal = nut.TrueObliquity Dim dPsi As Decimal = nut.dPsi Dim dEps As Decimal = nut.dEps Dim RA As Decimal = EC.RightAscension Dim Dec As Decimal = EC.Declination 'Correction in Right Ascension and declination (In radians) corr.RightAscension = (Math.Cos(Eps) + Math.Sin(Eps) * Math.Sin(RA) * Math.Tan(Dec)) * dPsi - (Math.Cos(RA) * Math.Tan(Dec)) * dEps corr.Declination = (Math.Sin(Eps) * Math.Cos(RA)) * dPsi + Math.Sin(RA) * dEps Return corr End Function End Class
3. Aberration of Light
The third correction we need to apply id the correction for Aberration of Light. Aberration is caused because Earth itself is not a stationary object. It revolves around the Sun with a spped of app. 30km/s. Solar system itself rotates around the center of our galaxy with a speed of about 220km/s (that is 792.000 kilometers per hour, or some 430.000 knots). our galaxy is not stationary either.
The best way to explain aberration of ligh is using the analogy of rain. Imagine we are standing outside on a windless albeit rainy day. If we are standing still and not moving forward. the rain drops appear to fall down vertically. However, if we move at a speed say 10 km/h (say we're running for shelter), the rain would APPEAR to fall at an angle. The faster we move, greater the angle. Driving really fast with a car, like Formula 1, and it appears that the rain is falling horizontally towards the car.
Same goes with light. Light has a velocity of about 300.000 km/s. constructing velocity triangles just for a moving Earth, and its 30 km/s speed of rotation around the sun, and try to construct a velocity triangle, we could calculate an angle as 89°59' 39,4", which is not 90°. the difference, 20,62648", is very close to accepted constnant for annual aberration κ which equals 20,49552" (arc-seconds).
To calculate corrections in RA and Declination for A given JDE, we proceed as follows:
Calculate Julian centuries from J2000:
T = (JDE - 2451545) / 36525
Calculate Eccentricity of Earth Orbit e:
e = 0.016708617 - 0.000042037 * T - 0.0000001236 * T ^ 2
Calculate Longitude of Perichelion of this orbit π:
π = 102.93735 + 1.71953 * T + 0.00046 * T ^ 2
Mean Longitude of the Sun is L0:
L0 = 280.46645 + 36000.76983 * T + 0.0003032 * T ^ 2
Mean Anomaly of the Sun M:
M = 357.5291 + 35999.0503 * T - 0.0001559 * T ^ 2 - 0.00000048 * T ^ 3
Sun's Equation of Center C:
C = (1.9146 - 0.004817 * T) * Sin(M) + (0.019993 - 0.000101 * T) * Sin(2 * M) + 0.00029 * Sin(3 * M)
True Longitude of the Sun θ:
θ = L0 + C
Corrections for RA and Declination are then calculated as Follows:
ΔRA = -k * ((Cos(RA) * Cos(θ) * Cos(ε) + Sin(RA) * Sin(θ)) / Cos(Dec)) + + e * k * ((Cos(RA) * Cos(π) * Cos(ε) + Sin(RA) * Sin(π)) / Cos(Dec))
ΔDec = -k * (Cos(θ) * Cos(ε) * (Tan(ε) * Cos(Dec) - Sin(RA) * Sin(Dec)) + Cos(RA) * Sin(Dec) * Sin(θ)) + + e * k * (Cos(π) * Cos(ε) * (Tan(ε) * Cos(Dec) - Sin(RA) * Sin(Dec) + Cos(RA) * Sin(Dec) * Sin(π)))
Where ε is true Obliquity of the ecliptic.
The code to calculate the correction for Aberration is inside a Aberration Class:
Public Class Aberration Public Function FindCorrection(ByVal JDE As Decimal, ByVal EC As Equatorial) As Equatorial Dim corr As New Equatorial Dim nut As New IAU2000BNutationCalculation Dim T As Decimal = (JDE - 2451545) / 36525 'Mean Longitude of the sun Dim L0 As Decimal = (280.46645 + 36000.76983 * T + 0.0003032 * T ^ 2) * Constants.dtr 'Mean Anomaly of the Sun Dim M As Decimal = (357.5291 + 35999.0503 * T - 0.0001559 * T ^ 2 - 0.00000048 * T ^ 3) * Constants.dtr 'Sun Equation of Center Dim C As Decimal = ((1.9146 - 0.004817 * T) * Math.Sin(M) + (0.019993 - 0.000101 * T) * Math.Sin(2 * M) + 0.00029 * Math.Sin(3 * M)) * Constants.dtr 'True Longitude of the sun Dim Theta As Decimal = L0 + C 'kappa Dim k As Decimal = 20.49552 * Constants.dtr / 3600 'Eccentricity of Earth Orbit Dim e As Decimal = 0.016708617 - 0.000042037 * T - 0.0000001236 * T ^ 2 'True Obliquity of the ecliptic Dim Eps As Decimal = nut.GetNutation(JDE).TrueObliquity 'Longitude of Perichelion of Earth's orbit Dim PI As Decimal = (102.93735 + 1.71953 * T + 0.00046 * T ^ 2) * Constants.dtr Dim RA As Decimal = EC.RightAscension Dim Dec As Decimal = EC.Declination corr.RightAscension = -k * ((Math.Cos(RA) * Math.Cos(Theta) * Math.Cos(Eps) + Math.Sin(RA) * Math.Sin(Theta)) / Math.Cos(Dec)) + e * k * ((Math.Cos(RA) * Math.Cos(PI) * Math.Cos(Eps) + Math.Sin(RA) * Math.Sin(PI)) / Math.Cos(Dec)) * Constants.dtr corr.Declination = -k * (Math.Cos(Theta) * Math.Cos(Eps) * (Math.Tan(Eps) * Math.Cos(Dec) - Math.Sin(RA) * Math.Sin(Dec)) + Math.Cos(RA) * Math.Sin(Dec) * Math.Sin(Theta)) + e * k * (Math.Cos(Theta) * Math.Cos(Eps) * (Math.Tan(Eps) * Math.Cos(Dec) - Math.Sin(RA) * Math.Sin(Dec) + Math.Cos(RA) * Math.Sin(Dec) * Math.Sin(Theta))) * Constants.dtr Return corr End Function End Class
In the next post, we will recap all 5 posts so far on some working examples. In order to get a better picture of where on our journey towards Nautical Almanac we got so far. Stay tuned!