Tutorial

Tutorial

This tutorial will walk you through the features and functionality of CALCEPH.jl

Ephemerides sources

The supported sources of ephemerides are:

Example:

download("https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de435.bsp","planets.dat")
# WARNING this next file is huge (Jupiter Moons ephemerides)
download("https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/jup310.bsp","jupiter_system.bsp")

Ephemerides context

The user first need to load the ephemerides files into an ephemerides context object that will be used later to retrieve position and velocities of celestial objects.

A context can be made from one or several files:

using CALCEPH

# load a single file in context eph1
eph1 = Ephem("planets.dat")
# load multiple files in context eph2
eph2 = Ephem(["planets.dat","jupiter_system.bsp"])

You must specify the relative or absolute path(s) of the file(s) to load.

You can prefetch the ephemerides data into main memory for faster access:

prefetch(eph2)

Epoch arguments

CALCEPH function takes the epoch as the sum of two double precision floating arguments jd1 and jd2. The sum jd1 + jd2 is interpreted as the julian date in the timescale of the ephemerides context (usually TDB or sometimes TCB).

For maximum accuracy, it is recommended to set jd2 to the fractional part of the julian date and jd1 to the difference: jd2 magnitude should be less than one while jd1 should have an integer value.

If a high accuracy in timetag is not needed, jd1 can be set to the full julian date and jd2 to zero.

Options

Many CALCEPH function takes an integer argument to store options. The value of this argument is the sum of the option to enable (each option actually corresponds to a single bit of that integer). Each option to enable can appear only once in the sum!

The following options are available:

The useNaifId option controls the identification scheme for the input arguments: target and center.

The units options controls the units of the outputs. It is compulsory to set the output units if the routine has the input argument options.

For example to compute the position and velocity in kilometers and kilometers per second of body target (given as its NAIF identification number) with respect to center (given as its NAIF identification number), the options argument should be set as such:

options = unitKM + unitSec + useNaifId

Body identification scheme

CALCEPH has the following identification scheme for bodies:

If target is 14, 15, 16 or 17 (nutation, libration, TT-TDB or TCG-TCB), center must be 0.

The more complete NAIF identification scheme can be used if the value useNaifId is added to the options argument.

NAIF body identification scheme

See https://naif.jpl.nasa.gov/pub/naif/toolkitdocs/C/req/naifids.html

CALCEPH uses this identification scheme only when the value useNaifId is added to the options argument.

The CALCEPH julia wrapper comes with the naifId object which contains the mapping between NAIF identification numbers and names:

julia> naifId.id[:sun]
10

julia> naifId.id[:mars]
499

julia> naifId.names[0]
Set(Symbol[:ssb, :solar_system_barycenter])

naifId also stores the following identifiers:

naifId is actually an instance of mutable struct BodyId. The user can also create its own identification scheme for its SPICE kernels:

const MyUniverseIds = CALCEPH.BodyId()
CALCEPH.add!(MyUniverseIds,:tatooine,1000001)
CALCEPH.add!(MyUniverseIds,:dagobah,1000002)
CALCEPH.add!(MyUniverseIds,:endor,1000003)
CALCEPH.add!(MyUniverseIds,:deathstar,1000004)
CALCEPH.add!(MyUniverseIds,:endor_deathstar_system_barycenter,1000005)
CALCEPH.add!(MyUniverseIds,:edsb,1000005)

You can also load identification data from an external file:

CALCEPH.loadData!(MyUniverseIds, "MyUniverseIds.txt")

See example: https://github.com/JuliaAstro/CALCEPH.jl/blob/master/data/NaifIds.txt

Names from the file are converted to lower case and have spaces replaced by underscores before being converted to symbols/interned strings.

Computing positions and velocities:

The following methods are available to compute position and velocity with CALCEPH:

compute(eph,jd1,jd2,target,center)
compute(eph,jd1,jd2,target,center,options)
compute(eph,jd1,jd2,target,center,options,order)

Those methods compute the position and its time derivatives of target with respect to center.

When order is not specified, position and velocity are computed.

Example:

Computing position only of Jupiter system barycenter with respect to the Earth Moon center in kilometers at JD=2456293.5 (Ephemeris Time).

options = useNaifId + unitKM + unitSec
jd1 = 2456293.0
jd2 = 0.5
center = naifId.id[:moon]
target = naifId.id[:jupiter_barycenter]
pos = compute(eph2, jd1, jd2, target, center, options,0)

Computing orientation:

The following methods are available to compute orientation angles with CALCEPH:

orient(eph,jd1,jd2,target,options)
orient(eph,jd1,jd2,target,options,order)

Those methods compute the Euler angles of target and their time derivatives.

Example:

JPL DE405 binary ephemerides contain Chebychev polynomials for the IAU 1980 nutation theory. Interpolating those is much faster than computing the IAU 1980 nutation series. Computing Earth nutation angles in radians at JD=2456293.5 (Ephemeris Time).

download("ftp://ssd.jpl.nasa.gov/pub/eph/planets/Linux/de405/lnxp1600p2200.405","DE405")
eph1 = Ephem("DE405")
options = useNaifId + unitRad + unitSec + outputNutationAngles
jd1 = 2456293.0
jd2 = 0.5
target = naifId.id[:earth]
angles = orient(eph1, jd1, jd2, target, options,0)

Note that the returned value is a vector of 3 even though there are only 2 nutation angles. The last value is zero and meaningless.

Computing angular momentum:

The following methods are available to compute body angular momentum with CALCEPH:

rotAngMom(eph,jd1,jd2,target,options)
rotAngMom(eph,jd1,jd2,target,options,order)

Those methods compute the angular momentum of target and their time derivatives.

Time ephemeris

The time ephemeris TT-TDB or TCG-TCB at the geocenter can be evaluated with a suitable source.

INPOP and some JPL DE ephemerides includes a numerically integrated time ephemeris for the geocenter which is usually more accurate than the analytical series: Moreover it is much faster to interpolate those ephemerides than to evaluate the analytical series. This is only for the geocenter but a simple correction can also be added for the location of the observer (and its velocity in case the observer is on a highly elliptical orbit).

Files that can be used to obtain the difference between TT and TDB are, e.g.:

Example:

Computing TT-TDB at geocenter in seconds at JD=2456293.5 (Ephemeris Time).

download("ftp://ftp.imcce.fr/pub/ephem/planets/inpop17a/inpop17a_TDB_m100_p100_tt.dat","INPOP17a")
eph1 = Ephem("INPOP17a")
options = useNaifId + unitSec
jd1 = 2456293.0
jd2 = 0.5
target = naifId.id[:ttmtdb]
center = naifId.id[:timecenter]
ttmtdb = compute(eph1, jd1, jd2, target, center, options,0)

Note that the returned value is a vector of 3 even though there is only one meaningful value. The last 2 values are zero and meaningless.

In place methods

In place versions of the methods described above are also available. Those are:

unsafe_compute!(result,eph,jd1,jd2,target,center)
unsafe_compute!(result,eph,jd1,jd2,target,center,options)
unsafe_compute!(result,eph,jd1,jd2,target,center,options,order)
unsafe_orient!(result,eph,jd1,jd2,target,options)
unsafe_orient!(result,eph,jd1,jd2,target,options,order)
unsafe_rotAngMom!(result,eph,jd1,jd2,target,options)
unsafe_rotAngMom!(result,eph,jd1,jd2,target,options,order)

Those methods do not perform any checks on their inputs. In particular, result must be a contiguous vector of double precision floating point number of dimension at least 6 when order is not specified or at least 3*(order+1) otherwise.

Constants

Ephemerides files may contain related constants. Those can be obtained by the constants method which returns a dictionary:

download("ftp://ftp.imcce.fr/pub/ephem/planets/inpop17a/inpop17a_TDB_m100_p100_tt.dat","INPOP17a")
eph1 = Ephem("INPOP17a")
# retrieve constants from ephemeris as a dictionary
con = constants(eph1)
# list the constants
keys(con)
# get the sun J2
J2sun = con[:J2SUN]

Introspection

Time scale

timeScale(eph)

returns the Ephemeris Time identifier:

Time span

timespan(eph)

returns the triplet:

Position records

positionRecords(eph)

retrieve position records metadata in ephemeris associated to handler eph. This is a vector of metadata about the ephemerides records ordered by priority. The compute methods use the highest priority ephemerides records when there are multiple records that could satisfy the target and epoch.

Each record metadata contains the following information:

Orientation records

orientationRecords(eph)

retrieve orientation records metadata in ephemeris associated to handler eph. This is a vector of metadata about the ephemerides records ordered by priority. The orient methods use the highest priority ephemerides records when there are multiple records that could satisfy the target and epoch.

Each record metadata contains the following information:

Cleaning up

Because, Julia's garbage collector is lazy, you may want to free the memory managed by the context before you get rid of the reference to the context with eg:

finalize(eph1)
eph1 = Nothing

or after with

eph1 = Nothing
GC.gc()

Error handling

By default, the CALCEPH C library prints error messages directly to the standard output but this can be modified.

The Julia wrapper provides the following interface for this purpose:

CALCEPH.setCustomHandler(f)

where f should be a user function taking a single argument of type String which will contain the CALCEPH error message. f should return Nothing.

To disable CALCEPH error messages printout to the console:

CALCEPH.setCustomHandler(s->Nothing)

To get back the default behavior:

CALCEPH.disableCustomHandler()