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
Correcting coordinates for Light time
Sun emits light by itself. All other bodies in our solar system don't. What we see is Sun's light reflection from them. The speed of light is roughly 300.000 km/s (299.792,458 km/s to be exact). Above sample can quickly show us, that the distance between Jupiter and Earth equals 848.056.265,11 kilometers. That is some 21200 times Earth circumference, which roughly equals 40.000 km. Light needs to travel from the planet to reach Earth. in the example above, light travel from Jupiter to Earth would last some 2828,81 seconds, or roughly 47 minutes. Even if Earth would be stationary, the calculation above defines coordinates as per the exact instance given. In other words. The position od Jupiter, calculated above, is not the position, which we, on Earth, see at January 17, 2006 at 12:00:00 TBD. It is the position, which we WILL see some 47 minutes later. What we need is XYZ coordinates for the jupiter at some earlier instance. The problem we face is that we do not know in advance, what is the correct Light travel time between the Celestial body and the Earth. We solve this by iteration
- Calculate XYZ coordinates for Body and the Earth, as shown above.
- Calculate the distance between Body and the Earth by using formula: Dist = SQRT(X*X + Y*Y + Z*Z)
- Get Light travel time (lt) for the distance calculated above.
- Subtract light from Julian day.
- Recalculate Cordinates XYZ for the Body. (but not Earth)
- Repeat steps 2 - 5 until a convergence in XYZ positions is achieved
I've noticed that, in practice, only a single pass is usually required, as more than sufficient accuracy is achieved for our purposes. Let us repeat the calculation and include the light time correction:
To calculate the distance between Jupiter and Earth, we will use Geocentric coordinates XB, YB, and ZB, as they were calculated above. The distance between the two bodies is then:
Dist = SQRT(XB * XB + YB * YB + ZB * ZB) = 848.056.265,12 km.
Light travel time for this distance (based on the light speed of 299.792,458 km/s equals 2828,81 seconds. There are 86400 seconds in each day, this would amount to 0,03274087045739 day. Subtracting this from Julian day in question ( 2453753,0) we get the instance at which we need to calculate SSB rectangular coordinates for Jupiter: 2453752,96725913. It is this Julian day we use to calculate new SSB coordinates for Jupiter:
- Xb = - 659.580.860,36
- Yb = - 442.214.938,44
- Zb = - 173.492.165,44
Implementing these values instead of the ones obtained before, our geocentric equatorial coordinates become:
- RA = 3,90077512348832rad or 223,4979514 deg or in proper notation 14h 53m 59,5s
- Declinatination = -0,269387177 rad or -15,43474831 deg, and in the proper notation S 15° 26' 5,09"
These values represent Astrometric Equatorial Coordinates for Jupiter on January 17,2006 at 12:00:00 TDB FOR THE REFERENCE FRAME IERS or J2000.
The complete code, which 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(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