Fortran from python – Python Binding for Fortran Code

As an aeronautical engineer, I received fortran as a legacy. From NAL to working for Rolls-Royce, work life revolves around fortran. This  is the language that I have been programming for 12+ years now. That explains this many fortran related posts here and here.

But from couple of years, I am falling in love with Python. Here in this post, I will walk through the steps to use fortran and python together.

There are many times, when we don’t want to use the methods included in Python and need to use the trusted, optimised, time tested fortran library. This post will show how to create a python binding of the fortran code and use it in python.

I will present the complete code first and then explain all the parts of it.

I will use this simple NACA code library.

naca4_sym(c,t,n) evaluates y(x) for a NACA symmetric 4-digit airfoil. C the chord length. T the maximum relative thickness. N the number of sample points.

file name: naca.f

This is the actual fortran code to be used from python.

subroutine naca4_symmetric ( t, c, n, x, y )

!*****************************************************************************80
!
!! NACA4_SYMMETRIC evaluates y(x) for a NACA symmetric 4-digit airfoil.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 21 May 2014
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Eastman Jacobs, Kenneth Ward, Robert Pinkerton,
! "The characteristics of 78 related airfoil sections from tests in
! the variable-density wind tunnel",
! NACA Report 460, 1933.
!
! Parameters:
!
! Input, real ( kind = 8 ) T, the maximum relative thickness.
!
! Input, real ( kind = 8 ) C, the chord length.
!
! Input, integer ( kind = 4 ) N, the number of sample points.
!
! Input, real ( kind = 8 ) X(N), points along the chord length.
! 0.0 <= X(*) <= C.
!
! Output, real ( kind = 8 ) Y(N), for each value of X, the corresponding
! value of Y so that (X,Y) is on the upper wing surface, and (X,-Y) is on the
! lower wing surface.
!
implicit none

integer ( kind = 4 ) n

real ( kind = 8 ) c
real ( kind = 8 ) t
real ( kind = 8 ) x(n)
real ( kind = 8 ) y(n)

y(1:n) = 5.0D+00 * t * c * ( &
0.2969D+00 * sqrt ( x(1:n) / c ) &
+ (((( &
– 0.1015D+00 ) * ( x(1:n) / c ) &
+ 0.2843D+00 ) * ( x(1:n) / c ) &
– 0.3516D+00 ) * ( x(1:n) / c ) &
– 0.1260D+00 ) * ( x(1:n) / c ) )

return
end

file name: aefoilnaca.py
This is the python wrapper using the nacalib created from fortran code

def naca4_sym(c,t,n):
	"""
	NACA4_SYMMETRIC evaluates y(x) for a NACA symmetric 4-digit airfoil.
	C the chord length.
	T the maximum relative thickness.
	N the number of sample points.

	python wrapper for the code by  John Burkardt

	by Sukhbinder Singh

	"""
	import numpy as np
	import nacalib

	x = np.linspace(0.0,c,n)
	y=np.zeros(n)
	y=nacalib.naca4_symmetric(t,c,x,n)
	xy=np.zeros((2,2*n))
	xy[0,0:n]=x
	xy[0,n:]=x[n::-1]
	xy[1,0:n]=-y
	xy[1,n:]=y[n::-1]
	return xy

file name: test.py
This file tests the python wrapper module.

"""
Testing aefoilnaca

"""

import aefoilnaca
import matplotlib.pyplot as plt

c = 10.0 # Chord length
t=0.15   # the maximum relative thickness
n=80     # number of point
m=0.02   # the maximum camber
p=0.4    # the location of maximum camber

y = aefoilnaca.naca4_sym(c,t,n)

plt.plot(y[0,:],y[1,:])
plt.show()

The steps:
First, we create a signature file from naca.f by running

f2py naca.f -m nacalib -h naca.pyf

The signature file is saved to naca.pyf

Second step edit naca.pyf file.

Next we’ll open the naca.pyf and tell F2PY that the argument t,c,x,n are input arguments (use intent(in) attribute) and that the result, i.e. the contents of y after calling Fortran function naca4_symmetric should be returned to Python (use intent(out) attribute). In addition, an array a should be created dynamically using the size given by the input argument n (use depend(n) attribute to indicate dependence relation).

This is the content of the modified naca.pyf

file name: naca.pyf

!    -*- f90 -*-
! Note: the context of this file is case sensitive.

python module nacalib
    interface
        subroutine naca4_symmetric(t,c,n,x,y)
            real(kind=8),intent(in) :: t
            real(kind=8),intent(in) :: c
            integer(kind=4), optional,check(len(x)>=n),depend(x) :: n=len(x)
            real(kind=8) dimension(n),intent(in) :: x
            real(kind=8) dimension(n),depend(n),intent(out) :: y
        end subroutine naca4_symmetric
    end interface
end python module nacalib

! This file was auto-generated with f2py (version:2).
! See http://cens.ioc.ee/projects/f2py2e/

And finally, we build the extension module by running

f2py -c naca.pyf naca.f

This creates the naca.so file and it can be used as shown in the aefoilnaca.py

	import nacalib

	x = np.linspace(0.0,c,n)
	y=np.zeros(n)
	y=nacalib.naca4_symmetric(t,c,x,n)

Hope this is clear and simple, if you have any questions, feel free to comment.

All this code is available at my github account

For further reading refer to this two excellent reference
1. Python as a glue
2. Using f2py – The smart way