#
# $Id: airfoil.dem,v 1.10 2006/06/30 02:17:21 sfeam Exp $
#
# This demo shows how to use bezier splines to define NACA four
# series airfoils and complex variables to define Joukowski
# Airfoils.  It will be expanded after overplotting in implemented
# to plot Coefficient of Pressure as well.
#		Alex Woo, Dec. 1992
#
# The definitions below follows: "Bezier presentation of airfoils",
# by Wolfgang Boehm, Computer Aided Geometric Design 4 (1987) pp 17-22.
#
#				Gershon Elber, Nov. 1992
#
# m = percent camber
# p = percent chord with maximum camber
print  "NACA four series airfoils by bezier splines"
print  "Will add pressure distribution later with Overplotting"
mm = 0.6
# NACA6xxx
thick = 0.09  
# nine percent  NACAxx09
pp = 0.4
# NACAx4xx
# Combined this implies NACA6409 airfoil
#
# Airfoil thickness function.
#
set xlabel "NACA6409 -- 9% thick, 40% max camber, 6% camber"
x0 = 0.0
y0 = 0.0
x1 = 0.0
y1 = 0.18556
x2 = 0.03571
y2 = 0.34863
x3 = 0.10714
y3 = 0.48919
x4 = 0.21429 
y4 = 0.58214
x5 = 0.35714
y5 = 0.55724
x6 = 0.53571
y6 = 0.44992
x7 = 0.75000
y7 = 0.30281
x8 = 1.00000
y8 = 0.01050
#
# Directly defining the order 8 Bezier basis function for a faster evaluation.
#
bez_d4_i0(x) =     (1 - x)**4
bez_d4_i1(x) = 4 * (1 - x)**3 * x
bez_d4_i2(x) = 6 * (1 - x)**2 * x**2
bez_d4_i3(x) = 4 * (1 - x)**1 * x**3
bez_d4_i4(x) =                  x**4

bez_d8_i0(x) =      (1 - x)**8
bez_d8_i1(x) =  8 * (1 - x)**7 * x
bez_d8_i2(x) = 28 * (1 - x)**6 * x**2
bez_d8_i3(x) = 56 * (1 - x)**5 * x**3
bez_d8_i4(x) = 70 * (1 - x)**4 * x**4
bez_d8_i5(x) = 56 * (1 - x)**3 * x**5
bez_d8_i6(x) = 28 * (1 - x)**2 * x**6
bez_d8_i7(x) =  8 * (1 - x)    * x**7
bez_d8_i8(x) =                   x**8


m0 = 0.0
m1 = 0.1
m2 = 0.1
m3 = 0.1
m4 = 0.0
mean_y(t) = m0 * mm * bez_d4_i0(t) + \
	    m1 * mm * bez_d4_i1(t) + \
	    m2 * mm * bez_d4_i2(t) + \
	    m3 * mm * bez_d4_i3(t) + \
	    m4 * mm * bez_d4_i4(t)

p0 = 0.0
p1 = pp / 2
p2 = pp
p3 = (pp + 1) / 2
p4 = 1.0
mean_x(t) = p0 * bez_d4_i0(t) + \
	    p1 * bez_d4_i1(t) + \
	    p2 * bez_d4_i2(t) + \
	    p3 * bez_d4_i3(t) + \
	    p4 * bez_d4_i4(t)

z_x(x) = x0 * bez_d8_i0(x) + x1 * bez_d8_i1(x) + x2 * bez_d8_i2(x) + \
	 x3 * bez_d8_i3(x) + x4 * bez_d8_i4(x) + x5 * bez_d8_i5(x) + \
	 x6 * bez_d8_i6(x) + x7 * bez_d8_i7(x) + x8 * bez_d8_i8(x)

z_y(x, tk) = \
   y0 * tk * bez_d8_i0(x) + y1 * tk * bez_d8_i1(x) + y2 * tk * bez_d8_i2(x) + \
   y3 * tk * bez_d8_i3(x) + y4 * tk * bez_d8_i4(x) + y5 * tk * bez_d8_i5(x) + \
   y6 * tk * bez_d8_i6(x) + y7 * tk * bez_d8_i7(x) + y8 * tk * bez_d8_i8(x)

#
# Given t value between zero and one, the airfoild curve is defined as
# 
#			c(t) = mean(t1(t)) +/- z(t2(t)) n(t1(t)),
#
# where n is the unit normal to the mean line. See the above paper for more.
#
# Unfortunately, the parametrization of c(t) is not the same for mean(t1)
# and z(t2). The mean line (and its normal) can assume linear function t1 = t,
#                                                     -1
# but the thickness z_y is, in fact, a function of z_x  (t). Since it is
# not obvious how to compute this inverse function analytically, we instead
# replace t in c(t) equation above by z_x(t) to get:
# 
#			c(z_x(t)) = mean(z_x(t)) +/- z(t) n(z_x(t)),
#
# and compute and display this instead. Note we also ignore n(t) and assumes
# n(t) is constant in the y direction,
#

airfoil_y1(t, thick) = mean_y(z_x(t)) + z_y(t, thick)
airfoil_y2(t, thick) = mean_y(z_x(t)) - z_y(t, thick)
airfoil_y(t) = mean_y(z_x(t))
airfoil_x(t) = mean_x(z_x(t))
unset grid
unset zeroaxis
set parametric
set xrange [-0.1:1.1]
set yrange [-0.1:.7]
set trange [ 0.0:1.0]
set title "NACA6409 Airfoil"
plot airfoil_x(t), airfoil_y(t) title "mean line" w l lt 2, \
     airfoil_x(t), airfoil_y1(t, thick) title "upper surface" w l lt 1, \
     airfoil_x(t), airfoil_y2(t, thick) title "lower surface" w l lt 1
pause -1 "Press Return"
mm = 0.0
pp = .5
thick = .12
set title "NACA0012 Airfoil"
set xlabel "12% thick, no camber -- classical test case"
plot airfoil_x(t), airfoil_y(t) title "mean line" w l lt 2, \
     airfoil_x(t), airfoil_y1(t, thick) title "upper surface" w l lt 1, \
     airfoil_x(t), airfoil_y2(t, thick) title "lower surface" w l lt 1
pause -1 "Press Return"
set title ""
set xlab ""
set key box
set parametric
set samples 100
set isosamples 10
set style data lines
set style function lines
print  "Joukowski Airfoil using Complex Variables" 
set title "Joukowski Airfoil using Complex Variables" offset 0,0
set time
set yrange [-.2 : 1.8]
set trange [0: 2*pi]
set xrange [-.6:.6]
zeta(t) = -eps + (a+eps)*exp(t*{0,1})
eta(t) = zeta(t) + a*a/zeta(t)
eps = 0.06
a =.250
set xlabel "eps = 0.06 real"
plot real(eta(t)),imag(eta(t))
pause -1 "Press Return"
eps = 0.06*{1,-1}
set xlabel "eps = 0.06 + i0.06"
plot real(eta(t)),imag(eta(t))
pause -1 "Press Return"
reset