jpl_ephemeris_module Module

For reading the JPL planetary and lunar ephemerides. This is an extensively modified version of the original FORTRAN 77 code from JPL.

Ephemeris files

Note that this module uses the JPL binary ephemeris files, which can be obtained using the instructions here. See also the comments in ephemeris_test for more details.

License

Original JPL License

THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE SOFTWARE AND RELATED MATERIALS, HOWEVER USED.

IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.

RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.

Modifications

Modifications for the Fortran Astrodynamics Toolkit are covered under the following license.

History

  • Original code from JPL Version : March 25, 2013
  • Extensive modifications by Jacob Williams for the Fortran Astrodynamics Toolkit.

Note

Warning: all calculations here are done with 64 bit reals. if using 128-bit reals the state is just copied into a 128-bit variable. Perhaps we should do the interpolations in 128-bit? (the data from the file must be read as 64-bit)


Uses

  • module~~jpl_ephemeris_module~~UsesGraph module~jpl_ephemeris_module jpl_ephemeris_module iso_fortran_env iso_fortran_env module~jpl_ephemeris_module->iso_fortran_env module~ephemeris_module ephemeris_module module~jpl_ephemeris_module->module~ephemeris_module module~kind_module kind_module module~jpl_ephemeris_module->module~kind_module module~ephemeris_module->module~kind_module module~celestial_body_module celestial_body_module module~ephemeris_module->module~celestial_body_module module~kind_module->iso_fortran_env module~celestial_body_module->module~kind_module module~base_class_module base_class_module module~celestial_body_module->module~base_class_module module~numbers_module numbers_module module~celestial_body_module->module~numbers_module module~numbers_module->module~kind_module

Used by

  • module~~jpl_ephemeris_module~~UsedByGraph module~jpl_ephemeris_module jpl_ephemeris_module module~fortran_astrodynamics_toolkit fortran_astrodynamics_toolkit module~fortran_astrodynamics_toolkit->module~jpl_ephemeris_module proc~transformation_module_test transformation_module::transformation_module_test proc~transformation_module_test->module~jpl_ephemeris_module

Variables

Type Visibility Attributes Name Initial
integer, private, parameter :: wp = real64

using double precision reals since that is how the ephemeris is stored.
integer, public, parameter :: nmax = 1000

Current maximum number of ephemeris constants used in the integration and listed in the header.xxx file for the ephemeris.
integer, private, parameter :: oldmax = 400

For earlier versions of the code, the maximum number of ephemeris constants used in the integration and listed in the header.xxx file for the ephemeris.
integer, private, parameter :: nrecl = 4

nrecl=1 if recl in the open statement is the record length in s.p. words nrecl=4 if recl in the open statement is the record length in bytes
character(len=*), private, parameter, dimension(15) :: list_of_bodies = ['mercury                ', 'venus                  ', 'earth                  ', 'mars                   ', 'jupiter                ', 'saturn                 ', 'uranus                 ', 'neptune                ', 'pluto                  ', 'moon                   ', 'sun                    ', 'solar-system barycenter', 'earth-moon barycenter  ', 'nutations              ', 'librations             ']

These correspond to the numbering convention for 'ntarg' and 'ncent' in get_state:

Derived Types

type, public, extends(ephemeris_class) ::  jpl_ephemeris

Main class for accessing a JPL ephemeris file.

Components

Type Visibility Attributes Name Initial
character(len=:), public, allocatable :: namfil

name of the binary ephemeris file
integer, public :: ksize = 2036

ksize must be set by the user according to the ephemeris to be read: for de200, set ksize=1652, for de405, set ksize=2036, for de406, set ksize=1456, for de414, set ksize=2036, for de418, set ksize=2036, for de421, set ksize=2036, for de422, set ksize=2036, for de423, set ksize=2036, for de424, set ksize=2036, for de430, set ksize=2036.
integer, public, dimension(3,13) :: ipt = 0

ipt(39)
real(kind=wp), public, dimension(nmax) :: cval = 0.0_wp
real(kind=wp), public, dimension(3) :: ss = 0.0_wp
real(kind=wp), public :: au = 0.0_wp
real(kind=wp), public :: emrat = 0.0_wp
integer, public :: ncon = 0
integer, public :: numde = 0
real(kind=wp), public, dimension(6) :: pvsun = 0.0_wp

dp 6-word array containing the barycentric position and velocity of the sun.
logical, public :: km = .true.

logical flag defining physical units of the output states. km = .true. : km and km/sec km = .false. : au and au/day for nutations and librations. angle unit is always radians.
logical, public :: bary = .false.

logical flag defining output center. only the 9 planets are affected. bary = .true. : center is solar-system barycenter bary = .false. : center is sun
character(len=6), public, dimension(14,3) :: ttl = ''
character(len=6), public, dimension(nmax) :: cnam = ''
logical, public :: initialized = .false.

is the ephemeris initialized?
integer, public :: nrfile = 0

file unit for the ephemeris file
integer, public :: nrl = -1

this was formerly in state
integer, public :: ncoeffs = 0
real(kind=wp), public, dimension(1500) :: buf = 0.0_wp
real(kind=wp), public, dimension(18) :: pc = 0.0_wp
real(kind=wp), public, dimension(18) :: vc = 0.0_wp
integer, public :: np = 2
integer, public :: nv = 3
real(kind=wp), public :: twot = 0.0_wp

Type-Bound Procedures

procedure, public :: get_rv => get_rv_from_jpl_ephemeris
procedure, public :: initialize => initialize_ephemeris
procedure, public :: get_state
procedure, public :: get_constants
procedure, public :: close => close_ephemeris
procedure, public :: interp
procedure, public :: state

Functions

private pure function spice_id_to_old_id(spice_id) result(old_id)

Author
Jacob Williams
Date
3/20/2016

Convert the NAIF SPICE ID code to the old one used by the JPL ephemeris. Returns 0 if the body was not found.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: spice_id

the ID code used by SPICE

Return Value integer

the ID code used by this module (old JPL ephemeris code)


Subroutines

private subroutine get_rv_from_jpl_ephemeris(me, et, targ, obs, rv, status_ok)

Interface for the ephemeris_module.

Arguments

Type IntentOptional Attributes Name
class(jpl_ephemeris), intent(inout) :: me
real(kind=fat_wp), intent(in) :: et

ephemeris time [sec]
type(celestial_body), intent(in) :: targ

target body
type(celestial_body), intent(in) :: obs

observer body
real(kind=fat_wp), intent(out), dimension(6) :: rv

state of targ w.r.t. obs [km,km/s] in ICRF frame
logical, intent(out) :: status_ok

true if there were no problems

private subroutine get_state(me, jd, ntarg, ncent, rrd, status_ok)

This subroutine reads the JPL planetary ephemeris and gives the position and velocity of the point ntarg with respect to ncent.

Read more…

Arguments

Type IntentOptional Attributes Name
class(jpl_ephemeris), intent(inout) :: me
real(kind=wp), intent(in) :: jd

d.p. Julian ephemeris date at which interpolation is wanted.
integer, intent(in) :: ntarg

integer number of 'target' point.
integer, intent(in) :: ncent

integer number of 'center' point.
real(kind=wp), intent(out), dimension(6) :: rrd

output 6-word d.p. array containing position and velocity of point ntarg relative to ncent. the units are AU and AU/day (or km and km/sec if me%km=.true.). For librations the units are radians and radians per day. In the case of nutations the first four words of rrd will be set to nutations and rates, having units of radians and radians/day.
logical, intent(out) :: status_ok

true if there were no problems

private subroutine interp(me, buf, t, ncf, ncm, na, ifl, pv)

this subroutine differentiates and interpolates a set of chebyshev coefficients to give position and velocity.

Arguments

Type IntentOptional Attributes Name
class(jpl_ephemeris), intent(inout) :: me
real(kind=wp), dimension(ncf,ncm,*) :: buf

1st location of array of d.p. chebyshev coefficients of position
real(kind=wp), intent(in), dimension(2) :: t
integer, intent(in) :: ncf Read more…
integer, intent(in) :: ncm Read more…
integer, intent(in) :: na

(i.e., # of sub-intervals in full interval)

Read more…
integer, intent(in) :: ifl

integer flag = 1 for positions only = 2 for pos and vel
real(kind=wp), dimension(ncm,*) :: pv

interpolated quantities requested. dimension expected is pv(ncm,ifl), dp.

private subroutine split(tt, fr)

this subroutine breaks a d.p. number into a d.p. integer and a d.p. fractional part.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: tt

d.p. input number
real(kind=wp), intent(out), dimension(2) :: fr

d.p. 2-word output array. fr(1) contains integer part fr(2) contains fractional part for negative input numbers, fr(1) contains the next more negative integer; fr(2) contains a positive fraction.

private subroutine initialize_ephemeris(me, filename, ksize, km, bary, status_ok)

Initialize the ephemeris. This routine may be called to load a different ephemeris file. Otherwise, it is called on the first call to get_state, and loads the file specified in the module header.

Read more…

Arguments

Type IntentOptional Attributes Name
class(jpl_ephemeris), intent(inout) :: me
character(len=*), intent(in) :: filename

ephemeris file name
integer, intent(in), optional :: ksize

corresponding ksize
logical, intent(in), optional :: km

defining physical units of the output states. km = .true. : km and km/sec [default], km = .false. : au and au/day.
logical, intent(in), optional :: bary

logical flag defining output center. only the 9 planets are affected. bary = .true. : center is solar-system barycenter, bary = .false. : center is sun [default].
logical, intent(out) :: status_ok

true if there were not problems.

private subroutine close_ephemeris(me)

Close the ephemeris.

Arguments

Type IntentOptional Attributes Name
class(jpl_ephemeris), intent(inout) :: me

private subroutine state(me, et2, list, pv, pnut, status_ok)

This subroutine reads and interpolates the JPL planetary ephemeris file.

Read more…

Arguments

Type IntentOptional Attributes Name
class(jpl_ephemeris), intent(inout) :: me
real(kind=wp), intent(in), dimension(2) :: et2

2-word Julian ephemeris epoch at which interpolation is wanted. any combination of et2(1)+et2(2) which falls within the time span on the file is a permissible epoch.

Read more…
integer, intent(in), dimension(12) :: list

12-word integer array specifying what interpolation is wanted for each of the bodies on the file. list(i) = 0 : no interpolation for body i, list(i) = 1 : position only, list(i) = 2 : position and velocity.

Read more…
real(kind=wp), intent(out), dimension(6,11) :: pv

dp 6 x 11 array that will contain requested interpolated quantities (other than nutation, stored in pnut). the body specified by list(i) will have its state in the array starting at pv(1,i). (on any given call, only those words in pv which are affected by the first 10 list entries, and by list(12) if librations are on the file, are set. the rest of the pv array is untouched.) the order of components starting in pv(1,i) is: x,y,z,dx,dy,dz.

Read more…
real(kind=wp), intent(out), dimension(4) :: pnut

dp 4-word array that will contain nutations and rates, depending on the setting of list(11). the order of quantities in pnut is:

Read more…
logical, intent(out) :: status_ok

true if there were no problems

private subroutine get_constants(me, nam, val, sss, n)

Obtain the constants from the ephemeris file.

Arguments

Type IntentOptional Attributes Name
class(jpl_ephemeris), intent(inout) :: me
character(len=6), intent(out), dimension(:) :: nam

array of constant names
real(kind=wp), intent(out), dimension(:) :: val

array of values of constants
real(kind=wp), intent(out), dimension(3) :: sss

jd start, jd stop, step of ephemeris
integer, intent(out) :: n

number of entries in nam and val arrays

public subroutine ephemeris_test()

Ephemeris test routine.

Read more…

Arguments

None