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.
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
The upcoming Fortran standard formerly known as Fortran 2015 has a new name: Fortran 2018. It was decided to change it in order to match the expected year of publication. This makes sense. The previous standard (Fortran 2008) was published in 2010.
Waiting for an updated Fortran standard is an exercise in Zen-like patience. Almost a decade after Fortran 2008, we'll get a fairly minor update to the core language. And it will be years after that before it's fully supported by any compiler that most users will have available (gfortran still doesn't have a bug-free implementation of all of Fortran 2008 or even 2003). Fortran was essentially reborn in the Fortran 2003 standard, which was an amazing update that brought Fortran into the modern world. It's a terrific programming language for scientific and technical computing. However, the limitations are all too clear and have been for a long time:
- We need better facilities for generic programming. It's impossible to do some things without having to duplicate code or use "tricks" like preprocessing or include files (see JSON-Fortran for examples).
- We need some kind of exception handling. Fortran does have a floating-point exception handling feature, but honestly, it's somewhat half-baked.
- We need a better implementation of strings. Allocatable strings (introduced in Fortran 2003) are great, but not enough, since they can't be used in all instances where strings are needed.
- We need the language to be generally less verbose. I'm tired to having to type multiple nested SELECT TYPE statements to do something that is a one liner in Python (sure I know Fortran will never be as succinct as Python, but some of the verbosity required for object-oriented Fortran is just perverse).
- We need any number of new features to make it easier to extend the language with third-party libraries (so we don't have to wait two decades for a feature we want).
- We also need the language development process to embrace a more open collaborative model using modern tools (Usenet is not the future). I guess the recent survey was unprecedented, but it's not enough.
Fortran is a programming language that needs a better PR department. Legacy Fortran codes are being rewritten in C++, Python, or even Julia. NASA frequently throws massive Fortran 77 libraries with decades of heritage (e.g., DPTRAJ/ODP, GTDS, SPICELIB) into the trash in order to rewrite it all from the ground up in C++, without ever considering Fortran 2003+ (or maybe not realizing it exists?). The information about modern Fortran on the internet is spotty at best, outdated, or downright wrong (what is the deal with REAL*8?). In popular consciousness Fortran is mostly a punchline (usually something to do with punchcards and your granddad). A language like Python (which was never designed for technical computing) is now seen by many as a superior solution for technical computing. Matlab refers to itself as "the only top programming language dedicated to mathematical and technical computing"! The Julia website lists somewhat misleading benchmarks than implies that C, Julia, and even Lua are faster than Fortran.
Now, get off my lawn, you kids!
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
