"SfR Fresh" - the SfR Freeware/Shareware Archive 
Member "gnuplot/demo/airfoil.dem" of archive gp423win32.zip:
As a special service "SfR Fresh" has tried to format the requested source page into HTML format using source code syntax highlighting with prefixed line numbers.
Alternatively you can here view or download the uninterpreted source code file.
That can be also achieved for any archive member file by clicking within an archive contents listing on the first character of the file(path) respectively on the according byte size field.
1 #
2 # $Id: airfoil.dem,v 1.10 2006/06/30 02:17:21 sfeam Exp $
3 #
4 # This demo shows how to use bezier splines to define NACA four
5 # series airfoils and complex variables to define Joukowski
6 # Airfoils. It will be expanded after overplotting in implemented
7 # to plot Coefficient of Pressure as well.
8 # Alex Woo, Dec. 1992
9 #
10 # The definitions below follows: "Bezier presentation of airfoils",
11 # by Wolfgang Boehm, Computer Aided Geometric Design 4 (1987) pp 17-22.
12 #
13 # Gershon Elber, Nov. 1992
14 #
15 # m = percent camber
16 # p = percent chord with maximum camber
17 print "NACA four series airfoils by bezier splines"
18 print "Will add pressure distribution later with Overplotting"
19 mm = 0.6
20 # NACA6xxx
21 thick = 0.09
22 # nine percent NACAxx09
23 pp = 0.4
24 # NACAx4xx
25 # Combined this implies NACA6409 airfoil
26 #
27 # Airfoil thickness function.
28 #
29 set xlabel "NACA6409 -- 9% thick, 40% max camber, 6% camber"
30 x0 = 0.0
31 y0 = 0.0
32 x1 = 0.0
33 y1 = 0.18556
34 x2 = 0.03571
35 y2 = 0.34863
36 x3 = 0.10714
37 y3 = 0.48919
38 x4 = 0.21429
39 y4 = 0.58214
40 x5 = 0.35714
41 y5 = 0.55724
42 x6 = 0.53571
43 y6 = 0.44992
44 x7 = 0.75000
45 y7 = 0.30281
46 x8 = 1.00000
47 y8 = 0.01050
48 #
49 # Directly defining the order 8 Bezier basis function for a faster evaluation.
50 #
51 bez_d4_i0(x) = (1 - x)**4
52 bez_d4_i1(x) = 4 * (1 - x)**3 * x
53 bez_d4_i2(x) = 6 * (1 - x)**2 * x**2
54 bez_d4_i3(x) = 4 * (1 - x)**1 * x**3
55 bez_d4_i4(x) = x**4
56
57 bez_d8_i0(x) = (1 - x)**8
58 bez_d8_i1(x) = 8 * (1 - x)**7 * x
59 bez_d8_i2(x) = 28 * (1 - x)**6 * x**2
60 bez_d8_i3(x) = 56 * (1 - x)**5 * x**3
61 bez_d8_i4(x) = 70 * (1 - x)**4 * x**4
62 bez_d8_i5(x) = 56 * (1 - x)**3 * x**5
63 bez_d8_i6(x) = 28 * (1 - x)**2 * x**6
64 bez_d8_i7(x) = 8 * (1 - x) * x**7
65 bez_d8_i8(x) = x**8
66
67
68 m0 = 0.0
69 m1 = 0.1
70 m2 = 0.1
71 m3 = 0.1
72 m4 = 0.0
73 mean_y(t) = m0 * mm * bez_d4_i0(t) + \
74 m1 * mm * bez_d4_i1(t) + \
75 m2 * mm * bez_d4_i2(t) + \
76 m3 * mm * bez_d4_i3(t) + \
77 m4 * mm * bez_d4_i4(t)
78
79 p0 = 0.0
80 p1 = pp / 2
81 p2 = pp
82 p3 = (pp + 1) / 2
83 p4 = 1.0
84 mean_x(t) = p0 * bez_d4_i0(t) + \
85 p1 * bez_d4_i1(t) + \
86 p2 * bez_d4_i2(t) + \
87 p3 * bez_d4_i3(t) + \
88 p4 * bez_d4_i4(t)
89
90 z_x(x) = x0 * bez_d8_i0(x) + x1 * bez_d8_i1(x) + x2 * bez_d8_i2(x) + \
91 x3 * bez_d8_i3(x) + x4 * bez_d8_i4(x) + x5 * bez_d8_i5(x) + \
92 x6 * bez_d8_i6(x) + x7 * bez_d8_i7(x) + x8 * bez_d8_i8(x)
93
94 z_y(x, tk) = \
95 y0 * tk * bez_d8_i0(x) + y1 * tk * bez_d8_i1(x) + y2 * tk * bez_d8_i2(x) + \
96 y3 * tk * bez_d8_i3(x) + y4 * tk * bez_d8_i4(x) + y5 * tk * bez_d8_i5(x) + \
97 y6 * tk * bez_d8_i6(x) + y7 * tk * bez_d8_i7(x) + y8 * tk * bez_d8_i8(x)
98
99 #
100 # Given t value between zero and one, the airfoild curve is defined as
101 #
102 # c(t) = mean(t1(t)) +/- z(t2(t)) n(t1(t)),
103 #
104 # where n is the unit normal to the mean line. See the above paper for more.
105 #
106 # Unfortunately, the parametrization of c(t) is not the same for mean(t1)
107 # and z(t2). The mean line (and its normal) can assume linear function t1 = t,
108 # -1
109 # but the thickness z_y is, in fact, a function of z_x (t). Since it is
110 # not obvious how to compute this inverse function analytically, we instead
111 # replace t in c(t) equation above by z_x(t) to get:
112 #
113 # c(z_x(t)) = mean(z_x(t)) +/- z(t) n(z_x(t)),
114 #
115 # and compute and display this instead. Note we also ignore n(t) and assumes
116 # n(t) is constant in the y direction,
117 #
118
119 airfoil_y1(t, thick) = mean_y(z_x(t)) + z_y(t, thick)
120 airfoil_y2(t, thick) = mean_y(z_x(t)) - z_y(t, thick)
121 airfoil_y(t) = mean_y(z_x(t))
122 airfoil_x(t) = mean_x(z_x(t))
123 unset grid
124 unset zeroaxis
125 set parametric
126 set xrange [-0.1:1.1]
127 set yrange [-0.1:.7]
128 set trange [ 0.0:1.0]
129 set title "NACA6409 Airfoil"
130 plot airfoil_x(t), airfoil_y(t) title "mean line" w l lt 2, \
131 airfoil_x(t), airfoil_y1(t, thick) title "upper surface" w l lt 1, \
132 airfoil_x(t), airfoil_y2(t, thick) title "lower surface" w l lt 1
133 pause -1 "Press Return"
134 mm = 0.0
135 pp = .5
136 thick = .12
137 set title "NACA0012 Airfoil"
138 set xlabel "12% thick, no camber -- classical test case"
139 plot airfoil_x(t), airfoil_y(t) title "mean line" w l lt 2, \
140 airfoil_x(t), airfoil_y1(t, thick) title "upper surface" w l lt 1, \
141 airfoil_x(t), airfoil_y2(t, thick) title "lower surface" w l lt 1
142 pause -1 "Press Return"
143 set title ""
144 set xlab ""
145 set key box
146 set parametric
147 set samples 100
148 set isosamples 10
149 set style data lines
150 set style function lines
151 print "Joukowski Airfoil using Complex Variables"
152 set title "Joukowski Airfoil using Complex Variables" offset 0,0
153 set time
154 set yrange [-.2 : 1.8]
155 set trange [0: 2*pi]
156 set xrange [-.6:.6]
157 zeta(t) = -eps + (a+eps)*exp(t*{0,1})
158 eta(t) = zeta(t) + a*a/zeta(t)
159 eps = 0.06
160 a =.250
161 set xlabel "eps = 0.06 real"
162 plot real(eta(t)),imag(eta(t))
163 pause -1 "Press Return"
164 eps = 0.06*{1,-1}
165 set xlabel "eps = 0.06 + i0.06"
166 plot real(eta(t)),imag(eta(t))
167 pause -1 "Press Return"
168 reset