Blog

Slovenian Nautical Almanac - Part 2 - JPL DE 405 Ephemerides

22.02.2016
Aljoša Kocen

In the process of obtaining the data for the Ephemerides calculation, I really only had two options available:

  1. 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.
  2. 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
    • 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
    • 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:

  1. Open ascpXXXX.405 file (XXXX being the "year" in the file name)
  2. Declare auxiliarry list, which will hold the data. I used System.Data.DataTable object. 
  3. Loop through the ASCII file line by line and:
    1. Trim excessive spaces from the line
    2. split the line into array of data
    3. If array has only two records, and the second one equals "1018", skip complete line
    4. In any other case, we have three numbers, each of which is inserted into the DataTable
  4. With the aid of EPPlus Excel library, create an empty Excel Workbook
  5. Load the resulting DataTable into an empty Excel Worksheet
  6. 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:

  1. 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
  2. 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.

  3. 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.
  4. 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
 
Now we simply create a list of these "Bodies" (In quotation marks as Nutations and Librations are really not Celestial bodies):
 
 
 
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
 
We need to declare a special object which holds the Ephemerides Coefficients: We shall name it 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
 
EphemDataBlock inherits the 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.
 
The routine which reads the data from the Excel file  is listed Next. It accepts Julian Day (A single decimal value) as input and returns full dataset from a single Excel file into the 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
 
We are using Epplus library here again. Why did we load entire Set (some 240.000 records) and not only a single 32 day EphemDataBlock? We will see later in the series, that calculating celestial positions for the range of dates (Almanacs tipically list Ephemerides for an entire year), would otherwise mean that we need to open the data file each time we cross 32 day border for each EphemDataBlock, which would result in massive performance penalty.
 
Given the basic DataBlock one can calculate  Cartesian coordinates as follows. First we create an object to hold our X,Y, and Z Coordinates. We shall call it 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 Class
Note 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
Overloading basic operators in this case makes our code much more readable and much less error-prone. I've included also overloaded operators for scalar multiplication and division. It will come handy while doing the barycentric for positions of Earth. EPhemerides are given for 
 
 
To Obtain Cartesian coordinates for a given celestial body, we create a class called 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
The function GetBodyPosition accepts CelestialBody object (described above) Julian day as an input and then:
 
  1. It loads the DataBlock with 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. 
  2. Depending on the  JulianDay input, it extracts the 1020 record SubBlock
  3. 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
  4. loads three Coefficient arrays with the data from the  SubBlock (lines 50 to 71)
  5. Calculates the Chebyshev time, which is always in the interval from -1 to 1, and uses it to get the Chebyshev polynomial terms
  6. 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