I have mentioned various kinds of configuration file formats used in Fortran here before. One that I haven't mentioned is the text PCK file format used by the NAIF SPICE Toolkit. This is a format that is similar in some ways to Fortran namelists, but with a better API for reading the file and querying the contents. Variables read from a PCK file are inserted into the SPICE variable pool. An example PCK file (pck00010.tpc) can be found here.
The PCK file format has some limitations. It doesn't allow inline comments (only block comments separate from the variable declarations). Integers are stored as rounded doubles, and logical variables are not supported at all (you have to use 0.0 or 1.0 for this). One particular annoyance is that the files are not cross platform (the line breaks must match the platform they are being used on). However, SPICELIB also provides routines to programmatically enter variables into the pool. This allows us to create files in other formats that can be used in a SPICE application. For example, we can use JSON-Fortran to read variables in JSON files and insert them into the SPICE pool.
First, let's declare the interfaces to the SPICE routines we need (since SPICELIB is straight up Fortran 77, they are not in a module):
interface
subroutine pcpool ( name, n, cvals )
implicit none
character(len=*) :: name
integer :: n
character(len=*) :: cvals ( * )
end subroutine pcpool
subroutine pdpool ( name, n, values )
import :: wp
implicit none
character(len=*) :: name
integer :: n
real(wp) :: values ( * )
end subroutine pdpool
subroutine pipool ( name, n, ivals )
implicit none
character(len=*) :: name
integer :: n
integer :: ivals ( * )
end subroutine pipool
subroutine setmsg ( msg )
implicit none
character(len=*) :: msg
end subroutine setmsg
subroutine sigerr ( msg )
implicit none
character(len=*) :: msg
end subroutine sigerr
end interface
Then, we can use the following routine:
subroutine json_to_spice(jsonfile)
implicit none
character(len=*),intent(in) :: jsonfile
!! the JSON file to load
integer :: i,j
type(json_core) :: json
type(json_value),pointer :: p,p_var,p_element
real(wp) :: real_val
integer :: var_type,int_val,n_vars,n_children,element_var_type
logical :: found
integer,dimension(:),allocatable :: int_vec
real(wp),dimension(:),allocatable :: real_vec
character(len=:),dimension(:),allocatable :: char_vec
character(len=:),allocatable :: char_val,name
integer,dimension(:),allocatable :: ilen
nullify(p)
! allow for // style comments:
call json%initialize(comment_char='/')
! read the file:
call json%parse(file=jsonfile, p=p)
if (json%failed()) then
! we will use the SPICE error handler:
call setmsg ( 'json_to_spice: Could not load file: '&
trim(jsonfile) )
call sigerr ( 'SPICE(INVALIDJSONFILE)' )
else
! how many variables are in the file:
call json%info(p,n_children=n_vars)
main : do i = 1, n_vars
call json%get_child(p, i, p_var, found)
! what kind of variable is it?:
call json%info(p_var,var_type=var_type,name=name)
if (var_type == json_array) then
! how many elements in this array?:
call json%info(p_var,n_children=n_children)
! first make sure all the variables are the same type
! [must be integer, real, or character]
do j = 1, n_children
call json%get_child(p_var, j, p_element, found)
call json%info(p_element,var_type=element_var_type)
if (j==1) then
var_type = element_var_type
else
if (var_type /= element_var_type) then
call setmsg ( 'json_to_spice: Invalid array ('&
trim(name)//') in file: '//trim(jsonfile) )
call sigerr ( 'SPICE(INVALIDJSONVAR)' )
exit main
end if
end if
end do
! now we know the var type, so get as a vector:
select case (var_type)
case(json_integer); call json%get(p_var,int_vec)
case(json_double ); call json%get(p_var,real_vec)
case(json_string ); call json%get(p_var,char_vec,ilen)
case default
call setmsg ( 'json_to_spice: Invalid array ('&
trim(name)//') in file: '//trim(jsonfile) )
call sigerr ( 'SPICE(INVALIDJSONVAR)' )
exit main
end select
else
! scalar:
n_children = 1
select case (var_type)
case(json_integer)
call json%get(p_var,int_val); int_vec = [int_val]
case(json_double )
call json%get(p_var,real_val); real_vec = [real_val]
case(json_string )
call json%get(p_var,char_val); char_vec = [char_val]
case default
call setmsg ( 'json_to_spice: Invalid variable ('&
trim(name)//') in file: '//trim(jsonfile) )
call sigerr ( 'SPICE(INVALIDJSONVAR)' )
exit main
end select
end if
! now, add them to the pool:
select case (var_type)
case(json_integer)
call pipool ( name, n_children, int_vec )
case(json_double )
call pdpool ( name, n_children, real_vec )
case(json_string )
call pcpool ( name, n_children, char_vec )
end select
end do main
end if
call json%destroy(p)
end subroutine json_to_spice
Here we are using the SPICE routines PIPOOL, PDPOOL, and PCPOOL to insert variables into the pool. We are also using the built-in SPICE error handling routines (SETMSG and SIGERR) to manage errors. The code works for scalar and array variables. It can be used to read the following example JSON file:
{
"str": "qwerty", // a scalar string
"strings": ["a", "ab", "abc", "d"], // a vector of strings
"i": 8, // a scalar integer
"ints": [255,127,0], // a vector of integers
"r": 123.345, // a scalar real
"reals": [999.345, 5671.8] // a vector of reals
}
After reading it, we can verify that the variables were added to the pool by printing the pool contents using the routine WRPOOL. This produces:
\begindata
str = 'qwerty'
i = 0.80000000000000000D+01
strings = ( 'a',
'ab',
'abc',
'd' )
ints = ( 0.25500000000000000D+03,
0.12700000000000000D+03,
0.00000000000000000D+00 )
r = 0.12334500000000000D+03
reals = ( 0.99934500000000003D+03,
0.56718000000000002D+04 )
\begintext
See also
JPL's SPICE Toolkit (SPICELIB) is the premier software library for computations related to solar system geometry. It is freely distributed, and is also one of the best-documented libraries I have ever come across. SPICELIB also includes a comprehensive set of routines for date and time conversions. An example is shown here:
program spice_test
use iso_fortran_env, only: wp => real64
implicit none
interface
! SPICELIB routines
subroutine timout ( et, pictur, output )
import :: wp
implicit none
real(wp),intent(in) :: et
character(len=*),intent(in) :: pictur
character(len=*),intent(out) :: output
end subroutine timout
subroutine str2et ( string, et )
import :: wp
implicit none
character(len=*),intent(in) :: string
real(wp),intent(out) :: et
end subroutine str2et
subroutine furnsh ( file )
implicit none
character(len=*),intent(in) :: file
end subroutine furnsh
end interface
character(len=*),parameter :: time_in = &
'2017 Aug 12 00:00:00 TDB'
character(len=*),parameter :: pictur = &
'Mon DD,YYYY HR:MN:SC.#### UTC ::UTC'
real(wp) :: et
character(len=100) :: time_out
! load the leap second kernel:
call furnsh('naif0012.tls')
! example conversion:
call str2et(time_in, et)
call timout(et, pictur, time_out)
write(*,*) 'time_in: ', time_in
write(*,*) 'et: ', et
write(*,*) 'time_out: ', time_out
end program spice_test
A few things to note:
- Here we are using the SPICE routines str2et and timout to convert a string from a TDB calendar date to ephemeris time and then to a UTC calendar date. These routines are very flexible and can convert a wide range of date formats. Other routines are available to do other transformations.
- The base time system of SPICE is Barycentric Dynamical Time (TDB). "Ephemeris time" is a count of TDB seconds since the J2000 epoch (Jan 1, 2000 12:00:00).
- We have to load the latest leap second kernel (naif0012.tls in this case), which is necessary to define UTC.
- The SPICE routines are not in a module (the code is Fortran 77), and so have no explicit interfaces. Thus it is good practice to specify them as I do here.
The output of this example is:
time_in: 2017 Aug 12 00:00:00 TDB
et: 555768000.00000000
time_out: Aug 11,2017 23:58:50.8169 UTC
See also
JPL recently released an update to their awesome SPICE Toolkit (it is now at version N66). The major new feature in this release is the Digital Shape Kernel (DSK) capability to define the shapes of bodies (such as asteroids) via tessellated plate models.
Unfortunately for Fortran users, they also announced that they have decided to reimplement the entire library in C++. SPICELIB is currently written in Fortran 77, which they f2c to provide a C version (which is also callable from IDL, Matlab, and Python, among others). Their reason for this "upgrade" is to provide thread safety and object oriented features. Of course, modern Fortran can be thread safe and object oriented, and upgrading the code to modern standards could be done in a fraction of the time it will take to rewrite everything from scratch in C++. SPICELIB is extremely well-written Fortran 77 code, and is not infested with COMMON blocks, EQUIVALENCE statements, etc. I actually don't think it would take much effort to modernize it. In addition, Fortran/C interoperability could be employed to easily provide an interface that is callable from C without source transformation.
However, I guess it isn't meant to be, and the science/engineering community will lose another Fortran code to C++ like many times before, in spite of C++ being a terrible language for scientists and engineers.
See also
The IAU Working Group on Cartographic Coordinates and Rotational Elements (WGCCRE) is the keeper of official models that describe the cartographic coordinates and rotational elements of planetary bodies (such as the Earth, satellites, minor planets, and comets). Periodically, they release a report containing the coefficients to compute body orientations, based on the latest data available. These coefficients allow one to compute the rotation matrix from the ICRF frame to a body-fixed frame (for example, IAU Earth) by giving the direction of the pole vector and the prime meridian location as functions of time. An example Fortran module illustrating this for the IAU Earth frame is given below. The coefficients are taken from the 2009 IAU report [Reference 1]. Note that the IAU models are also available in SPICE (as "IAU_EARTH", "IAU_MOON", etc.). For Earth, the IAU model is not suitable for use in applications that require the highest possible accuracy (for that a more complex model would be necessary), but is quite acceptable for many applications.
module iau_orientation_module
use, intrinsic :: iso_fortran_env, only: wp => real64
implicit none
private
!constants:
real(wp),parameter :: zero = 0.0_wp
real(wp),parameter :: one = 1.0_wp
real(wp),parameter :: pi = acos(-one)
real(wp),parameter :: pi2 = pi/2.0_wp
real(wp),parameter :: deg2rad = pi/180.0_wp
real(wp),parameter :: sec2day = one/86400.0_wp
real(wp),parameter :: sec2century = one/3155760000.0_wp
public :: icrf_to_iau_earth
contains
!Rotation matrix from ICRF to IAU_EARTH
pure function icrf_to_iau_earth(et) result(rotmat)
implicit none
real(wp),intent(in) :: et !ephemeris time [sec from J2000]
real(wp),dimension(3,3) :: rotmat !rotation matrix
real(wp) :: w,dec,ra,d,t
real(wp),dimension(3,3) :: tmp
d = et * sec2day !days from J2000
t = et * sec2century !Julian centuries from J2000
ra = ( - 0.641_wp * t ) * deg2rad
dec = ( 90.0_wp - 0.557_wp * t ) * deg2rad
w = ( 190.147_wp + 360.9856235_wp * d ) * deg2rad
!it is a 3-1-3 rotation:
tmp = matmul( rotmat_x(pi2-dec), rotmat_z(pi2+ra) )
rotmat = matmul( rotmat_z(w), tmp )
end function icrf_to_iau_earth
!The 3x3 rotation matrix for a rotation about the x-axis.
pure function rotmat_x(angle) result(rotmat)
implicit none
real(wp),dimension(3,3) :: rotmat !rotation matrix
real(wp),intent(in) :: angle !angle [rad]
real(wp) :: c,s
c = cos(angle)
s = sin(angle)
rotmat = reshape([one, zero, zero, &
zero, c, -s, &
zero, s, c],[3,3])
end function rotmat_x
!The 3x3 rotation matrix for a rotation about the z-axis.
pure function rotmat_z(angle) result(rotmat)
implicit none
real(wp),dimension(3,3) :: rotmat !rotation matrix
real(wp),intent(in) :: angle !angle [rad]
real(wp) :: c,s
c = cos(angle)
s = sin(angle)
rotmat = reshape([ c, -s, zero, &
s, c, zero, &
zero, zero, one ],[3,3])
end function rotmat_z
end module iau_orientation_module
References
- B. A. Archinal, et al., "Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009", Celest Mech Dyn Astr (2011) 109:101-135.
- J. Williams, Fortran Astrodynamics Toolkit - iau_orientation_module [GitHub]
The SPICE Toolkit software is an excellent package of very well-written and well-documented routines for a variety of astrodynamics applications. It is produced by NASA's Navigation and Ancillary Information Facility (NAIF). Versions are available for Fortran 77, C, IDL, and Matlab.
To speed up the execution of SPICE-based programs, there are a few things you can do:
- For the Fortran SPICELIB, recompile it with optimization enabled (say, -O2). The default library released by NAIF is not compiled with optimization.
- Turn off the SPICE traceback system (
call TRCOFF()). If you are using a compiler that has a built-in stack trace routine (for example TRACEBACKQQ in the Intel compiler), just include a call to it in the SPICE BYEBYE.F routine. That will give you a stack trace for any fatal errors.
- Converting a non-native binary PCK to native form will also speed up data access somewhat.
- Calls to ephemeris routines where the target and observer bodies are input as strings will be slightly slower than the ones where the inputs are the NAIF ID codes (which are integers). For example, for the geometric state (position and velocity) of one body relative to another, use SPKGEO instead of SPKEZR. Also, if all you require is position and not velocity, use SPKGPS instead of SPKGEO, since it will be a bit faster.
Also note that for applications only requiring the ephemerides of the solar system major bodies, there is an older code from JPL (PLEPH) which is simpler and faster than SPICE, but uses a different format for the ephemeris files.
JPL just released a new version of the SPICE Toolkit (N65). According to the What's New page, changes include:
- support for some new environments and termination of some old environments
- new Geometry Finder (GF) interfaces -- illumination angle search, phase angle search, user-defined binary quantity search
- new high-level SPK APIs allowing one to specify the observer or target as a location with a constant position and velocity rather than as an ephemeris object
- new high-level APIs that check for occultation and in-Field-Of-View (FOV) conditions
- new high-level illumination angle routines
- new high-level reference frame transformation routine
- new high-level coordinate system transformation routine for states
- many Icy and Mice APIs that were formerly available only in SPICELIB and CSPICE
- new data types: SPK (19, 20, 21), PCK (20), and CK (6)
- performance improvements in the range of 10-50 percent for some types of use
- ability to load up to 5000 kernels
- increased buffers in the POOL and in SPK and CK segment search subsystems
- new built-in body name/ID code mappings and body-fixed reference frames
- a significant upgrade of utility SPKDIFF including the ability to sample SPK data
- updates to other toolkit utilities -- BRIEF, CKBRIEF, FRMDIFF, MKSPK, MSOPCK
- bug fixes in Toolkit and utility programs
