In Part 2 We calculated SSB cetred Rectangular coordinates for planets. We are now able to calculate X, Y, and Z coordinates for any instance in TDB for any planet and Sun. In this post, we will convert these coordinates to Geocentric Equatorial Coordinates, namely:

• Right Ascension, which is an angle, measured along celestial equator, between celestial body and First point of Aries (i.e. point, where plane of Ecliptic crosses the plane of Equator). it is measured in Hours, minutes, and seconds in eastward direction
• Declination, an angle between Body and celestial equator. It is measured in angular degrees, minutes, and seconds, either Noth or South of Celestial Equator

Given that we already know rectangular coordinates on the same plane, the calculation is rather straightforward:

• RA = Atan2(Y, X)
• Dec = Atan(Z / Sqrt(X * X + Y * Y)

Note that we are using Atan2 function in the calculation of Right Ascension, instead of a simple Atan function. The rational behind this is that Atan2 function automatically registers the correct quadrant, in which the angle is.  The formula for computing tangens of an angle in rectangular coordinate system is tan(φ) = Y / X.  Now let's say that we measure angle in clock-wise direction, with point, say (0,1) being Zero. If we have coordinates (X,Y) = (1,1), formula above would yield 1, and corresponding angle would be 45°. Bit if the coordinates would be (-1,-1), equation would also return 1, and Angle measured would also be 45°. In reality, the angle would be 225°. 

Net Framework Math.Atan2 function returns the angle in radians between -π/2 and +π/2. However, because RA is measured in HMS notation, going around full circle, from 0h to 24h, we need to translate this into a positive number. Putting it simply, if RA is negative by using equation above, add 2π to it.

Reducing SSB centered coordinates to Earth centered coordinates.

The coordinates, given by JPL DE Ephemerides are based on Solar System Barycenter (Solar system Center of Mass). To get RA and Declination, we need to reduce these coordinates to coordinate system, based on center of the Earth.  The problem we are facing now is that JPL DE Ephemerides enable us getting XYZ values for Earth - Moon Barycenter.  The moon itself does not revolve strictly around the center of the Earth. Instead, it revolves around the center of mass of the  Earth - Moon Couple. To get true Earth XYZ Coordinates from tabulated XYZ coordinates of EMB, we need to perform what is called Barycentring.

• From JPL DE Ephemeris, obtain Celestial Body coordinates Xb, Yb, and Zb
• From JPL DE Ephemerides, obtain Earth-Moon-Barycenter coordinates Xemb, Yemb, and Zemb
• From JPL DE Ephemerides, obtain Geocentric rectangular coordinates of the Moon, namely Xm, Ym, and Zm.

All coordinates will be returned in kilometers. We also need to know the Earth-Moon Mass ratio (EMR), which is 81,3005800000000044

To obtain Rectangular coordinates for the Earth (Xe, Ye, and Ze), we look at the solution here: 

• Xe = Xemb - (Xm / (1 + EMR))
• Ye = Yemb - (Ym / (1 + EMR))
• Ze = Zemb - (Zm / (1 + EMR))

Now we have rectangular coordinates, based on Solar system BaryCenter for Earth.  Let us define geocentric coordinates for the Body as XB, YB, and ZB. We proceed by calculating them as follows:

• XB = Xb - Xe
• YB = Yb - Ye
• ZB = Zb - Ze

Now to calculate Right ascension and Declination for the body, we use the two formulas above. 

• RA = Atan2(YB / XB)
• Dec = Atan(ZB / Sqrt(YB * YB + XB *XB)) 

An Example:

Let us calculate Right Ascension and Declination for Jupiter on January 17, 2006 at 12:00:00 TDB:

Julian Day for the set date and time of day equals 2453753,0. Calculated coordinates for EMB, Jupiter, and Moon for this instance are as follows:

First we will calculate Earth coordinates:

• Xe = Xemb - (Xm / (1 + EMR)) = -66.596.740.39 - (-354.436,33 / (1 + 81,30058)) = -66.592.433,78
• Ye = Yemb - (Ym / (1 + EMR)) = -120.471.863,66 - (172.075.13 / (1 + 81,30058)) = 120.469.772,85
• Ze = Zemb - (Zm / (1 + EMR))=  52.20.9.922.09 - (97.343.28 / (1 + 81,30058)) = 52.208.739,31

To obtain Gocentric coordinates for Jupiter:

• XB = Xb - Xe = -659.559.726,77 - (-66.592.433,78) = -   592.967.292,99
• YB = Yb - Ye = -442.240.720.21 - 120.469.772,85 =  -   562.710.493,06
• ZB = Zb - Ze = -173.503.731,00 - 52.208.739,78 = -   225.712.470,31

Calculated Geocentric coordinates for Jupiter would therefore yield:

• RA: -2,38236951 Radians, truncated to positive values equal 3,90081580 or in Hour-Minute-Second (HMS) Notation:  14h 54m 0,07s
• Declination: -0,26939956 Radians, which is  -15,43545806 degrees, or in proper notation:  S 15° 26' 07,65"

Our VB.Net code to calculate this is shown below:

    Private Sub CalculatePosition(ByVal TBD As DateTime, ByVal CelBody As CelestialBody)
        'PositionCalculation Object
        Dim PosCalc As New PositionCalculation

        'CartesianCoordinates Objects for Body, EMB and Moon
        Dim Body As New CartesianCoordinates
        Dim EMBPos As New CartesianCoordinates
        Dim Moon As New CartesianCoordinates

        'Julian Day Number
        Dim JD As Decimal = New DateTime(dtpdateTime.Value.Year,
                                         dtpdateTime.Value.Month,
                                         dtpdateTime.Value.Day,
                                         12, 0, 0).ToJulianDay

        'SSB based Body Coordinates
        Body = PosCalc.GetBodyPosition(CelBody, JD)

        'SSB EMB Coordinates
        EMBPos = PosCalc.GetBodyPosition(bodies(2), JD)

        'Geocentric position of the moon
        Moon = PosCalc.GetBodyPosition(bodies(9), JD)

        'SSB based Earth coordinates
        Dim Earth As New CartesianCoordinates


        Earth.X = EMBPos.X - Moon.X / (1 + Constants.EarthMoonMassRatio)
        Earth.Y = EMBPos.Y - Moon.Y / (1 + Constants.EarthMoonMassRatio)
        Earth.Z = EMBPos.Z - Moon.Z / (1 + Constants.EarthMoonMassRatio)


        'Earth based rectangular coordinates of the body.
        Dim BodyEarth As CartesianCoordinates = Body - Earth

        Dim X As Decimal = BodyEarth.X
        Dim Y As Decimal = BodyEarth.Y
        Dim Z As Decimal = BodyEarth.Z

        'Astrometric Coordinates
        Dim Dec As Decimal = Math.Atan(Z / Math.Sqrt(X ^ 2 + Y ^ 2))
        Dim RA As Decimal = Math.Atan2(Y, X)

        'Truncate Right Ascension if negative, or if larger than 2*PI. 
        RA = MathFunctions.TruncateRad(RA)

    End Sub

	Private Sub CalculatePosition(ByVal TBD As DateTime, ByVal CelBody As CelestialBody)
		'PositionCalculation Object
		Dim PosCalc As New PositionCalculation

		'CartesianCoordinates Objects for Body, EMB and Moon
		Dim Body As New CartesianCoordinates
		Dim EMBPos As New CartesianCoordinates
		Dim Moon As New CartesianCoordinates

		'Julian Day Number
		Dim JD As Decimal = New DateTime(TBD.Year,
										 TBD.Month,
										 TBD.Day,
										 12, 0, 0).ToJulianDay

		'SSB based Body Coordinates
		Body = PosCalc.GetBodyPosition(CelBody, JD)

		'SSB EMB Coordinates
		EMBPos = PosCalc.GetBodyPosition(bodies(2), JD)

		'Geocentric position of the moon
		Moon = PosCalc.GetBodyPosition(bodies(9), JD)

		'SSB based Earth coordinates
		Dim Earth As New CartesianCoordinates


		Earth.X = EMBPos.X - Moon.X / (1 + Constants.EarthMoonMassRatio)
		Earth.Y = EMBPos.Y - Moon.Y / (1 + Constants.EarthMoonMassRatio)
		Earth.Z = EMBPos.Z - Moon.Z / (1 + Constants.EarthMoonMassRatio)


		'Earth based rectangular coordinates of the body.
		Dim BodyEarth As CartesianCoordinates = Body - Earth

		Dim X As Decimal = BodyEarth.X
		Dim Y As Decimal = BodyEarth.Y
		Dim Z As Decimal = BodyEarth.Z

		'Astrometric Coordinates
		Dim Dec As Decimal = Math.Atan(Z / Math.Sqrt(X ^ 2 + Y ^ 2))
		Dim RA As Decimal = Math.Atan2(Y, X)

		'Truncate Right Ascension if negative, or if larger than 2*PI. 
		RA = MathFunctions.TruncateRad(RA)

		'In order to correct for light travel time, we add the following:

		'Distance between earth and Jupiter in kilometers
		Dim Distance As Decimal = BodyEarth.Distance

		'Light travel time:
		Dim lt As Decimal = Distance / (299792.458 * 86400)

		'Get a modified Julian Day (MJD)

		Dim MJD As Decimal= JD - lt

		'Get New Body position for Modified Julian Day
		Body=PosCalc.GetBodyPosition(CelBody,MJD)
		
		'Geocentric XYZ coordinates:
		BodyEarth= Body - Earth

		X=BodyEarth.X
		Y=BodyEarth.Y
		Z=BodyEarth.Z

		'Astrometric geocentric coordinates (in radians)
		Dec = Math.Atan(Z / Math.Sqrt(X ^ 2 + Y ^ 2))
		RA = MathFunctions.TruncateRad(Math.Atan2(Y, X))

	End Sub

WHat the hell is TDB, and on which clock is measured.