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 )
!! NACA4_SYMMETRIC evaluates y(x) for a NACA symmetric 4-digit airfoil.
! This code is distributed under the GNU LGPL license.
! 21 May 2014
! John Burkardt
! 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.
! 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.
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 ) )
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()
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