Libration Points in the Earth-Moon System. From Farquhar, 1971.
The Circular Restricted Three-Body Problem (CR3BP) is a continuing source of delight for astrodynamicists. If you think that, after 200 years, everything that could ever be written about this subject has already been written, you would be wrong. The basic premise is to describe the motion of a spacecraft in the vicinity of two celestial bodies that are in circular orbits about their common barycenter. The "restricted" part means that the spacecraft does not affect the motion of the other two bodies (i.e., its mass is negligible).
The normalized equations of motion for the CR3BP are very simple (see the Fortran Astrodynamics Toolkit):
subroutine crtbp_derivs(mu,x,dx)
implicit none
real(wp),intent(in) :: mu
!! CRTBP parameter
real(wp),dimension(6),intent(in) :: x
!! normalized state [r,v]
real(wp),dimension(6),intent(out) :: dx
!! normalized state derivative [rdot,vdot]
real(wp),dimension(3) :: r1,r2,rb1,rb2,r,v,g
real(wp) :: r13,r23,omm,c1,c2
r = x(1:3) ! position
v = x(4:6) ! velocity
omm = one - mu
rb1 = [-mu,zero,zero] ! location of body 1
rb2 = [omm,zero,zero] ! location of body 2
r1 = r - rb1 ! body1 -> sc vector
r2 = r - rb2 ! body2 -> sc vector
r13 = norm2(r1)**3
r23 = norm2(r2)**3
c1 = omm/r13
c2 = mu/r23
!normalized gravity from both bodies:
g(1) = -c1*(r(1) + mu) - c2*(r(1)-one+mu)
g(2) = -c1*r(2) - c2*r(2)
g(3) = -c1*r(3) - c2*r(3)
! derivative of x:
dx(1:3) = v ! rdot
dx(4) = two*v(2) + r(1) + g(1) ! vdot
dx(5) = -two*v(1) + r(2) + g(2) !
dx(6) = g(3) !
end subroutine crtbp_derivs
There are many types of interesting periodic orbits that exist in this system. I've mentioned Distant Retrograde Orbits (DROs) here before. Probably the most well-known CR3BP periodic orbits are called "halo orbits". Halo orbits can be computed in a similar manner as DROs (see previous post). For halos, we can get a good initial guess using an analytic approximation, and then use this guess to target the real thing. See the halo orbit module in the Fortran Astrodynamics Toolkit for an example of this approach. We can also use a continuation type method to generate families of orbits (say, for the L2 libration point, we have a range of halo orbits from near-planars that "orbit" L2, to near-rectilinears with close approaches to the secondary body). In a realistic force model (with the real ephemeris of the Earth, Moon, and Sun), the orbits generated using CR3BP assumptions are not periodic, but we can use the CR3BP solution as an initial guess in a multiple-shooting type approach with a numerical nonlinear equation solver to generate quasi-periodic multi-rev solutions. See the references below for details.
I generated the following image of a bunch of planar periodic CR3BP orbits using a program incorporating the Fortran Astrodynamics Toolkit and OpenFrames (an interface to the open source OpenSceneGraph library that can be used to visualize spacecraft trajectories). The initial states for the various orbits were taken from this database. Eventually, I will make this program open source when it is finished. There's also an online CR3BP propagator here.
See also
- M. Henon, "Numerical Exploration of the Restricted Problem. V. Hill’s Case: Periodic Orbits and Their Stability," Astronomy and Astrophys. Vol 1, 223-238, 1969.
- R. W. Farquhar, "The Utilization of Halo Orbits in Advanced Lunar Operation", NASA TN-D-6365, July 1971.
- D. L. Richardson, "Analytic Construction of Periodic Orbits About the Collinear Points", Celestial Mechanics 22 (1980)
- R. Mathur, C. A. Ocampo, "An Architecture for Incorporating Interactive Visualizations Into Scientific Simulations", 57th International Astronautical Congress, Valencia, Spain. IAC-06-D1.P.1.6.
- R. L. Restrepo, R. P. Russell, "A Database of Planar Axi-Symmetric Periodic Orbits for the Solar System," Paper AAS 17-694, AAS/AIAA Astrodynamics Specialist Conference, Stevenson, WA, Aug 2017.
- J. Williams, D. E. Lee, R. J. Whitley, K. A. Bokelmann, D. C. Davis, and C. F. Berry. "Targeting cislunar near rectilinear halo orbits for human space exploration", AAS 17-267
- Nabla Zero Labs, The Circular Restricted Three Body Problem.
- D. Davis, Near Rectilinear Halo Orbits Explained and Visualized, a.i. solutions, Jul 28, 2017 [YouTube]
There are various ways to interface Fortran and Python. However, many of the online examples on this subject seem to be content with a 40 year old variant of the language, and do not tend to focus on modern Fortran. One particular thing that is hard to figure out is how to send arbitrary or dynamic data from Python to Fortran (and vice versa)? This could be variable-sized arrays, variable-sized strings (or even a variable number of variable sized strings), or just any random data structure that contains different types of variables that aren't necessarily known at compile time. It is not easy to find out how to do this, but one option presented here is to use JSON as an intermediate format.
Sending a dict from Python to Fortran
Say we want to communicate an arbitrary Python dictionary structure to Fortran. We can do this fairly easily by dumping the dict to a JSON string, sending the string to Fortran, and then parsing the string in Fortran using JSON-Fortran. Then we have access to all the data in the dictionary in the Fortran code. An example of this is presented here.
First we need a function python_str_to_fortran() to convert a Python dict object into a c_char_p string that can be sent via the C interface to a Fortran routine:
from ctypes import *
import json
def python_str_to_fortran(s):
return cast(create_string_buffer(s.encode()),c_char_p)
def python_dict_to_fortran(d):
return python_str_to_fortran(json.dumps(d))
Note that we are using Python's ctypes module here to create a c_char_p string that can be passed to Fortran. I admit, this bit of code is the result of some trial and error to figure out what would work (the documentation is less than helpful and of course doesn't mention Fortran at all). Ctypes requires that the Fortran code be compiled as a shared library (a .DLL on Windows, or .so on Linux). We can then define the interfaces to the routines in the DLL we want to call from Python. As long as we get the types of the arguments right (which can be tricky for strings as will be shown) it works great. Consider this example:
# create a dict:
a = {'Generated in Python': True,
'scalar': 1,
'vector': [1,2,3],
'string': 'hello'}
cp = python_dict_to_fortran(a) # creates the c_char_p
dll = CDLL ('lib/json.so') # load the shared library
# define the interface to one of the Fortran routines:
send_json_to_fortran = dll.send_json_to_fortran
send_json_to_fortran.argtypes = [POINTER(c_char_p)]
send_json_to_fortran.restype = None
send_json_to_fortran(byref(cp)) # send the string to fortran
Here we are creating a dictionary, converting it into a string, loading our Fortran DLL, defining the interface to the Fortran function we want to call (send_json_to_fortran()), and then calling it. On the Fortran side, we first define a routine to convert the input string into a JSON structure (using JSON-Fortran):
subroutine c_ptr_to_json(cp,json)
type(c_ptr),intent(in) :: cp
!! a `c_char_p` from python containing a JSON string.
type(json_file),intent(inout) :: json
!! the JSON data structure
character(len=:),allocatable :: fstr
!! string containing the JSON data
if (c_associated(cp)) then
fstr = c_ptr_to_f_string(cp)
call json%load_from_string(fstr) ! parse JSON
deallocate(fstr)
end if
end subroutine c_ptr_to_json
Note that we are using the following routine to convert the Python c_char_p string into a normal Fortran string:
function c_ptr_to_f_string(cp) result(fstr)
type(c_ptr),intent(in) :: cp !! a `c_char_p` from python
character(len=:),allocatable :: fstr
!! the corresponding fortran string
integer :: ilen !! string length
ilen = strlen(cp) ! C library function
block
! convert the C string to a Fortran string
character(kind=c_char,len=ilen+1),pointer :: s
call c_f_pointer(cp,s)
fstr = s(1:ilen)
nullify(s)
end block
end function c_ptr_to_f_string
The example function we are calling from Python is defined below. It is important to note that Fortran is expecting the cp argument to be passed by reference. That is why we use byref(cp) when we call it from Python.
subroutine send_json_to_fortran(cp) &
bind (c, name='send_json_to_fortran')
type(c_ptr),intent(in) :: cp
!! a `c_char_p` from python containing a JSON string.
type(json_file) :: json
call c_ptr_to_json(cp,json)
! do something with the data
! (in this case, just print it)
call json%print_file()
call json%destroy() ! free memory
end subroutine send_json_to_fortran
This prints the JSON structure (from Fortran) like so:
{
"Generated in Python": true,
"scalar": 1,
"vector": [
1,
2,
3
],
"string": "hello"
}
Note that or send_json_to_fortran() routine has the BIND(c) attribute to make it callable from C (or in this case, Python). Once we have the JSON data on the Fortran side, we can also access the variables, query what's in the structure, get the values, etc., like we can any JSON data. If you want to modify the data on the Fortran side, keep reading.
Sending a JSON structure from Fortran to Python
The inverse operation (sending arbitrary data using JSON from Fortran back to Python) is also possible, but a little bit more complicated. First, in Fortran, we will define a container type that will hold a deferred-length (allocatable) string:
type,public :: container
character(len=:),allocatable :: str
end type container
The reason we need this is so that we can point to it with a pointer, which will be come apparent shortly. We then define a Fortran function to convert a JSON structure into a string, and then allocate a pointer to a container to hold the string:
subroutine json_to_c_ptr(json,cp,destroy)
implicit none
type(json_file),intent(inout) :: json !! JSON data
type(c_ptr) :: cp
!! a pointer to a container
!! containing the JSON data as a string
logical,intent(in) :: destroy
!! to also destroy the JSON file
!! (must be destroyed on the fortran
!! side somewhere to prevent memory leak)
character(len=:),allocatable :: str
type(container),pointer :: c
call json%print_to_string(str)
if (destroy) call json%destroy()
allocate(c)
c%str = str
cp = c_loc(c)
end subroutine json_to_c_ptr
We also will need a routine that is callable from Python to get the length of a string within a container:
function get_string_length(cp) result(ilen) &
bind (c, name='get_string_length')
type(c_ptr),intent(in) :: cp !! pointer to a container
integer(c_int) :: ilen !! the length of the string
type(container),pointer :: c
ilen = 0
if (c_associated(cp)) then
call c_f_pointer (cp, c)
if (allocated(c%str)) ilen = len(c%str)
end if
end function get_string_length
Now, we also need another routine callable from Python to populate a c_char_p string buffer:
subroutine populate_char_string(cp,string) &
bind(c,name='populate_char_string')
type(c_ptr),intent(in) :: cp !! pointer to a container
type(c_ptr),intent(inout) :: string
!! a preallocated string buffer that
!! the string will copied into
type(container),pointer :: c
if (c_associated(cp)) then
call c_f_pointer (cp, c)
if (allocated(c%str)) then
call f_string_to_c_ptr(c%str, string)
end if
end if
end subroutine populate_char_string
Finally, we also need a way to destroy the string that we are creating in Fortran. This has to be done in Fortran, since the container pointer was allocated in Fortran. Python will garbage collect its own variables, but it doesn't know how to garbage collect a Fortran pointer, so we'll have to call this routine from Python when we no longer need the variable.
subroutine destroy_string(cp) &
bind (c, name='destroy_string')
type(c_ptr),intent(inout) :: cp
!! pointer to a container
type(container),pointer :: c
if (c_associated(cp)) then
call c_f_pointer (cp, c)
if (allocated(c%str)) deallocate(c%str)
deallocate(c)
end if
cp = c_null_ptr
end subroutine destroy_string
Now, say we have the following Fortran routine that generates some JSON data and sends it to Python:
subroutine get_json_from_fortran(cp) &
bind (c, name='test_send_json_to_python')
type(c_ptr) :: cp
!! pointer to a container containing a json string
type(json_file) :: json
! sample data:
call json%add('Generated in Fortran', .true.)
call json%add('scalar', 1)
call json%add('vector', [1,2,3])
call json%add('string', 'hello')
call json%add('string array', ['1','2','3'])
! convert it to a c_ptr (and destroy JSON structure)
call json_to_c_ptr(json,cp,destroy=.true.)
end subroutine get_json_from_fortran
On the Python side, we have:
# declare some more interfaces to the Fortran routines:
get_string_length = dll.get_string_length
get_string_length.argtypes = [POINTER(c_void_p)]
get_string_length.restype = c_int
populate_char_string = dll.populate_char_string
populate_char_string.argtypes = [POINTER(c_void_p),POINTER(c_char_p)]
populate_char_string.restype = None
destroy_string = dll.destroy_string
destroy_string.argtypes = [POINTER(c_void_p)]
destroy_string.restype = None
def fortran_str_to_python_dict(cp,destroy=True):
# get the length of the string:
length = c_int()
length = get_string_length(byref(cp))
# preallocate a string buffer of the correct size to hold it:
s = c_char_p()
s = cast(create_string_buffer(b' '.ljust(length)),c_char_p)
# convert it to a normal python string:
populate_char_string(byref(cp),byref(s))
string = s.value.decode()
if (destroy):
# destroy it on the Fortran side:
destroy_string(byref(cp))
return json.loads(string)
The fortran_str_to_python_dict() routine is the main routine that converts the string that comes from Fortran into a Python dictionary. First we get the string length, then allocate a string buffer to hold it, we populate the buffer, and then parse the JSON. Note that we also have the option to destroy the Fortran string (which if you recall was allocated as a pointer in json_to_c_ptr() so it needs to be deallocated to prevent a memory leak).
Now we can use these routines from Python like so:
get_json_from_fortran = dll.get_json_from_fortran
get_json_from_fortran.argtypes = [POINTER(c_void_p)]
get_json_from_fortran.restype = None
cp = c_void_p()
get_json_from_fortran(byref(cp))
d = fortran_str_to_python_dict(cp,destroy=True)
print(json.dumps(d, indent=2))
Which prints the result:
{
"Generated in Fortran": true,
"scalar": 1,
"vector": [
1,
2,
3
],
"string": "hello",
"string array": [
"1",
"2",
"3"
]
}
From Python to Fortran and Back
Finally, another use case would be to generate some data in Python, pass it to Fortran where it is modified, and then return it back to Python. In this case, we need to use the container approach, with a few additional routines. On the Python side:
c_ptr_to_container_c_ptr = dll.c_ptr_to_container_c_ptr
c_ptr_to_container_c_ptr.argtypes = [POINTER(c_char_p),POINTER(c_void_p)]
c_ptr_to_container_c_ptr.restype = None
def python_str_to_container(s):
cp = python_str_to_fortran(s)
ccp = c_void_p()
c_ptr_to_container_c_ptr(byref(cp),byref(ccp))
return ccp
def python_dict_to_container(d):
return python_str_to_container(json.dumps(d))
And we also need a routine on the Fortran side to convert from the string buffer format to the container format:
subroutine c_ptr_to_container_c_ptr(cp,ccp) &
bind (c, name='c_ptr_to_container_c_ptr')
type(c_ptr),intent(in) :: cp
!! a `c_char_p` from python
type(c_ptr),intent(out) :: ccp
!! pointer to a container that contains the string
character(len=:),allocatable :: str
!! fortran version of the string
type(container),pointer :: c
!! container to hold the string
if (c_associated(cp)) then
! get the fortran string:
str = c_ptr_to_f_string(cp)
! return a pointer to the container:
allocate(c)
c%str = str
ccp = c_loc(c)
else
ccp = c_null_ptr
end if
end subroutine c_ptr_to_container_c_ptr
Here is an example use case:
subroutine modify_json_in_fortran(cp) &
bind (c, name='modify_json_in_fortran')
type(c_ptr),intent(inout) :: cp !! a pointer to a container
type(json_file) :: json
type(container),pointer :: c
if (c_associated(cp)) then
call c_f_pointer (cp, c)
if (allocated(c%str)) then
! parse JSON:
call json%load_from_string(c%str)
!do something with the data:
call json%print_file()
call json%add('Added in Fortran', [9,10])
! convert it to a c_ptr (and destroy JSON structure)
call json_to_c_ptr(json,cp,destroy=.true.)
end if
end if
end subroutine modify_json_in_fortran
Which we can call from Python like so:
modify_json_in_fortran = dll.modify_json_in_fortran
modify_json_in_fortran.argtypes = [POINTER(c_void_p)]
modify_json_in_fortran.restype = None
cp = python_dict_to_container(a)
modify_json_in_fortran(byref(cp))
d = fortran_str_to_python_dict(cp,destroy=True)
print('')
print('modified by Fortran:')
print(json.dumps(d, indent=2))
Which prints:
{
"Generated in Python": true,
"scalar": 1,
"vector": [
1,
2,
3
],
"string": "hello"
}
modified by Fortran:
{
"Generated in Python": true,
"scalar": 1,
"vector": [
1,
2,
3
],
"string": "hello",
"Added in Fortran": [
9,
10
]
}
Summary
Once all the interface routines are created, this technique is fairly easy to use and seems to work great. It is probably fine for a lot of use cases where extreme performance is not required. For example, it may be a good enough technique for exchanging data from a Python GUI and a core application written in Fortran. However, for very large data sets, this may not be very efficient since it does involve converting everything to and from strings.
See also
Updated on 7 June 2018: the code is now available on GitHub.
Altair's ascent stage docking with Orion [NASA]
The lunar gravity field is bumpy. Thus a spacecraft in a low lunar orbit for an extended period of time must perform periodic correction maneuvers or the orbit eventually becomes unstable and can crash into the Moon. When the Orion spacecraft was being designed during the late Constellation Program, the requirements of global surface access sortie missions, six-month polar outpost missions, anytime aborts, and autonomous operation of the vehicle while the crew was on the surface, required lots of analysis to figure out how to make the architecture work. And then the program was canceled.
The reference low-lunar orbit maintenance algorithm we used during Constellation was fairly simple. Orion initially inserted into a retrograde circular parking orbit with a 100 km altitude. During orbit propagation, if the altitude was reduced by a specified deadband altitude (say 10 km) the vehicle would perform a periapsis raise maneuver at the next apoapsis to bring the periapsis back up to 100 km. The apoapsis radius was not controlled, and indeed all the other orbital elements were also allowed to float. This would continue until the crew was ready to return (on the Altiar lander ascent stage), at which point Orion would circularize the orbit and perform a plane change to allow Altair to perform an in-plane ascent.
The \(\Delta v\) maneuver to perform at apoapsis of an orbit with a periapsis radius \(r_{pi}\) and apoapsis radius \(r_{ai}\), to target a periapsis radius \(r_{pf}\) can be computed by:
- Initial apoapsis velocity: \(v_{ai} = \sqrt{ 2 \mu \frac{r_{pi}}{r_{ai} (r_{ai}+r_{pi})} }\)
- Final apoapsis velocity: \(v_{af} = \sqrt{ 2 \mu \frac{r_{pf}}{r_{ai} (r_{ai}+r_{pf})} }\)
- Maneuver magnitude: \(\Delta v = v_{af} - v_{ai}\)
We can implement this method using the Fortran Astrodynamics Toolkit fairly easily. We also need an integrator with an event finding capability (we can use DDEABM, a variable step size variable order Adams method). We propagate from the initial insertion state until the deadband altitude is violated, then propagate to the next apoapsis where the periapsis raise maneuver is performed, then repeat as necessary. The code implementing this scheme can be found on GitHub here. For the force model, we are using a 20 x 20 GRAIL geopotential model of the Moon (not available during Constellation, when we used the Lunar Prospector derived LP150Q model). A contour plot of the results for 100 km retrograde orbits with a 10 km altitude deadband and 10 days of propagation time are shown below. Each maneuver ends up being about 2 m/s for this example, and the bumpiness of the gravity field is evident. Some insertion orbits require no maneuvers, while others require about 18 m/s over the 10 day period.
See also
- J. Williams, S. M. Stewart, D. E. Lee, E. C. Davis, G. E. Condon, J. S. Senent. "The Mission Assessment Post Processor (MAPP): A New Tool for Performance Evaluation of Human Lunar Missions", 20th AAS/AIAA Space Flight Mechanics Meeting, 14-17 February 2010, San Diego, CA.
- F.G. Lemoine, S. Goossens, T.J. Sabaka, et al. GRGM900C: A degree 900 lunar gravity model from GRAIL primary and extended mission data. Geophysical Research Letters. 2014;41(10):3382-3389.
- R. B. Roncoli, Lunar Constants and Models Document, JPL D-32296, Sept 23, 2005.
- M. Beckman, L. Rivers, Stationkeeping for the Lunar Reconnaissance Orbiter (LRO), 20th International Symposium on Space Flight Dynamics, 24-28 Sep. 2007.
A standard way to access a program's command line arguments wasn't added to Fortran until 2003. Command line arguments were one of those things (like environment variables, the file system, and the operating system shell) that the Fortran standards committee pretended didn't exist or wasn't important. Of course, that was and is ridiculous. So for decades, all Fortran compilers provided a non-standard way to get the command line arguments (for example, GETARG in gfortran). Currently, the standards committee still pretends that the file system doesn't exist (for example, there's no standard way to change an application's working directory, or list the files in a directory, or even to check for the existence of directories). But I digress.
The two routines you need to know are COMMAND_ARGUMENT_COUNT and GET_COMMAND_ARGUMENT. Frequently in examples on the internet (even the gfortran example) you will see people just declare a really long string and use GET_COMMAND_ARGUMENT to get the command line argument and hope for the best. That's terrible of course, since this routine does have the ability to return the actual length of the argument, so you should use that to allocate a string with the correct size to hold it.
Consider the following example:
module command_line_module
implicit none
private
type,public :: string
!! a variable-length string type
character(len=:),allocatable :: str
end type
public :: get_command_arguments
contains
subroutine get_command_arguments(nargs, args)
!! return all the command line arguments.
!! (each one is the correct size)
implicit none
integer,intent(out) :: nargs !! number of arguments
type(string),dimension(:),allocatable,intent(out) :: args
integer :: i, ilen
nargs = command_argument_count()
allocate(args(0:nargs))
do i = 0, nargs
call get_command_argument(i,length=ilen)
allocate(character(len=ilen) :: args(i)%str)
call get_command_argument(i,value=args(i)%str)
end do
end subroutine get_command_arguments
end module command_line_module
Here we have to define our own string type (string), so we can create an array of deferred-length strings (since a useful varying-length string type is another thing Fortran pretends is not important). See StringiFor for a full-featured variable-string type, but for our purposes here the simple one is fine. The procedure in the get_command_arguments() routine is to check for the size of each argument, allocate a string large enough to hold it, and then retrieve the value. It returns an array of string variables containing the command line arguments. With a wrapper module like this, getting command line arguments is effortless, and there is no reason anymore to use the old non-standard routines.
It looks like the Fortran community has lost another venerable NASA Fortran library to C++. I recently noticed that the latest release of NASA's Global Reference Atmosphere Model (GRAM) is now C++. From the release link:
Earth-GRAM 2016 is now available as an open-source C++ computer code that can run on a variety of platforms including PCs and UNIX stations. The software provides a model that offers values for atmospheric parameters such as density, temperature, winds, and constituents for any month and at any altitude and location within the Earth's atmosphere. Earth-GRAM 2010 is available in FORTRAN. Similar computer models of Mars, Venus, Titan, and Neptune are available for diverse mission applications through NASA’s Space Environments and Effects (SEE) Program, a partnership with industry, academia, and other government agencies that seeks to develop more reliable, more effective spacecraft.
The Earth-GRAM model, originally written in Fortran 77, is produced by MSFC and has been around for a while, with many revisions (i.e., GRAM-86, GRAM-88, GRAM-90, GRAM-95, GRAM-99, GRAM 2007, and GRAM 2010). In the 2010 revision, they converted the old Fortran 77 code to Fortran 90, which I thought was a step in the right direction (although maybe 20 years too late). It seems like they could have used Fortran 2003 and provided a C interface using the standard Fortran/C Interoperability feature, which would have enabled calling it from a variety of languages (including C++).
I remain utterly unconvinced that C++ is the programming language that engineers and scientists should be using. Are CS people taking over programming jobs once done by engineers? Is this a failure of university engineering departments? Although it seems like even CS people don't want to be using C++ anymore. Every time I turn around a new programming language (Java, Rust, Go, D, C#, Swift, ...) is created by CS people fed up with the C++ dumpster fire. So I don't know who is telling engineers that C++ is the future. And the perception that Fortran is an obsolete programming language is very strong, even within the engineering/technical disciplines that it was designed for. Sometimes this is due to sheer ignorance, but sometimes it is due to ancient Fortran 77 legacy code that people have to deal with. I sympathize with them, Fortran 77 code is terrible and needs to go away, but it should be modernized, not thrown out and replaced with C++. Think of the children! (and the memory leaks!)
Parting words of wisdom
Security tips when programming in C (2017 edition):
1) Stop typing
2) Delete what you've already typed
— ryan huber (@ryanhuber) June 21, 2017
See also
- Earth-GRAM 2016 is now available!
- Earth Global Reference Atmospheric Model (Earth-GRAM 2010) -- The last Fortran 90 release. [request here]
- F.W. Leslie, C.G. Justus, "The NASA Marshall Space Flight Center Earth Global Reference Atmospheric Model—2010 Version", NASA/TM—2011–216467, 2011
- ModelWeb Catalogue and Archive -- other atmosphere models listed on GSFC ModelWeb
- D. L. Johnson, R. Blocker, C. G. Justus, Global Reference Atmosphere Model (GRAM), NASA Tech Brief, MFS-28397, Jan 01, 1993
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
The IAU Standards of Fundamental Astronomy (SOFA) library implements standard models used in fundamental astronomy. Version 14 ( 2018-01-30) has just been released. According to the release notes, this update includes the following:
- Change in the copyright status of the iau_DAT routine. This is to provide for a user-supplied mechanism for updating the number of leap seconds. SOFA doesn't provide any such mechanism, but they are now allowing you to replace or modify this routine to provide your own without having to rename the routine (all other SOFA routines cannot be modified without renaming them).
- Implementation of two new categories of routines:
- Three new routines for the horizon/equatorial plane coordinates.
- iau_AE2HD — (azimuth, altitude) to (hour angle, declination)
- iau_HD2AE — (hour angle, declination) to (azimuth, altitude)
- iau_HD2PA — parallactic angle
- Six new routines dealing with gnomonic (tangent plane) projections.
- iau_TPORS — solve for tangent point, spherical
- iau_TPORV — solve for tangent point, vector
- iau_TPSTS — project tangent plane to celestial, spherical
- iau_TPSTV — project tangent plane to celestial, vector
- iau_TPXES — project celestial to tangent plane, spherical
- iau_TPXEV — project celestial to tangent plane, vector
- The Astrometry Tools Cookbook, the test program and other supporting files have also been updated.
- Other minor documentation/typographical corrections to various files.
IAU SOFA has a Fortran 77 and a C version. While this is a great resource, I do wonder why the IAU insists on using a source format for the Fortran version that was rendered obsolete 30 years ago? There's no reason for this. This is new code that is regularly kept up to date, so why not modernize it for those of us living in the 21st century?
Consider this SOFA routine:
SUBROUTINE iau_TR ( R, RT )
IMPLICIT NONE
DOUBLE PRECISION R(3,3), RT(3,3)
DOUBLE PRECISION WM(3,3)
INTEGER I, J
DO 2 I=1,3
DO 1 J=1,3
WM(I,J) = R(J,I)
1 CONTINUE
2 CONTINUE
CALL iau_CR ( WM, RT )
END SUBROUTINE iau_TR
All this is doing is transposing a 3x3 matrix. Well, since the late 1990s, to do that in Fortran all you need is:
The original routine is just awful: all caps, line numbers, CONTINUE statements, and 100% unnecessary. Unfortunately, it's code like this that gives Fortran a bad name (reinforcing the perception that Fortran is an arcane and obsolete programming language). Anyone who is compiling this code is probably using a compiler that is at least Fortran 95 compatible. Who on earth is using straight up Fortran 77 compilers at this point? So why restrict such a useful library to a programming style that is so totally out of date? I think it is even probably affecting the efficiency of the code, since the example above required using a temporary array and a function call to do what can now be done with a (presumably much more efficient) intrinsic routine. It wouldn't even be that hard to convert this code to modern Fortran style. It's just low-level math routines, it doesn't have to be object-oriented or anything fancy like that. A while back, I wrote a little Python script to merge all the SOFA files into a single Fortran module (another innovation from the 1990s) so they can take advantage of the automatic interface checking that modules provide. There are also any number of tools out that could be used to automate the fixed to free-form conversion.
See also
