airfoil.py

Dateiinhalt

#!/usr/bin/env python2
# -*- coding: utf-8 -*-

# Philipp Maierhöfer, 2017

import sys
import os
import argparse

allowed_fileexts = ['.pdf', '.ps', '.eps', '.png', '.svg']

parser = argparse.ArgumentParser(description=
    """Use the Joukowski function to plot an airfoil as the image of a circle
    around z0=-a+i*b and corresponding streamlines. If no output file is
    specified, print the plot on the screen. This only works when the matplotlib
    backend is configured properly (set 'backend' in your matplotlibrc file to
    an available output driver). Requires matplotlib and numpy.""")
parser.add_argument('-o', '--output', dest='outfile', metavar='outfile',
                    default=None,
                    help='output file name, supported file formats are {}'
                         .format(' '.join(allowed_fileexts)))
parser.add_argument('-v', dest='vstep', metavar='vstep', type=float,
                    default=0.1,
                    help='the distance of the streamlines')
parser.add_argument('a', type=float, help='the negative real part of z0')
parser.add_argument('b', type=float, help='the imaginary part of z0')
args = parser.parse_args()

z0 = -args.a + 1j*args.b

if args.outfile:
    fileext = os.path.splitext(args.outfile)[1].lower()
    if fileext not in allowed_fileexts:
        print 'Error: Unsupported file extension \'' + fileext + '\'.'
        print '       Supported formats are', allowed_fileexts
        sys.exit(1)

try:
    import matplotlib
except ImportError:
    print "Error: matplotlib not found."
    sys.exit(1)
import numpy
import matplotlib.pyplot as plt

# List backends
#print matplotlib.rcsetup.all_backends
# Print the current backend
#print matplotlib.get_backend()
# Print the currently used matplotlibrc file path
#print matplotlib.matplotlib_fname()

def hinvpart(z, z0):
    return (z0 + numpy.absolute(1-z0) *
            (z + 1j*numpy.sqrt(1-z**2)*numpy.sign(numpy.imag(z))))

def hinv(x, y, z0):
    return (hinvpart(x+1j*y, z0) + 1/hinvpart(x+1j*y, z0))/2

# a small positive regulator
ep = 0.00001
# the number of points to sample in each curve
npoints = 1000

# prepare the data for the airfoil as a polygon;
# need to reverse one of the data sets
t = numpy.linspace(-1,1,npoints)
airfoil_xdata = numpy.concatenate(
    [numpy.real(hinv(t,ep,z0)), numpy.real(hinv(t,-ep,z0))[::-1]])
airfoil_ydata = numpy.concatenate(
    [numpy.imag(hinv(t,ep,z0)), numpy.imag(hinv(t,-ep,z0))[::-1]])

plotrangex = [min(1.3*min(airfoil_xdata), -1.5),
              max(1.3*max(airfoil_xdata), 1.5)]
plotrangey = [min(1.3*min(airfoil_ydata), -1),
              max(1.3*max(airfoil_ydata), 1)]
plt.figure()
plt.title(u'z₀ = {}{:+}i'.format(-args.a, args.b))
plt.xlim(*plotrangex)
plt.ylim(*plotrangey)
plt.axes().set_aspect('equal')

# draw the streamlines; the range in x direction will be larger than plotrangex
u = numpy.linspace(plotrangex[0], plotrangex[1], npoints)
for vn in range(int((plotrangey[1]-plotrangey[0])/args.vstep)):
    v = plotrangey[0] + vn*args.vstep + ep
    plt.plot(numpy.real(hinv(u,v,z0)),
             numpy.imag(hinv(u,v,z0)),
             linestyle='dotted',
             color='blue')

# draw the airfoil as a filled polygon
plt.fill(airfoil_xdata, airfoil_ydata, color='red')

if args.outfile:
    plt.savefig(args.outfile, format=fileext[1:])
else:
    plt.show()