cubic spline interpolation for airfoil rather than linear. also a sma…

…ll amount of smoothing

acmultiple.jl

@@ -231,7 +231,7 @@ type Turbine
     twist::Float64
     delta::Float64
     B::Int64
-    af::AirfoilData
+    af
     Omega::Float64
     centerX::Float64
     centerY::Float64
@@ -272,8 +272,9 @@ function radialforce(uvec::Array{Float64,1}, vvec::Array{Float64,1}, thetavec::A
     # Re = rho*W*chord/mu  # currently no Re dependence
 
     # airfoil
-    cl = turbine.af.cl[alpha]
-    cd = turbine.af.cd[alpha]
+    # cl = turbine.af.cl[alpha]
+    # cd = turbine.af.cd[alpha]
+    cl, cd = turbine.af(alpha)
 
     # rotate force coefficients
     cn = cl.*cos(phi) + cd.*sin(phi)

airfoilread.jl

@@ -2,18 +2,19 @@
 
 # export AirfoilData, readaerodyn
 
-import Interpolations: interpolate, Gridded, Linear, GriddedInterpolation
+# import Interpolations: interpolate, Gridded, Linear, GriddedInterpolation
+import Dierckx: Spline1D, evaluate
 
 # type AirfoilData
 #     alpha::Array{Float64,1}
 #     cl::Array{Float64,1}
 #     cd::Array{Float64,1}
 # end
 
-type AirfoilData
-    cl::GriddedInterpolation
-    cd::GriddedInterpolation
-end
+# type AirfoilData
+#     cl::GriddedInterpolation
+#     cd::GriddedInterpolation
+# end
 
 function readaerodyn(filename)
     """currently only reads one Reynolds number if multiple exist"""
@@ -45,18 +46,35 @@ function readaerodyn(filename)
     # af = AirfoilData(alpha*pi/180.0, cl, cd)
 
     # 1D interpolations for now.  ignoring Re dependence (which is very minor)
-    afcl = interpolate((alpha*pi/180.0,), cl, Gridded(Linear()))
-    afcd = interpolate((alpha*pi/180.0,), cd, Gridded(Linear()))
-    af = AirfoilData(afcl, afcd)
+    # afcl = interpolate((alpha*pi/180.0,), cl, Gridded(Linear()))
+    # afcd = interpolate((alpha*pi/180.0,), cd, Gridded(Linear()))
+    # af = AirfoilData(afcl, afcd)
+
+    afcl = Spline1D(alpha*pi/180, cl, s=0.1)
+    afcd = Spline1D(alpha*pi/180, cd, s=0.001)
+
+    function af(alpha)
+
+        cl = evaluate(afcl, alpha)
+        cd = evaluate(afcd, alpha)
+
+        return cl, cd
+    end
 
-    return af
+    return af#, alpha*pi/180, cl, cd
 end
 
 # end
 
-# filename = "airfoils/NACA_0012_mod.dat"
-# alpha, cl, cd = readaerodyn(filename)
+# filename = "airfoils/naca0015-wt.dat"
+# af, a0, cl0, cd0 = readaerodyn(filename)
+# alpha = linspace(-pi, pi, 2000)
+# cl, cd = af(alpha)
 # using PyPlot
+# figure()
 # plot(alpha, cl)
+# plot(a0, cl0, "o")
+# figure()
 # plot(alpha, cd)
+# plot(a0, cd0, "o")
 # show()