Degenerate Conic

Fortran and Pyplot

Python's matplotlib.pyplot is a very nice collection of functions that provide an easy Matlab-like interface for data plotting. It can be used to generate quite professional looking plots. There is a lot of information on the internet about calling Fortran from Python, but what if all you want to do is generate some plots from your (modern) Fortran program? With this in mind, I've created a very simple Fortran wrapper for Pyplot, allowing you to make pretty good plots by writing only Fortran code. Consider the following example:

program test

use,intrinsic :: iso_fortran_env, only: wp => real64
use pyplot_module

implicit none

real(wp),dimension(100) :: x,sx,cx,tx
type(pyplot) :: plt
integer :: i

!generate some data:
x = [(real(i,wp), i=0,size(x)-1)]/5.0_wp
sx = sin(x)
cx = cos(x)
tx = sx * cx

!plot it:
call plt%initialize(grid=.true.,xlabel='angle (rad)',&
                    title='python_plot test',legend=.true.)
call plt%add_plot(x,sx,label='\$\sin (x)\$',&
                  linestyle='b-o',markersize=5,linewidth=2)
call plt%add_plot(x,cx,label='\$\cos (x)\$',&
                  linestyle='r-o',markersize=5,linewidth=2)
call plt%add_plot(x,tx,label='\$\sin (x) \cos (x)\$',&
                  linestyle='g-o',markersize=2,linewidth=1)
call plt%savefig('test.png')

end program test

The main user interface is the pyplot class, which has methods such as initialize, add_plot, and savefig. This example produces the following plot:

All the module does is generate a Python script, and then runs it and saves the result. The Python call is completely transparent to the Fortran user.

I posted the project to GitHub (pyplot-fortran). Maybe someone else will find it useful, or help to expand it.

See also

1. matplotlib
2. F2PY Users Guide and Reference Manual --Fortran to Python interface generator

json-fortran 4.0.0

I just tagged the 4.0.0 release of json-fortran. This is the first release with Unicode support (thanks to Izaak Beekman). Who says there are no good open source Fortran projects on the internet?

The Unicode build of the library is optional (and only enabled using the preprocessor directive USE_UCS4). Currently, this only works with Gfortran. It doesn't yet work with the Intel Fortran Compiler, which is lagging behind on Unicode support. The Fortran standard supports Unicode via the selected_char_kind function, which can be used to specify the character set used for a character string, like so:

integer,parameter :: u = selected_char_kind('ISO_10646')
character(kind=u,len=11) :: string = u_'Hello World'

Pikaia

I've started a new project on Github for a new modern object-oriented version of the Pikaia genetic algorithm code. This is a refactoring and upgrade of the Pikaia code from the High Altitude Observatory. The original code is public domain and was written by Paul Charbonneau & Barry Knapp (version 1.2 was released in 2003). The new code differs from the old code in the following respects:

• The original fixed-form source (FORTRAN 77) was converted to free-form source.
• The code is now object-oriented Fortran 2003/2008. All user interaction is now through a class, which is the only public entity in the module.
• All real variables are now 8 bytes (a.k.a. "double precision"). The precision can also be changed by the user by modifying the wp parameter.
• The random number generator was replaced with calls to the intrinsic Fortran random_number function.
• There are various new options (e.g., a convergence window with a tolerance can be specified as a stopping condition, and the user can specify a subroutine for reporting the best solution at each generation).
• Mapping the variables to be between 0 and 1 now occurs internally, rather than requiring the user to do it.
• Can now include an initial guess in the initial population.

This is another example of how older Fortran codes can be refactored to bring them up to date with modern computing concepts without doing that much work. In my opinion, the object-oriented Fortran 2003 interface is a lot better than the old FORTRAN 77 version. For one, it's easier to follow. Also, data can now easily be passed into the user fitness function by extending the pikaia class and putting the data there. No common blocks necessary!

References

1. Original Pikaia code and documentation [High Altitude Observatory]

Complex Step Differentiation

Complex step differentiation is a method for estimating the derivative of a function \(f(x)\) using the equation:

• \(f'(x) \approx \frac{ Im[f(x + ih)] }{h}\)

Unlike finite-difference formulas, the complex-step formula does not suffer from roundoff errors due to subtraction, so a very small value of the step size \(h\) can be used, producing a much more accurate derivative. The example shown below is for the function: \(f(x) = e^x + \sin(x)\). The derivative error is compared to three finite-difference formulas:

• \(f'(x) \approx \frac{f(x+h) - f(x)}{h}\) (Two-point forward difference)
• \(f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}\) (Two-point central difference)
• \(f'(x) \approx \frac{f(x-2h) -8 f(x-h) + 8 f(x+h) - f(x+2h) }{12 h}\) (Four-point central difference)

For large values of \(h\), truncation error dominates. The finite-difference error decreases as \(h\) is decreased, until roundoff error begins to dominate (around \(10^{-8}\) for the forward difference, \(10^{-5}\) for the 2-point central difference, and \(10^{-3}\) for the 4-point central difference). The complex-step estimate, however, can produce a derivative estimate to machine precision for any value of \(h < 10^{-8}\).

See also

1. J.R.R.A. Martins, I.M. Kroo, J.J. Alonso, "An Automated Method for Sensitivity Analysis using Complex Variables", AIAA-2000-0689.
2. complex_step.f90 in the Fortran Astrodynamics Toolkit.

Revenge of the Animated GIF

Yes, now animated GIFs can be created using only Fortran code. I present FGIF, now on GitHub. It is a modified and updated version of the public domain code found on the Fortran Wiki. It was used to create the image shown here (which is an example of a circle illusion).

For more information on the GIF format, see:

• What's In A GIF
• GIF Format

It was interesting to discover that animated GIFs are actually a Netscape extension to the GIF file format. If you open any animated GIF in a hex editor, you will find the string "NETSCAPE2.0".

json-fortran 1.0.0

I just tagged the 1.0.0 release of json-fortran. You're welcome, interwebs.

Fortran Astrodynamics Toolkit

I'm starting a new project on GitHub: the Fortran Astrodynamics Toolkit. Hardly anyone is developing open source orbital mechanics software for modern Fortran, so the time has come. Most of the code from this blog will eventually find its way there. My goal is to include modern Fortran implementations of all the standard orbital mechanics algorithms such as:

• Lambert solvers
• Kepler propagators
• ODE solvers (Runge-Kutta, Nyström, Adams)
• Orbital element conversions

Then we'll see where it goes from there. The code will be released under a permissive BSD style license.

JSON + Fortran

Introducing json-fortran, an easy-to-use JSON API written in modern Fortran. As far as I know, it's currently the only publicly-available production-ready JSON API for Fortran that does not involve an interface to C libraries. The code is hosted at GitHub, and is released under a BSD-style license.

Fortran users may find JSON quite useful as a configuration file format. Typically, large and complicated Fortran codes (say, in the fields of climate modeling or trajectory optimization) require data to be read from configuration files. JSON has many advantages over a roll-your-own file format, or a Fortran namelist, namely:

• It's a standard.
• It's human-readable and human-editable.
• API's exist for many other programming languages.

Consider the example of a program to propagate a spacecraft state vector. The required information is the initial time t0, step size dt, final time tf, central body gravitational parameter mu, and the 6-element state vector x0. The JSON configuration file might look something like this:

{  
  "t0": 0.0,  
  "dt": 1.0,  
  "tf": 86400.0,  
  "mu": 398600.4418,  
  "x0": [ 10000.0,  
         10000.0,  
         10000.0,  
        1.0,  
         2.0,  
         3.0  
  ]  
}  

The code would look something like this:

program propagate

use json_module

implicit none

real(wp) :: t0, dt, tf, mu  
real(wp),dimension(:),allocatable :: x0  
type(json_file) :: config  
logical :: found

!load the file:  
call config%load_file('config.json')

!read in the data:  
call config%get('t0',t0,found)  
call config%get('dt',dt,found)  
call config%get('tf',tf,found)  
call config%get('mu',mu,found)  
call config%get('x0',x0,found)

!propagate:  
! ...

end program propagate  

Of course, real production code would have more graceful error checking. For example, to make sure the file was properly parsed and that each variable was really present. These features are included in the API, including helpful error messages if something goes wrong.