In the process of obtaining the data for the Ephemerides calculation, I really only had two options available:
- VSOP87 data, constructed by Pierre Bretagnon et. al. This data is based on an elaborate curve fitting to results, based on calculations, provided by NASA JPL DE200 Ephemerides data.
- NASA JPL Developement Ephemerides. These are based on laser, radio and telescopic measurements of Sun, Moon, planets, Asteroids and other celestial bodies.
Knowing that the first option is really based on the second one, albeit with a different edition of it (VSOP87 uses numerical Integration of JPL DE200 Ephemerides), and it does not incorporate Moon data, I decided to go forward with NASA JPL DE405 Ephemerides data.
Description of JPL DE 405 Ephemerides data
NASA JPL Ephemerides can be found at this publicly available FTP address. The format I chose, is ASCII, as it is universally readable, regardless of the hardware one uses. They all list Coefficients for computing Solar System Barycentric (SSB) rectangular planetary coordinates for the reference frame J2000:
- Mercury
- Venus
- Earth-Moon Barycenter (EMB)
- Mars
- Jupiter
- Saturn
- Uranus
- Pluto
- Moon (Geocentric)
- Sun
- Nutations
- Librations
The last two (Nutations and Librations) are really not needed in Ephemerides calculations, but they are included in the datasets. As you will see, the nutation is an important part at calculating the apparent position of the celestial body. We will get to that later.
JPL DE Ephemerides provide Coefficients for calculating Rectangular coordinates based on Solar system barycenter.
- X-axis points toward Vernal Equinox or First point of Aries
- Y-axis points toward RA = 6h or 90°
- Z-axis is perpendicular to both X, and Y axis
All three coordinates are given in kilometers. The X and Y axis lie on the plane of Celestial Equator. Z axis is perpendicular to it.
N.B.: We shall leave Nutations and Librations out of our calculations for now.
There are a number of different versions of data. I went with DE405 version. It covers dates from 9.12.1599 AD to 20.2.2201 AD. Complete series is divided into a number of text files, namely "ascp1600.405" or "ascp2000.405", each covering roughly 20 years. Files do not have ".txt" extension, but no need to worry, they can be opened in any text editor. Here's a screenshot of Notepad:
Each file consists of several data blocks, lasting 32 days each. in the pic above, you can see "1 1018" in the first line. It denotes the sequence number of the Data Block inside a file. And the second number denotes the number of coefficients inside DataBlock. The above means that we are looking at first DataBlock in the file and that there are 1018 coefficients in the DataBlock.
Each DataBlock begins with two numbers, which determine the date validity of the DataBlock. The first two numbers in the second row we have: 0.247345650000000000D+07 and 0.247348850000000000D+07. These denote the first and the last JULIAN DAY instances for which coefficients should be used. The third annd all consecutive numbers denote 1018 Coefficients. These can be obtained for individual celestial bodies as follows:
To achieve desired accuracy, DataBlock is sub-divided into intervals. How many intervals there are for different celestial body, and how many coefficients there are in each sub interval is defined in "header.405" file. I'll list these numbers for you in the following table:
Coefficients per SubInterval are the same for each of the coordinates (X,Y, or Z). To represent DataBloch Schematically in a "tree", one could write:
- DataBlock (1018 Records)
- Start Julian day for the DataBlock
- End Julian day for the DataBlock
- Mercury Data
- SubInterval 1 of 4
- 14 Coefficients for X Coordinate
- 14 Coefficients for Y Coordinate
- 14 Coefficients for Z Coordinarte
- SubInterval 2 of 4
- 14 Coefficients for X Coordinate
- 14 Coefficients for Y Coordinate
- 14 Coefficients for Z Coordinate
- SubInterval 3 of 4
- 14 Coefficients for X Coordinate
- 14 Coefficients for Y Coordinate
- 14 Coefficients for Z Coordinate
- SubInterval 4 of 4
- 14 Coefficients for X Coordinate
- 14 Coefficients for Y Coordinate
- 14 Coefficients for Z Coordinate
- SubInterval 1 of 4
- Venus Data
- SubInterval 1 of 2
- 10 Coefficients for X Coordinate
- 10 Coefficients for Y Coordinate
- 10 Coefficients for Z Coordinarte
- SubInterval 2 of 2
- 10 Coefficients for X Coordinate
- 10 Coefficients for Y Coordinate
- 10 Coefficients for Z Coordinate
- SubInterval 1 of 2
- Earth - Moon BaryCenter Data
- ......
And so fort for all the other bodies, and other data. As the data is so beautifully structured, we can easily obtain Coefficients set for any Julian Day and for any body in question.
But first let's convert this somewhat awkward three column layout into a more descriptive single column layout. I will show you how to convert the data inside ".405" files in a single column Excel 2007/2010/2013/2016 files. Note however, that you don't need to do that, you can write your own logic to extract the aprropriate coefficients.
Converting ASCII files to Excel Format
Looking at the Notepad screenshot above, you see, that ASCII files are formatted in a somewhat awkward three-column layout. It' did not bother me really, as it is quite easy to traverse the columns and rows quite easily. But I must admit it frustrated me a bit, when trying to look-up a certain number in the great mess of hudredrs of thousands of records. Our calculation of SetIntervals and SetIndexes are better visualised right?
What we basically need to do is:
- Open ascpXXXX.405 file (XXXX being the "year" in the file name)
- Declare auxiliarry list, which will hold the data. I used System.Data.DataTable object.
- Loop through the ASCII file line by line and:
- Trim excessive spaces from the line
- split the line into array of data
- If array has only two records, and the second one equals "1018", skip complete line
- In any other case, we have three numbers, each of which is inserted into the DataTable
- With the aid of EPPlus Excel library, create an empty Excel Workbook
- Load the resulting DataTable into an empty Excel Worksheet
- Save the Workbook and clean up the memory allocated
The following routine does just that. It is pretty well documented in itself:
Private Sub ProcessFile(ByVal FileName As String)
'Auxilliary Results DataTable
Dim ResultsTable As New DataTable
ResultsTable.Columns.Add("Value", Type.GetType("System.String"))
Dim Line As String = ""
Dim objreader As System.IO.StreamReader
Dim arr() As String
Try
If System.IO.File.Exists(FileName) Then
Dim Label As Integer
objreader = New System.IO.StreamReader(FileName)
Do While objreader.Peek <> -1
'load Line of text from ASCII File
Line = objreader.ReadLine()
'Trim leading and trailing spaces, if any
Line = Line.Trim(" ")
'Remove multiple spaces between data
Line = Regex.Replace(Line, " {2,}", " ")
'Splot the line into array
arr = Line.Split(" ")
If arr(1) = "1018" Then
'DataBlock Header
Label = CInt(arr(0))
Else
'Actual data inside a DataBlock
ResultsTable.Rows.Add(arr(0).ToString.Replace("D", "E"))
ResultsTable.Rows.Add(arr(1).ToString.Replace("D", "E"))
ResultsTable.Rows.Add(arr(2).ToString.Replace("D", "E"))
End If
Loop
'Save to Excel
Dim ep As New ExcelPackage
Dim wb As ExcelWorkbook = ep.Workbook
Dim sheet As ExcelWorksheet = wb.Worksheets.Add("Data")
'Load the data into the Sheet without headers of the DataTable ("Value")
sheet.Cells(1, 1).LoadFromDataTable(ResultsTable, False)
Dim OutFileName As String = ""
'Save the Excel file inside the same folder as original ASCII Ephemerides file.
OutFileName = Path.GetDirectoryName(FileName) & "\"
OutFileName = OutFileName & Path.GetFileNameWithoutExtension(FileName) & ".xlsx"
Dim fi As System.IO.FileInfo = New FileInfo(OutFileName)
ep.SaveAs(fi)
'Clean up
sheet =Nothing
wb=Nothing
ep=Nothing
ResultsTable = Nothing
GC.Collect()
End If
Catch ex As Exception
Throw New NotImplementedException
End Try
End Sub
One thing to note in this code is, that we are literally writting String values to Excel Worksheet. There is a good reasoning on that one. If we tried to write actual numbers inside Excel worksheet, Excel would not register ALL decimal numbers. For example value " -0.116228947221808049E-05" would literally be truncated into: "-0.116228947221808000E-05" The last three significant digits get converted into zeros. As we are dealing with some pretty large nimbers here, I would not recomend zeroing out anything.
In case you want to convert 3-column layout in original .405 Ascii files into a single column text files, you can use similar routine, except for the fact that you are writing into the text file:
Public Sub ProcessFileToText(ByVal FileName As String)
Dim Line As String = ""
Dim objreader As System.IO.StreamReader
Dim OutFileName As String = ""
'Save in same Dir as input .405 file
OutFileName = Path.GetDirectoryName(FileName) & "\"
OutFileName = OutFileName & Path.GetFileNameWithoutExtension(FileName) & ".txt"
Dim objWriter As New System.IO.StreamWriter(OutFileName)
Dim arr() As String
Try
If System.IO.File.Exists(FileName) Then
Dim Label As Integer
objreader = New System.IO.StreamReader(FileName)
Do While objreader.Peek <> -1
'load Line of text from ASCII File
Line = objreader.ReadLine()
'Trim leading and trailing spaces, if any
Line = Line.Trim(" ")
'Remove multiple spaces between data
Line = Regex.Replace(Line, " {2,}", " ")
'Splot the line into array
arr = Line.Split(" ")
If arr(1) = "1018" Then
'DataBlock Header
Label = CInt(arr(0))
Else
'Actual data inside a DataBlock
objWriter.WriteLine(arr(0).ToString)
objWriter.WriteLine(arr(1).ToString)
objWriter.WriteLine(arr(2).ToString)
End If
Loop
objreader.Close()
objWriter.Close()
End If
Catch ex As Exception
End Try
End Sub
As you can see, the code is nearly the same. except for writing to DataTable object and then to Excel WorkSheet, we write to output stream while reading and parsing the input stream.
Calculating Solar System barycenter Rectagular coordinates of a body.
Let us Find Coefficients for calclationg position of Earth-Moon baryCenter for say 20.2.2016 at 00:00:00 TDB (More about different time scales in the next chapter):
Julian Day for 20.2.2016 00:00:00: 2457438,5
Solution:
- First let's find the appropriate Ascii file. From the filename, we can safely assume, that ascp2000.405 is the one we're looking for, as it starts at around year 2000 and covers a 20-year period. the file covers Julian days from 2451536,5 to 2458864,5.
- Inside this File we need to find the exact DataBlock (from 1 to 229), that we need to look for. We know, that each DataBlock covers 32 day period. There are exactly 5904 days. If each DataBlock covers 32, our data would be in DataBlock number:
DataBlockNumber = INT((JD-StartJD) / 32)+1 = 185
We check this, by searching for the 185th DataBlock Inside a file, and find, that DataBlock Number 185 covers Julian days from 2457438,5 to 2457456,5. Our JD falls inside this range.
- The first two numbers in the DataBlock are Julian Day numbers. So we skip them. Then we must skip all the data for Mercury, and Venus. Looking at the table above, Mercury has 3 coordinates (X,Y,Z), 4 SubIntervals, and 14 Coefficients for each coordinate per SubInterval. Therefore we need to skip
NumberToSkip = NumberOfSubIntervals * NumberOfCoordinates * NumberofCoefficients
for Mercury and Vens, to get to the start of EMB data. For Mercry we need to skip next 168 records, and for Venus Further 60, bringing total to 168+60+2 (Julian days in the begining of the DataBlock) = 230 records. - Given that Earh Moon Barycenter Data is sub-divided into 2 subintervals, we need to define into which SubInterval our Julian Day (2457438,5) fits. Three SubIntervals mean, that each covers exactly 16 days. JD-JD-Block equals to exactly 14 days. We calculate the SubIntervalIndex as:
SubIntervalIndex = FLOOR(JD-StartJD)/32
In our case using 2457438,5 and 2457424,5 for our JD and DataBlock starting JD respectivelly, this yields zero. So we grab the first 13 numbers for coordinate X, 13 numbers for coordinatesY, and the 13 numbers for the Z coordinate. As shown in the table, we got the following set of coefficients:
Now let's translate this into code. We will use Excel version of the data. First we'll declare CelestialBody object, which holds the following data:
- ID
- Name of the Body
- Number of Coefficients Sets
- Number of Coefficients per Set
- Number of Coordinates
- Starting index for each Body
Public Class CelestialBody Public Property ID As Integer Public Property Name As String Public Property NumberOfSets As Integer Public Property CoefficientsPerSet As Integer Public Property NumberOfCoordinates As Integer Public Property StartRead As Integer Public Sub New() End Sub Public Sub New(ByVal ID As Integer, ByVal Name As String, ByVal NumberOfSets As Integer, ByVal CoefficientsPerSet As Integer, ByVal NumberOfCoordinates As Integer) Me.ID = ID Me.Name = Name Me.NumberOfSets = NumberOfSets Me.CoefficientsPerSet = CoefficientsPerSet Me.NumberOfCoordinates = NumberOfCoordinates End Sub End Class
Public Class BodiesList Inherits List(Of CelestialBody) Public Sub New() 'Add the data for all Bodies, Nutations, and Librations Me.Add(New CelestialBody(1, "Mercury ", 4, 14, 3)) Me.Add(New CelestialBody(2, "Venus ", 2, 10, 3)) Me.Add(New CelestialBody(3, "Earth Moon BaryCenter", 2, 13, 3)) Me.Add(New CelestialBody(4, "Mars ", 1, 11, 3)) Me.Add(New CelestialBody(5, "Jupiter ", 1, 08, 3)) Me.Add(New CelestialBody(6, "Saturn ", 1, 07, 3)) Me.Add(New CelestialBody(7, "Uranus ", 1, 06, 3)) Me.Add(New CelestialBody(8, "Neptune ", 1, 06, 3)) Me.Add(New CelestialBody(9, "Pluto ", 1, 06, 3)) Me.Add(New CelestialBody(10, "Moon - Geocentric ", 8, 13, 3)) Me.Add(New CelestialBody(11, "Sun ", 2, 11, 3)) Me.Add(New CelestialBody(12, "Nutations ", 4, 10, 2)) Me.Add(New CelestialBody(13, "Librations ", 4, 10, 3)) 'Mercury Data Starts at index 2 of Each DataBlock Me(0).StartRead = 2 'All subsequent StartReadValues are calculated: For i As Integer = 1 To Me.Count - 1 Me(i).StartRead = Me(i - 1).StartRead + (Me(i - 1).NumberOfCoordinates * Me(i - 1).NumberOfSets * Me(i - 1).CoefficientsPerSet) Next End Sub End Class
EphemDataBlock:''' <summary>
''' Represents JPL DE-XXX DataBlock
''' </summary>
Public Class EphemDataBlock
Inherits List(Of Decimal)
''' <summary>
''' StartingJulian Day for the dataBlock
''' </summary>
''' <returns>Decimal value</returns>
Public Property StartBlockJD As Decimal
''' <summary>
''' Ending day day for the dataBlock
''' </summary>
''' <returns></returns>
Public Property EndBlockJD As Decimal
Public Property NumberOfBlocks As Integer
Public Function GetSubBlock(ByVal JDE As Decimal) As EphemDataBlock
Dim SubBlock As New EphemDataBlock
If Me.Count = 0 Then
Throw New Exception("Nothing in the dataBlock, cannot extract")
ElseIf (JDE > EndBlockJD Or JDE < StartBlockJD)
Throw New Exception("Julian day out of bounds for the current block")
Else
'Get the block number for the JDE in question
'Each SubBlock has 1020 records (First 2 records for JulianDays and 1018 coefficients)
'Get the DataBlock Index
Dim DataBlockIndex As Long = Math.Floor((JDE - Me.StartBlockJD) / 32)
'StartIndex
Dim StartIndex As Long = DataBlockIndex * 1020
For i As Long = StartIndex To StartIndex + 1019
SubBlock.Add(Me(i))
Next
SubBlock.StartBlockJD = SubBlock(0)
SubBlock.EndBlockJD = SubBlock(1)
Return SubBlock
End If
End Function
End Class
List(Of Decimal) object, and has one Function, GetSubBlock, which accepts JulianDay as parameter and returns a single 1020 record long basic DataBlock, which covers basic 32 day period.EphemDataBlock object:Imports System.Globalization
Imports System.IO
Imports OfficeOpenXml
Public Class DataFileRead
Public Shared Function ReadDataFile(ByVal JulTime As Decimal) As EphemDataBlock
Dim Block As New EphemDataBlock
Dim records As Integer
Dim StartingJD As Decimal
Dim EndingJD As Decimal
Dim FileName As String = ""
Dim FolderName As String="PathToFolderWithExcelData"
'Find the appropriate Excel File, based on Input Julian Day number - JulTime
If ((JulTime >= 2414992.5) _
AndAlso (JulTime < 2422320.5)) Then
StartingJD = 2414992.5
EndingJD = 2422320.5
FileName = "ASCP1900.xlsx"
records = 230
ElseIf ((JulTime >= 2422320.5) _
AndAlso (JulTime < 2429616.5)) Then
StartingJD = 2422320.5
EndingJD = 2429616.5
FileName = "ASCP1920.xlsx"
records = 229
ElseIf ((JulTime >= 2429616.5) _
AndAlso (JulTime < 2436912.5)) Then
StartingJD = 2429616.5
EndingJD = 2436912.5
FileName = "ASCP1940.xlsx"
records = 229
ElseIf ((JulTime >= 2436912.5) _
AndAlso (JulTime < 2444208.5)) Then
StartingJD = 2436912.5
EndingJD = 2444208.5
FileName = "ASCP1960.xlsx"
records = 229
ElseIf ((JulTime >= 2444208.5) _
AndAlso (JulTime < 2451536.5)) Then
StartingJD = 2444208.5
EndingJD = 2451536.5
FileName = "ASCP1980.xlsx"
records = 230
ElseIf ((JulTime >= 2451536.5) _
AndAlso (JulTime < 2458832.5)) Then
StartingJD = 2451536.5
EndingJD = 2458832.5
FileName = "ASCP2000.xlsx"
records = 229
ElseIf ((JulTime >= 2458832.5) _
AndAlso (JulTime < 2466128.5)) Then
StartingJD = 2458832.5
EndingJD = 2466128.5
FileName = "ASCP2020.xlsx"
records = 229
ElseIf ((JulTime >= 2466128.5) _
AndAlso (JulTime < 2473456.5)) Then
StartingJD = 2466128.5
EndingJD = 2473456.5
FileName = "ASCP2040.xlsx"
records = 230
ElseIf ((JulTime >= 2473456.5) _
AndAlso (JulTime < 2480752.5)) Then
StartingJD = 2473456.5
EndingJD = 2480752.5
FileName = "ASCP2060.xlsx"
records = 229
ElseIf ((JulTime >= 2480752.5) _
AndAlso (JulTime < 2488048.5)) Then
StartingJD = 2480752.5
EndingJD = 2488048.5
FileName = "ASCP2080.xlsx"
records = 229
ElseIf ((JulTime >= 2488048.5) _
AndAlso (JulTime < 2495344.5)) Then
StartingJD = 2488048.5
EndingJD = 2495344.5
FileName = "ASCP2100.xlsx"
records = 229
ElseIf ((JulTime >= 2495344.5) _
AndAlso (JulTime < 2502672.5)) Then
StartingJD = 2495344.5
EndingJD = 2502672.5
FileName = "ASCP2120.xlsx"
records = 230
ElseIf ((JulTime >= 2502672.5) _
AndAlso (JulTime < 2509968.5)) Then
StartingJD = 2502672.5
EndingJD = 2509968.5
FileName = "ASCP2140.xlsx"
records = 229
ElseIf ((JulTime >= 2509968.5) _
AndAlso (JulTime < 2517264.5)) Then
StartingJD = 2509968.5
EndingJD = 2517264.5
FileName = "ASCP2160.xlsx"
records = 229
ElseIf ((JulTime >= 2517264.5) _
AndAlso (JulTime < 2524624.5)) Then
StartingJD = 2517264.5
EndingJD = 2524624.5
FileName = "ASCP2180.xlsx"
records = 230
End If
'Load Appropriate File in ExcelPackage
FileName = Path.Combine(FolderName, FileName)
Dim ep As New ExcelPackage(New FileInfo(FileName))
Dim wb As ExcelWorkbook = ep.Workbook
'Load The Worksheet
Dim Sheet As ExcelWorksheet = wb.Worksheets("Data")
'Aux variable to be parsed from Sheet String values
Dim Number As Decimal
'Declare a Worksheet row pointer
Dim iPointer As Long=1
'Loop through the worksheet until all the data is collected into the EphemDataBlock
Do
'Parse the string into the Decimal Value
Decimal.TryParse(Sheet.Cells(iPOinter, 1).Text.Replace("D","E"), NumberStyles.Float, CultureInfo.InvariantCulture, Number)
'Increase the pointer
iPointer=iPointer+1
'Add value to EphemDataBlock
Block.Add(Number)
'Loop while we still have any values inside an Excel Worksheet
Loop While Sheet.Cells(iPointer,1).Text<>""
'Julian days border values for the data block
Block.StartBlockJD=StartingJD
Block.EndBlockJD=EndingJD
'Clean up
Sheet.Dispose()
wb.Dispose()
ep.Dispose()
GC.Collect()
'Return Result
Return Block
End Function
End Class
CartesianCoordinates:''' <summary> ''' Represents a point in Cartesian Coordinate System. ''' </summary> Public Class CartesianCoordinates Public Property X As Decimal Public Property Y As Decimal Public Property Z As Decimal ''' <summary> ''' Gets the distance of CartesianCoordinates Object from Coordinate System Origin ''' </summary> ''' <returns></returns> Public Readonly property Distance Get Return Math.Sqrt(X ^ 2 + Y ^ 2 + Z ^ 2) End Get End Property ''' <summary> ''' Default, parameter-less constructor ''' </summary> Public Sub New() 'X,Y,Z default to 0 Me.X = 0 Me.Y = 0 Me.Z = 0 End Sub ''' <summary> ''' Constructor with parameters ''' </summary> ''' <param name="X">X value</param> ''' <param name="Y">Y value</param> ''' <param name="Z">Z value</param> Public Sub New(ByVal X As Decimal, ByVal Y As Decimal, ByVal Z As Decimal) Me.X = X Me.Y = Y Me.Z = Z End Sub ''' <summary> ''' Performs Vector summation of CartesianCoordinates objects ''' </summary> ''' <param name="First">First CartesianCoordinates object</param> ''' <param name="Second">Second CartesianCoordinates object</param> ''' <returns>Vector sum of both input CartesianCoordinates </returns> Public Shared Operator +(ByVal First As CartesianCoordinates, Second As CartesianCoordinates) Dim res As New CartesianCoordinates res.X = First.X + Second.X res.Y = First.Y + Second.Y res.Z = First.Z + Second.Z Return res End Operator ''' <summary> ''' Performs Vector subtraction of CartesianCoordinates objects ''' </summary> ''' <param name="First">First CartesianCoordinates object</param> ''' <param name="Second">Second CartesianCoordinates object</param> ''' <returns>Difference between First and Second CartesianCoordinates objects</returns> Public Shared Operator -(ByVal First As CartesianCoordinates, Second As CartesianCoordinates) Dim res As New CartesianCoordinates res.X = First.X - Second.X res.Y = First.Y - Second.Y res.Z = First.Z - Second.Z Return res End Operator ''' <summary> ''' Performs Scalar division of CartesianCoordinates Object by a Scalar. It's an inverse of scalar multiplication ''' </summary> ''' <param name="Cartesian">CartesianCoordinates object</param> ''' <param name="Scalar">A number, used to divide a CartesianCoordinates object with</param> ''' <returns></returns> Public Shared Operator /(ByVal Cartesian As CartesianCoordinates, ByVal Scalar As Decimal) Dim res As New CartesianCoordinates res.X = Cartesian.X / Scalar res.Y = Cartesian.Y / Scalar res.Z = Cartesian.Z / Scalar Return res End Operator ''' <summary> ''' Performs Scalar multiplication of CartesianCoordinates Object by a Scalar ''' </summary> ''' <param name="Cartesian">CartesianCoordinates object</param> ''' <param name="Scalar">A number, used to multiply a CartesianCoordinates object</param> ''' <returns></returns> Public Shared Operator *(ByVal Cartesian As CartesianCoordinates, ByVal Scalar As Decimal) Dim res As New CartesianCoordinates res.X = Cartesian.X * Scalar res.Y = Cartesian.Y * Scalar res.Z = Cartesian.Z * Scalar Return res End Operator End ClassNote from the code, that I used operator overloading for basic manipulation with the CartesianCoordinate object. That is not necessary, but look at the following code snippet:
Dim A As CartesianCoordinates, B As CartesianCoordinates, Sum As CartesianCoordinates '1. Without operator overloading Sum.X=A.X+B.X Sum.Y=A.Y+B.Y Sum.Z=A.Z+B.Z '2. With operator overloading Sum=A+B
PositionCalculation:Public Class PositionCalculation ''' <summary> ''' Full DataBlock representing Full 20-Year data, read from a file ''' </summary> ''' <returns></returns> Public Property DataBlock As New EphemDataBlock ''' <summary> ''' 1020 records for a given Julian Day in question ''' </summary> ''' <returns></returns> Public Property SubBlock As New EphemDataBlock ''' <summary> ''' Calculates SSB Rectangular coordinates for Celestial Body. ''' </summary> ''' <param name="SelectedBody">CelestialBody</param> ''' <param name="JulTime">Julian day (in Barycentric Dynamic Time)</param> ''' <returns>CartesianCoordinates object, with coordinate values in kilometers</returns> Public Function GetBodyPosition(ByVal SelectedBody As CelestialBody, ByVal JulTime As Decimal) As CartesianCoordinates Dim Data As New CartesianCoordinates 'Check For Data Call CheckDataSet(JulTime) 'Get SubBlock, based on a required Julian day, where we only need 1020 records SubBlock = DataBlock.GetSubBlock(JulTime) 'Declare auxilliary local variables Dim NumberOfSets As Integer = SelectedBody.NumberOfSets Dim CoefficientsPerSet As Integer = SelectedBody.CoefficientsPerSet Dim NumberOfCoordinates As Integer = SelectedBody.NumberOfCoordinates 'Get pointer Dim iPointer As Integer = SelectedBody.StartRead 'Duration of each SubInterval is based on the number of SubIntervals inside a DataBlock for each Body Dim SubIntervalDuration As Integer = 32 / NumberOfSets Dim SubIntervalIndex As Integer = Math.Floor((JulTime - SubBlock.StartBlockJD) / SubIntervalDuration) 'Move the pointer to the iPointer = iPointer + (SubIntervalIndex * NumberOfCoordinates * CoefficientsPerSet) 'Declare Coefficient arrays for X,Y,Z Dim CoeffsX(CoefficientsPerSet - 1) As Decimal Dim CoeffsY(CoefficientsPerSet - 1) As Decimal Dim CoeffsZ(CoefficientsPerSet - 1) As Decimal 'Load X-axis Coefficients For j As Integer = 0 To CoefficientsPerSet - 1 CoeffsX(j) = SubBlock(iPointer) iPointer += 1 Next 'Load Y-axis Coefficients For j As Integer = 0 To CoefficientsPerSet - 1 CoeffsY(j) = SubBlock(iPointer) iPointer += 1 Next 'Load Z-axis Coefficients For j As Integer = 0 To CoefficientsPerSet - 1 CoeffsZ(j) = SubBlock(iPointer) iPointer += 1 Next 'Calculate Chebyshev Time in the interval [-1...0...1] Dim ChebTime As Decimal = 2 * (JulTime - (SubIntervalIndex * SubIntervalDuration + SubBlock.StartBlockJD)) / SubIntervalDuration - 1 'Calculate Position Chebyshev Polynomial Dim PositionPoly(CoefficientsPerSet - 1) As Decimal PositionPoly(0) = 1 PositionPoly(1) = ChebTime For j As Integer = 2 To CoefficientsPerSet - 1 PositionPoly(j) = 2 * ChebTime * PositionPoly(j - 1) - PositionPoly(j - 2) Next 'Calculate BodyPosition in CartesianCoordinates Data.X = 0 For j As Integer = 0 To CoefficientsPerSet - 1 Data.X = Data.X + CoeffsX(j) * PositionPoly(j) Next Data.Y = 0 For j As Integer = 0 To CoefficientsPerSet - 1 Data.Y = Data.Y + CoeffsY(j) * PositionPoly(j) Next Data.Z = 0 For j As Integer = 0 To CoefficientsPerSet - 1 Data.Z = Data.Z + CoeffsZ(j) * PositionPoly(j) Next 'Return the result Return Data End Function ''' <summary> ''' Loads coefficients from dta file into DataBlock ''' </summary> ''' <param name="JulTime">Julian day (TDB)</param> Private Sub CheckDataSet(ByVal JulTime As Decimal) 'If there is no data in DataBlock, load appropriate data from the file If DataBlock.Count = 0 Then DataBlock = DataFileRead.ReadDataFile(JulTime) ElseIf JulTime < DataBlock.StartBlockJD Or JulTime > DataBlock.EndBlockJD Then 'Requested JulianDay is outside of the current DataBlock, so reload new set of data from the file DataBlock = DataFileRead.ReadDataFile(JulTime) End If End Sub End Class
CelestialBody object (described above) Julian day as an input and then:- It loads the
DataBlockwith full set from the Excel file, if there are no data to work with, or reloads the data, if JulianDay input falls outside the bounds of previously loaded data. - Depending on the JulianDay input, it extracts the 1020 record
SubBlock - CelestialBody object already contains the parameters needed to move the read pointer to the required position inside the DataBlock, depending on the CelestialBody, and SubInterval as described above. See lines 34 to 48
- loads three Coefficient arrays with the data from the SubBlock (lines 50 to 71)
- Calculates the Chebyshev time, which is always in the interval from -1 to 1, and uses it to get the Chebyshev polynomial terms
- Combining Chebyshev polynomial temrs with each array of Coefficients it calculates all three rectangular coordinates, and returns the result.
A sample, calculating Earth-Moon Barycenter rectangular coordinates for complete Year 2016
I've setup a small example, on how to calculate SSB Rectangular coordinates for Earth Moon Barycenter. It will return the X, Y and Z data for each day from 1.1.2016 to 1.1.2017 at 12:00 TBD and save the result in excel File
'Declare a new PositionCalculation Instance
Dim PosCalc As New PositionCalculation
'Declare Bodies List
Dim Bodies As New BodiesList
'DateTimeBrackets
Dim FirstDate As DateTime = New DateTime(2016, 1, 1, 12, 0, 0)
Dim LastDate As DateTime = New DateTime(2017, 1, 1, 12, 0, 0)
Dim CurrentDate As DateTime=FirstDate
'Resulting dataTable
Dim ResultTable As New DataTable
ResultsTable.Columns.Add("JulianDay", Type.GetType("System.DateTime"))
ResultsTable.Columns.Add("X", Type.GetType("System.Decimal"))
ResultsTable.Columns.Add("Y", Type.GetType("System.Decimal"))
ResultsTable.Columns.Add("Z", Type.GetType("System.Decimal"))
'Create new CartesianCoordinates object to hold the values
Dim Coords As New CartesianCoordinates
'Loop through the calculation
Do
'Get Coordinates
Coords = PosCalc.GetBodyPosition(Bodies(2), CurrentDate.ToJulianDay)
ResultsTable.Rows.Add(CurrentDate, Coords.X, Coords.Y, Coords.Z)
CurrentDate = CurrentDate.AddDays(1)
Loop While CurrentDate <= LastDate
'Save to Excel
'FileInfo Object
Dim FileName as String=""
'Declare SaveFileDialg
Dim sfd As New SaveFileDialog
sfd.Filter = "Excel files (*.xlsx)|*.xlsx"
If sfd.ShowDialog = DialogResult.OK Then
FileName=sfd.FileName
End If
Dim fi As New FileInfo(FileName)
Dim Ep As New ExcelPackage()
Dim wb As ExcelWorkbook=Ep.Workbook
Dim Sheet As ExcelWorksheet=wb.Worksheets.Add("EMB_Data")
Sheet.Cells(1,1).LoadFromDataTable(ResultsTable,True)
Ep.SaveAs(FI)
Here is the result. Note that if we calculate the distance based on these coordinates, we really get very very close to distance between Sun and Earth. The number is not correct as our center of Coordinate system is not the center of the Sun, and our target body is not Earth, but Earth Moon Barycenter. Moon itself does not rotate strictly around the Earth. It rotates around the Center of Mass of Earth-Moon system. Reducing the coordinates to center of the Earth is called barycentring and will be described in the next post, along with the method of calculating Astrometric coordinates of Sun and planets for a given True Barycentric Time.
Note the line 25, or better CurrentDate.ToJulianDay part. Here I used something called an Extension method. I have extended the DateTime object to calculate JulianDay number from an existing DateTime object. The implementation follows:
Imports System.Runtime.CompilerServices Public Module ExtensionMethods ''' <summary> ''' Adds a method to calculate Julian Day from an existing DateTime Strure ''' </summary> ''' <param name="mDatetIime"></param> ''' <returns></returns> <Extension()> Public Function ToJulianDay(ByVal mDatetIime As DateTime) As Decimal Dim JD As Decimal Dim year As Integer = mDatetIime.Year Dim month As Integer = mDatetIime.Month Dim day As Integer = mDatetIime.Day If month < 3 Then month += 12 year -= 1 End If Dim A As Integer = Math.Floor(year/100) Dim B As Integer = 2 - A + Math.Floor(A/4) JD = Math.Floor(365.25*(year + 4716)) + Math.Floor(30.6001*(month + 1)) + day + B - 1524.5 JD = JD + (mDatetIime.Hour + mDatetIime.Minute/60 + mDatetIime.Second/3600 + mDatetIime.Millisecond/3600000)/24 Return JD End Function End Module
At the begining I tried using a special function and using DateTime as an input parameter. Using extension method just seems more elegant way of doing it.
Wrap up
In this post, we learned a ton of things. Description of NASA JPL DE405 epherides was shown. We learned how to convert raw .405 files to a more readable one column versions either in Excel or text format. Next we used Excel files generated to read and parse Coefficients sets for a given Julian day and Celestial body. We then calculated Solar System Barycenter Rectangular Coordinates for a given Celestial Body, and Julian day. We tested it for Eart-moon Barycenter.
From Programmatic point of view, I have shown, how to use operator overloading and Extension Methods.
That's it for now. In the next post I will go deeper in different time scales used in Astronomical calculations, we will calculate Geocentric equatorial coordinates for celestial bodies, Namely Right Ascension and Declination.
Best Regards