Time Conversions with SPICE

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

Posted in Programming Tagged with: , ,

Another One Bites the Dust

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

Posted in Programming Tagged with: , , ,

The New Features of Fortran 2015

The glacially slow pace of Fortran language development continues! The next standard, Fortran 2015, mainly consists of updates for Fortran/C interoperability and new coarray features such as teams. In addition, there are a bunch of minor changes and discrepancy fixes. A few of the new features are:

  • The venerable implicit none statement has been updated to allow for some additional use related to external procedures.
  • The stop code in error stop can now be any integer or character expression.
  • An out_of_range intrinsic was added to allow for testing whether a real or
    integer value can be safely converted to a different real or integer type and kind.
  • You can now declare the kind of the loop variable inside an implied do loop. For example: iarray = [(2*i, integer :: i=1,n)].
  • All procedures are now recursive by default. This is an interesting change. Ever since recursion was added to the language in Fortran 90, a procedure has had to be explicitly declared as recursive. Now, you have to use non_recursive if you don’t want to allow a procedure to be used recursively.
  • Some new syntax to indicate the locality status of variables within a do concurrent loop.
  • There are a lot of new IEEE intrinsic routines for some reason.

Fortran 2015 is expected to be published in 2018.

See also

Posted in Programming Tagged with:

Fortran in the Cloud

You can now try OpenCoarrays and Gfortran in the cloud, courtesy of Zaak Beekman and the Sourcery Institute. Just navigate to http://bit.ly/TryCoarrays and then click “Launch”. This awesome project is enabled by various other awesome tools like Binder, Jupyter, and GitHub. Truly, we are living in the future.

Coarrays are the parallel processing component built into the Fortran language (standardized in Fortran 2008). It uses the Partitioned Global Address Space (PGAS) and Single-Program-Multiple-Data (SPMD) programming models. OpenCoarrays is an open source library to enable coarray usage in Gfortran.

See also

Posted in Programming Tagged with: ,

Xojo and Memory Lane

I happened across this article on a programming language and IDE called Xojo. I never used it under this name, but I have fond memories of using it when it was called RealBASIC. I learned it shortly after it came out in the late 1990s. I found it very intuitive and it was super easy to create fairly nice looking GUI applications for the Mac. In fact, the first optimization code I ever wrote was in RealBASIC, for solving the traveling salesman problem.

In the mid-2000’s when I was a grad student, I used it to build a GUI interface to JPL’s BlackJack GPS receiver. This included reading data from the receiver, as well as sending commands to it. There were numerous floating windows that displayed data and plots (see screenshot at right). At the time, JPL’s interface code only ran on pre-MacOS X systems, which were rapidly becoming obsolete. My RealBASIC app ran great on MacOS X. In fact, I believe JPL requested a copy when they saw it (I have no idea if it lives on somewhere at JPL).

In 2006, I traveled to Germany with a blue clamshell iBook to deliver some updated software to the TerraSAR-X satellite, which included the same type of GPS receiver. I remember during the update process, something unexpected occurred (the details escape me now). In any event, it was required to generate some binary command files to upload to the system in order to set things right. My RealBASIC interface app was already set up to do that, so I fired up my iBook and was able to generate the files in short order. I distinctly remember some surprise from the German engineers that I could do such a thing with this silly looking machine! After uploading the files, things got back on track. TerraSAR-X has been in orbit since June 2007.

I haven’t used it in years, but I’m glad to know that RealBASIC is still around.

See also

Posted in Programming Tagged with:

Flang

It appears that Flang, a new open source Fortran front-end for LLVM, has appeared on GitHub recently with little fanfare. This is apparently the result of NVIDIA’s previously-announced plan to open source the PGI Fortran compiler. Unfortunately, they decided to give it the same name as another earlier attempt to create a Fortran/LLVM compiler (more confusion for poor Fortran programmers). I don’t really know how it compares to Gfortran or Intel (PGI appears to be lagging behind on support for the Fortran 2008 standard). Initial tests by Usenet denizens (yes Usenet still exists) indicate that maybe Flang isn’t quite ready for prime time. Hopefully it will improve with time. I think it’s great news to potentially have another free Fortran compiler available.

See also

Posted in Programming Tagged with: , ,

Fortran Configuration File Formats

String and file manipulation in Fortran isn’t as bad as you’ve heard (OK, it’s bad, but it’s getting better). Sure, modern Fortran only provides the bare minimum of features for this sort of thing out of the box, but various libraries are starting to appear that the modern Fortran programmer can use to deal with a range of different file formats including JSON, INI, CSV, XML, etc. The Fortran user community has been late to the game on this opensource thing. But it does exist and I encourage all Fortraners to use it. You don’t have to roll your own file format or write your own parser. If the existing libraries don’t have what you need, then help to improve them.

The following is a look at a few different options for configuration files in your Fortran program. It is limited to text files. For binary files, I would highly recommend a standard cross-platform format (such as HDF5). For goodness’ sake, do not invent your own binary file format.

Namelists

Fortran namelists are very old-school, but are still part of the standard (indeed this is the only high-level file read/write feature that is built into Fortran). Namelist have a few advantages, but they are generally outweighed by their many disadvantages. Some advantages include:

  • It is super easy to use. The reading and writing of the file is automatically handled by the compiler. A namelist is declared like so:
integer :: a
real :: b
character(len=5) :: c
integer :: iunit,istat

namelist /blah/ a,b,c

!set values:
a = 1
b = 2.0
c = 'three'

and can be written by:

 write(unit=iunit, nml=blah, iostat=istat) 

Which produces a file like this:

&BLAH
 A=          1,
 B=  2.00000000    ,
 C="three",
 /

The file can be read by:

 read(unit=iunit, nml=blah, iostat=istat) 
  • It automatically handles any Fortran variables you throw at it (including multidimensional arrays and derived data types).
  • Variables in the file can be optional. Any variable not present in the file retains the value it had before the namelist was read.

Disadvantages include:

  • Different compilers can and will output slightly different formats, and you have absolutely no control over that.
  • If you want to generate or manipulate your input files using a different language, you’ll probably have to write your own parser to do it. Hardly any other programming language will have a decent namelist parser or writer available. A notable exception is f90nml, a great Python library for reading and writing namelists.
  • The variables in the file must exactly correspond (in name, type, and size) to variables in the code. This turn out to be very restrictive for a lot of reasons. For one thing, it is quite annoying for strings, since it requires you to specify a maximum string length in the variable declaration.
  • It can be difficult to diagnose errors when reading the file (e.g, a syntax error, an unrecognized variable, or an array size or string length larger than the Fortran variable). The compiler will usually just throw an I/O error if anything goes wrong. However, it is possible on some compilers to retrieve the line where the error occurred by rewinding one line and then reading it with a normal READ statement.
  • It is difficult to tell if a variable is actually present in the file (if one isn’t there, then the corresponding variable in the code will just retain whatever value it had before the read was attempted. Thus, the only real way to check if a variable was present or not is to initialize them all to unusual values, and then check for this after the read.)

In my view, there is little reason to start using namelists at this point, since there are now better options available. Unless you have a legacy code that is using them which can’t be changed, I would recommend using something else.

JSON

JSON stands for JavaScript Object Notation, and is a very popular and lightweight data-interchange format. There is an interface for practically all known programming languages, including modern Fortran. My JSON-Fortran library was built on an older Fortran 95 library, which I forked in order to be able to use some of the newer Fortran 2003/2008 features. With JSON-Fortran, the example above would be:

type(json_file) :: json
call json%add('blah.a',1)
call json%add('blah.b',2.0)
call json%add('blah.c','three')

Which produces the JSON file:

{
    "blah": {
        "a": 1, 
        "b": 2.0, 
        "c": "three"
    }
}

Note that f90nml can also be used to transition away from namelists by converting them to JSON using Python. To go from a namelist to JSON, all you have to do is this:

import f90nml
import json
n = f90nml.read('my_namelist.nml')
with open('my_namelist.json', 'w') as outfile:
  json.dump(n, outfile, sort_keys = True, indent = 4)

Using my JSON-Fortran library, you can even convert a JSON file back into a Namelist if that sort of thing is your bag (as shown in a previous post). One disadvantage of JSON is that the standard does not allow for comments in the file. However, some parsers include support for comments (including JSON-Fortran).

INI

INI is a simple, easy-to-write format. Although not really standardized, it is based on the old MS-DOS initialization files. For modern Fortran, the FiNeR library provides a Fortran interface. The INI version of the above file would be:

[blah]
a=1
b=2.0
c=three

Here’s an example of converting an INI file to a JSON file (using FiNeR and JSON-Fortran):

subroutine ini2json(ini_filename, json)

  use json_module
  use finer

  implicit none

  character(len=*),intent(in) :: ini_filename
  type(json_file),intent(out) :: json

  type(file_ini) :: ini
  integer :: ierr
  integer :: i
  character(len=:),dimension(:),allocatable :: secs
  character(len=:),dimension(:),allocatable :: opts

  call ini%load(filename=ini_filename,error=ierr)
  call ini%get_sections_list(secs)
  do i=1,size(secs)
    do
      if (.not. ini%loop(trim(secs(i)),opts)) exit
      call json%add(trim(secs(i))//'.'//&
                    trim(opts(1)),trim(opts(2)))
    end do
  end do
  call ini%free()

end subroutine ini2json

You could call this routine like so:

program main
  use ini2json_module
  use json_module
  implicit none
  type(json_file) :: json
  call ini2json('test.ini',json)
  call json%print_file()
  call json%destroy()
end program main

And it would produce the following result for the file above:

{
  "blah": {
    "a": "1",
    "b": "2.0",
    "c": "three"
  }
}

Note that, in this example, all the variables are returned as strings (getting the numerical values as numbers is left as an exercise to the reader).

Others

[From XKCD]

  • XML (Extensible Markup Language) — I never really liked this one because it is not easily human readable or editable. There are some Fortran XML parsers out there, though.
  • YAML (YAML Ain’t Markup Language) — I’m not as familiar with this one (also I hate recursive acronyms). There is a Fortran interface to it called fortran-yaml, but I haven’t tried it.
  • TOML (Tom’s Obvious, Minimal Language) – I’m not aware of a Fortran parser for this one. This one was created by one of the founders of GitHub.
  • Lua – Another interesting possibility is for your configuration file to be written in a full-up scripting language. There are a couple of Lua-Fortran interfaces out there (e.g., AOTUS).

See also

Posted in Programming Tagged with: , , , , , ,

GFortran 7.1

GCC (which includes GFortran) version 7.1 has just been released. Various features from Fortran 2003-2015 have been added, including:

  • User-defined derived-type input/output (UDTIO).
  • Partial support for derived type coarrays with allocatable and pointer components.
  • Non-constant stop codes and error stop codes.
  • Derived types with allocatable components of recursive type.
  • Intrinsic assignment to polymorphic variables.

Partial support for OpenMP 4.5 is also included. It looks like they also added support for several non-standard extensions that are present in legacy code. GFortran still hasn’t quite made it to full Fortran 2003 support, over a decade after that standard was published. But they are getting there. Thanks to the small number of volunteers who work on this. You are the real heroes!

See also

Posted in Programming Tagged with: ,

JSON-Fortran 5.3.0

JSON-Fortran 5.3.0 is out. This release includes a few minor bug fixes and a couple of new features. For one thing, it is now possible to easily build or modify a JSON structure without using pointer variables. Here’s an example:

program test

use json_module, rk=>json_rk

implicit none

type(json_file) :: f

call f%initialize()

call f%add('name',         'test')
call f%add('inputs.t',     0.0_rk)
call f%add('inputs.x',     [1.0_rk,2.0_rk,3.0_rk])
call f%add('outputs.t',    10.0_rk)
call f%add('outputs.x(1)', 11.0_rk)
call f%add('outputs.x(2)', 22.0_rk)
call f%add('outputs.x(3)', 33.0_rk)

call f%print_file()

call f%destroy()

end program test

As shown, the structure can be built by specifying the variable paths. This way, objects, vectors, vector elements, etc. can all be added. The result of this example is:

{
  "name": "test",
  "inputs": {
    "t": 0.0E+0,
    "x": [
      0.1E+1,
      0.2E+1,
      0.3E+1
    ]
  },
  "outputs": {
    "t": 0.1E+2,
    "x": [
      0.11E+2,
      0.22E+2,
      0.33E+3
    ]
  }
}

So now, in addition to the previous methods of building a JSON structure (reading it from a file, reading it from a string, or building it yourself using pointer variables), there is now another way. Each method has its uses for different applications. Building it from the path strings can be very convenient and also opens up the possibility of easy conversion between a Fortran namelist and a JSON file. I’ll show an example of this in a future post. In a previous post, I showed how to do the reverse (JSON to namelist).

JSON-Fortran is my most popular GitHub project (70 stars as of today). It usually makes it to GitHub’s monthly “trending” Fortran projects page. I’m always interested to hear about what people are using it for. I added a page on the wiki to list some of the projects that use it, so if you’re using it, feel free to add your project to the list.

Posted in Programming Tagged with: ,

Binary Search

Binary search is a search algorithm that finds the position of a target value within a sorted array. The following binary search routine is taken from NASTRAN, a venerable finite element analysis program from NASA:

      SUBROUTINE BISLOC (*,ID,ARR,LEN,KN,JLOC)
C-----
C     BINARY SEARCH - LOCATE KEY WORD 'ID' IN ARRAY 'ARR', 1ST ENTRY
C     IF FOUND, 'JLOC' IS THE MATCHED POSITION IN 'ARR'
C     IF NOT FOUND, NON-STANDARD RETURN
C                                                     I.E.
C     ID  = KEY WORD TO MATCH IN ARR.      MATCH AGAINST 1ST COL OF ARR
C     ARR = ARRAY TO SEARCH.                          ARR(ROW,COL)
C     LEN = LENGTH OF EACH ENTRY IN ARRAY.            LEN=ROW
C     KN  = NUMBER OF ENTRIES IN THE ARR.             KN =COL
C     JLOC= POINTER RETURNED - FIRST WORD OF ENTRY.   MATCHED ROW
C-----
C
      INTEGER  ARR(1)
      DATA     ISWTCH / 16 /
C
      JJ = LEN - 1
      IF (KN .LT. ISWTCH) GO TO 120
      KLO = 1
      KHI = KN
   10 K   = (KLO+KHI+1)/2
   20 J   = K*LEN - JJ
      IF (ID-ARR(J)) 30,90,40
   30 KHI = K
      GO TO 50
   40 KLO = K
   50 IF (KHI-KLO -1) 100,60,10
   60 IF (K .EQ. KLO) GO TO 70
      K   = KLO
      GO TO 80
   70 K   = KHI
   80 KLO = KHI
      GO TO 20
   90 JLOC = J
      RETURN
  100 JLOC = KHI*LEN - JJ
      J    = KN *LEN - JJ
      IF (ID .GT.ARR(J)) JLOC = JLOC + LEN
  110 RETURN 1
C
C     SEQUENTIAL SEARCH MORE EFFICIENT
C
  120 KHI = KN*LEN - JJ
      DO 130 J = 1,KHI,LEN
      IF (ARR(J)-ID) 130,90,140
  130 CONTINUE
      JLOC = KHI + LEN
      GO TO 110
  140 JLOC = J
      GO TO 110
      END

As you can see, it is a scene of FORTRAN 66 horrors! However, it can be refactored into a nice modern Fortran routine like so:

pure function bisloc(id,arr) result(jloc)

!! binary search of a sorted array.

implicit none

integer,intent(in) :: id
   !! key word to match in `arr`
integer,dimension(:),intent(in) :: arr
   !! array to search (it is
   !! assumed to be sorted)
integer :: jloc
   !! the first matched index in 'arr'
   !! (if not found, 0 is returned)

integer :: j,k,khi,klo,n
integer,parameter :: iswtch = 16

n = size(arr)
jloc = 0

if ( n<iswtch ) then

    ! sequential search more efficient
    
    do j = 1 , n
        if ( arr(j)==id ) then
            jloc = j
            return
        elseif (arr(j)>id) then
            return ! error
        endif
    enddo
    return ! error

else

    klo = 1
    khi = n
    k = (klo+khi+1)/2
    do
        j = k
        if ( id<arr(j) ) then
            khi = k
        elseif ( id==arr(j) ) then
            jloc = j
            return
        else
            klo = k
        endif
        if ( khi-klo<1 ) then
            return ! error
        elseif ( khi-klo==1 ) then
            if ( k==klo ) then
                k = khi
            else
                k = klo
            endif
            klo = khi
        else
            k = (klo+khi+1)/2
        endif
    end do

endif

end function bisloc

I did change a couple of things:

  • The original routine was using an obsolete alternate return (note the * argument) for an error condition. The new one just returns 0 if there is an error.
  • The original routine was passing in the array as a vector and treating it like a matrix (it’s either searching a row or column of this matrix, I admit I had a hard time telling). The new version assumes the input is a rank-1 array. You can pass in whatever matrix slices you want (say a column: a(:,1) or a row: a(1,:)).
  • I also made it a pure function, because that’s how I roll.

Note that NASTRAN was released under the NASA Open Source Agreement version 1.3, a somewhat unusual open source license. It isn’t considered a “free” license by the FSF since it includes a statement requiring changes to be your “original creation”.

Unfortunately, the trend at NASA (and elsewhere) is to declare old Fortran code like this to be totally unmaintainable. So, instead of refactoring them into modern Fortran standards, they are thrown away and rewritten from scratch using a more “modern” language like C++ (whose roots go back to 1972!).

See also

Posted in Programming Tagged with: , , , ,