# Slovenian Nautical Almanac - Part 5 - From Astrometric to Apparent coordinates

11.03.2016
Aljoša Kocen

Wellcome to the 6th Port of Call on our journey to Slovenian Nautical Almanac. Things we've covered so far:

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 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. 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
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):

1. 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`
2. 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`
3. 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`
4. 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`
5. 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:

1. 5 Delaunay argument multipliers (M1 to M5)
2. 3 Longitude coefficients (Lon1 to Lon3), measured in 10000000 of Arc seconds
3. 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```
In this class, we define Mean obliquity of the Ecliptic, True Obliquity of the Ecliptic, and Δε and Δψ parameters.

Next, we define `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```
To create an array of `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```

In order to get proper corrections, we define `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!