--- /dev/null
+#
+# $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
projects / gnuplot / blobdiff