# set terminal svg size 600,400 dynamic enhanced fname 'arial' fsize 10 mousing name "airfoil_3" butt dashlength 1.0
# set output 'airfoil.3.svg'
set bar 1.000000 front
set style circle radius graph 0.02, first 0.00000, 0.00000
set style ellipse size graph 0.05, 0.03, first 0.00000 angle 0 units xy
set dummy t, y
set key inside right top vertical Right noreverse enhanced autotitle box lt black linewidth 1.000 dashtype solid
set style textbox transparent margins 1.0, 1.0 border
unset logscale
set parametric
set style data lines
unset paxis 1 tics
unset paxis 2 tics
unset paxis 3 tics
unset paxis 4 tics
unset paxis 5 tics
unset paxis 6 tics
unset paxis 7 tics
set title "Joukowski Airfoil using Complex Variables"
set timestamp "%a %b %d %H:%M:%S %Y"
set trange [ 0.00000 : 6.28319 ] noreverse nowriteback
set xlabel "eps = 0.06 real"
set xrange [ -0.600000 : 0.600000 ] noreverse nowriteback
set yrange [ -0.200000 : 1.80000 ] noreverse nowriteback
set paxis 1 range [ * : * ] noreverse nowriteback
set paxis 2 range [ * : * ] noreverse nowriteback
set paxis 3 range [ * : * ] noreverse nowriteback
set paxis 4 range [ * : * ] noreverse nowriteback
set paxis 5 range [ * : * ] noreverse nowriteback
set paxis 6 range [ * : * ] noreverse nowriteback
set paxis 7 range [ * : * ] noreverse nowriteback
set colorbox vertical origin screen 0.9, 0.2, 0 size screen 0.05, 0.6, 0 front noinvert bdefault
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
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)
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)
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))
zeta(t) = -eps + (a+eps)*exp(t*{0,1})
eta(t) = zeta(t) + a*a/zeta(t)
mm = 0.0
thick = 0.12
pp = 0.5
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.75
y7 = 0.30281
x8 = 1.0
y8 = 0.0105
GPFUN_bez_d4_i0 = "bez_d4_i0(x) = (1 - x)**4"
GPFUN_bez_d4_i1 = "bez_d4_i1(x) = 4 * (1 - x)**3 * x"
GPFUN_bez_d4_i2 = "bez_d4_i2(x) = 6 * (1 - x)**2 * x**2"
GPFUN_bez_d4_i3 = "bez_d4_i3(x) = 4 * (1 - x)**1 * x**3"
GPFUN_bez_d4_i4 = "bez_d4_i4(x) = x**4"
GPFUN_bez_d8_i0 = "bez_d8_i0(x) = (1 - x)**8"
GPFUN_bez_d8_i1 = "bez_d8_i1(x) = 8 * (1 - x)**7 * x"
GPFUN_bez_d8_i2 = "bez_d8_i2(x) = 28 * (1 - x)**6 * x**2"
GPFUN_bez_d8_i3 = "bez_d8_i3(x) = 56 * (1 - x)**5 * x**3"
GPFUN_bez_d8_i4 = "bez_d8_i4(x) = 70 * (1 - x)**4 * x**4"
GPFUN_bez_d8_i5 = "bez_d8_i5(x) = 56 * (1 - x)**3 * x**5"
GPFUN_bez_d8_i6 = "bez_d8_i6(x) = 28 * (1 - x)**2 * x**6"
GPFUN_bez_d8_i7 = "bez_d8_i7(x) = 8 * (1 - x) * x**7"
GPFUN_bez_d8_i8 = "bez_d8_i8(x) = x**8"
m0 = 0.0
m1 = 0.1
m2 = 0.1
m3 = 0.1
m4 = 0.0
GPFUN_mean_y = "mean_y(t) = m0 * mm * bez_d4_i0(t) + \t m1 * mm * bez_d4_i1(t) + \t m2 * mm * bez_d4_i2(t) + \t m3 * mm * bez_d4_i3(t) + \t m4 * mm * bez_d4_i4(t)"
p0 = 0.0
p1 = 0.2
p2 = 0.4
p3 = 0.7
p4 = 1.0
GPFUN_mean_x = "mean_x(t) = p0 * bez_d4_i0(t) + \t p1 * bez_d4_i1(t) + \t p2 * bez_d4_i2(t) + \t p3 * bez_d4_i3(t) + \t p4 * bez_d4_i4(t)"
GPFUN_z_x = "z_x(x) = x0 * bez_d8_i0(x) + x1 * bez_d8_i1(x) + x2 * bez_d8_i2(x) + \t x3 * bez_d8_i3(x) + x4 * bez_d8_i4(x) + x5 * bez_d8_i5(x) + \t x6 * bez_d8_i6(x) + x7 * bez_d8_i7(x) + x8 * bez_d8_i8(x)"
GPFUN_z_y = "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)"
GPFUN_airfoil_y1 = "airfoil_y1(t, thick) = mean_y(z_x(t)) + z_y(t, thick)"
GPFUN_airfoil_y2 = "airfoil_y2(t, thick) = mean_y(z_x(t)) - z_y(t, thick)"
GPFUN_airfoil_y = "airfoil_y(t) = mean_y(z_x(t))"
GPFUN_airfoil_x = "airfoil_x(t) = mean_x(z_x(t))"
eps = 0.06
a = 0.25
GPFUN_zeta = "zeta(t) = -eps + (a+eps)*exp(t*{0,1})"
GPFUN_eta = "eta(t) = zeta(t) + a*a/zeta(t)"
plot real(eta(t)),imag(eta(t))