# smooth spline through airfoil coordinates

I am trying to produce a smooth spline curve through a set of airfoil coordinates. An example for the data I am trying to use is here, and I am reading the file like so to get rid of some extraneous stuff:

airfoildata = Drop[Take[
Drop[
Import["n64015.dat"], 3], 53], {27}];
airfoildata = Join[Reverse[Drop[Take[airfoildata, 26], 1]],
Drop[airfoildata, 26]];


ListLinePlot gives a plot of the airfoil that's alright, but you can see some "kinks" in the plot near the leading edge, where data are sparse relative to the curvature. So, I tried getting a real smooth curve via B-Splines, using this code:

pts = Take[airfoildata, 20];
n = Length[pts];
dist = Accumulate[
Table[EuclideanDistance[pts[[i]], pts[[i + 1]]], {i, Length[pts] - 1}]];
param = N[Prepend[dist/Last[dist], 0]];
deg = 3;
knots = Join[ConstantArray[0, deg], Range[0, 1, 1/(n - deg)],
ConstantArray[1, deg]];
m = Table[BSplineBasis[{deg, knots}, j - 1, param[[i]]], {i, n}, {j, n}];
ctrlpts = LinearSolve[m, pts];


and then plotting it via

ListPlot[airfoildata, Prolog -> BSplineCurve[ctrlpts],
PlotStyle -> Directive[Red, PointSize[Medium]]]


Using only the first 20 points as in the above code fragment, things look good, but if I use 21 points I can see numerical instability creeping in, for 22 points things get truly awful, and beyond that the linear system for the control points has no solution.

So here's my question: How can I obtain a smooth curve (with a nicely rounded leading edge near the origin) that respects the given airfoil coordinates? It looks like a standard B-spline will not work, and somehow this seems to be related to the fact that I have a relatively large number of points. However, I have used various CAD packages in the past that had no problems creating a smooth curve through data like the ones I have.

• Since the data you linked corresponds to a NACA airfoil, you might be interested in an implementation of the formulae describing the profile that came up previously on the site. See: Proper use of a formula to produce a cambered airfoil. – MarcoB Jul 27 '16 at 20:53

After permuting the points so that the cusp is the first point (just like in mns's answer), the methods of this answer can be used:

airfoil = Flatten[MapAt[Most @* Reverse, DeleteCases[Split[Drop[
# =!= {} && #2 =!= {} &], {{}}], 1], 1];

(* Lee's method, http://dx.doi.org/10.1016/0010-4485(89)90003-1 *)
parametrizeCurve[pts_List, a : (_?NumericQ) : 1/2] :=
FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]] /;
MatrixQ[pts, NumericQ]

tvals = parametrizeCurve[airfoil];
m = 3; (* degree of the B-spline *)
(* knots for interpolating B-spline *)
knots = Join[ConstantArray[0, m + 1], MovingAverage[ArrayPad[tvals, -1], m],
ConstantArray[1, m + 1]];
(* basis function matrix *)
bas = Table[BSplineBasis[{m, knots}, j - 1, tvals[[i]]],
{i, Length[airfoil]}, {j, Length[airfoil]}];
ctrlpts = LinearSolve[bas, airfoil];

Graphics[{{Directive[AbsoluteThickness[3], ColorData[97, 1]],
BSplineCurve[ctrlpts, SplineDegree -> 3, SplineKnots -> knots]},
{AbsolutePointSize[5], Point[airfoil]}}, Axes -> None, Frame -> True]


Replace BSplineCurve[] with BSplineFunction[] if needed.

• Very nice solution, answering my (implied) question on a "good" parameterization of the curve for the purpose of interpolation. Thanks! I'll mark this as my best answer. – Pirx Jul 29 '16 at 17:02
• The code in your question used chord-length parametrization (corresponding to setting a to 1 in parametrizeCurve[]), which should also work; the snag you hit into was that you did not permute the points so that the cusp is at the start and end. But Lee shows that centripetal parametrization is often a more "natural" choice for many configurations. – J. M. is away Jul 29 '16 at 17:06

In an airfoil you probably want a sharp trailing edge and a rounded leading edge. Therefore I would choose the trailing edge as first and last point for your spline:

path = "http://m-selig.ae.illinois.edu/ads/coord/n64015.dat";
airfoilData = Drop[Take[Drop[Import@path,3], 53], {27}];
upperSide = airfoilData[[1 ;; Length@airfoilData/2]];
lowerSide = airfoilData[[Length@airfoilData/2 + 1 ;; -1]];
airfoilDataOrdered = Reverse@upperSide~Join~lowerSide[[2 ;; -1]];


Create spline:

f = BSplineFunction[airfoilDataOrdered];


Plot:

ListPlot[{airfoilDataOrdered, f /@ Range[0, 1, 0.01]}]


EDIT 1

To go through the points you can do a parametric interpolation.

f2 = Interpolation[Table[{i, airfoilDataOrdered[[i]]}, {i, Length@airfoilDataOrdered}],
Method -> "Spline", InterpolationOrder -> 5];
ParametricPlot[f2[x], {x, 1, Length@airfoilDataOrdered},
PlotRange -> {{0, 0.03}, {-0.03, 0.03}}]


The nose should now be ok.

EDIT 2

Polar approach:

iPoints =
Table[{π - ArcTan[-airfoilDataOrdered[[i, 1]] + 1/4, airfoilDataOrdered[[i, 2]]],
airfoilDataOrdered[[i]] - {1/4, 0}}, {i, Length@airfoilDataOrdered}];
iPoints[[-1, 1]] = 2 π;
f3 = Interpolation[iPoints, Method -> "Spline", InterpolationOrder -> 2];
ParametricPlot[f3[x], {x, 0, 2 π}, PlotRange -> All]

• Hmm, this kind of works, except for the fact that this takes the airfoil coordinates as the control points of the spline. As a consequence, the resulting function does not, in general, go through the given points. For example, in your solution above you can see that the leading edge (the very front of the airfoil) now sits around x=0.0025, not at x=0 as it should. – Pirx Jul 27 '16 at 20:42
• @Pirx Answer edited with a parametric interpolation. – Gypaets Jul 27 '16 at 22:07
• Awesome, that's exactly what I was looking for. Thanks! – Pirx Jul 27 '16 at 22:15
• One more comment: Your solution above uses the index of the airfoil coordinates as the independent variable in a parametric definition of the airfoil coordinates. I had this idea that it might be better to preserve the information on the spacing of the data points, doing something like f3 = Interpolation[ Table[{Sign[26 - i] airfoildata[[i, 1]], airfoildata[[i]]}, {i, Length@airfoildata}], Method -> "Spline", InterpolationOrder -> 3]; It turns out this works also, but if you increase the interpolation order this approach breaks down completely. – Pirx Jul 27 '16 at 22:39
• @Pirx In that case keep it low. You will almost always have problems with higher Interpolation orders. I have also added a polar approach where the nose seems to be a bit more rounded than with your spacing. – Gypaets Jul 28 '16 at 7:30