From e67a786f05151438f2e2d3a109813b4f194e19c1 Mon Sep 17 00:00:00 2001
From: Michael Fairchild <mjfairch@gmail.com>
Date: Mon, 4 Jan 2016 14:24:42 -0500
Subject: [PATCH 1/2] Fix NACA airfoil.  Change numpoints to num_points for
 consistency.

The thickness adjustment must be normal, rather than vertical, to the camber
line.  For uncambered foils (camber line corresponds with x-axis), there is no
difference, but for cambered foils there is.
---
 pysces/body.py            | 190 ++++++++++++++++++++++++++++++--------
 pysces/tests/test_body.py |   8 +-
 scripts/animate_foil.py   |   7 +-
 3 files changed, 158 insertions(+), 47 deletions(-)

diff --git a/pysces/body.py b/pysces/body.py
index 7351b16..808d3a6 100644
--- a/pysces/body.py
+++ b/pysces/body.py
@@ -56,11 +56,11 @@ def flat_plate(num_points):
     y = np.zeros_like(x)
     return Body(np.array([x, y]).T)
 
-def joukowski_foil(xcenter=-.1, ycenter=.1, a=1, numpoints=32):
+def joukowski_foil(xcenter=-.1, ycenter=.1, a=1, num_points=32):
     """Return a Joukowski foil Body.
     
     The foil has its trailing edge at (2a,0).  The foil has a total of
-    numpoints along the boundary.  Refer to chapter 4 of [1]_ for details.
+    num_points along the boundary.  Refer to chapter 4 of [1]_ for details.
 
     Parameters
     ----------
@@ -74,7 +74,7 @@ def joukowski_foil(xcenter=-.1, ycenter=.1, a=1, numpoints=32):
     a : float
         radius of the Joukowski preimage circle
 
-    numpoints : int
+    num_points : int
         number of points along the boundary
 
     References
@@ -82,7 +82,7 @@ def joukowski_foil(xcenter=-.1, ycenter=.1, a=1, numpoints=32):
     .. [1] Acheson, D. J., "Elementary Fluid Dynamics", Oxford, 1990.
     """
 
-    t = np.linspace(0,2*np.pi,numpoints)
+    t = np.linspace(0,2*np.pi,num_points)
     r = np.sqrt((a-xcenter)**2+ycenter**2)
     chi = xcenter + r*np.cos(t)
     eta = ycenter + r*np.sin(t)
@@ -91,7 +91,7 @@ def joukowski_foil(xcenter=-.1, ycenter=.1, a=1, numpoints=32):
     y = eta*(1-a**2/mag2)
     return Body(np.array([x,y]).T)
 
-def karman_trefftz_foil(xcenter=-.1, ycenter=0, a=.1, angle_deg=10, numpoints=32):
+def karman_trefftz_foil(xcenter=-.1, ycenter=0, a=.1, angle_deg=10, num_points=32):
     """Return a Karman-Trefftz foil Body.
     
     The Karman-Trefftz foil is a modified version of the Joukowski
@@ -106,7 +106,7 @@ def karman_trefftz_foil(xcenter=-.1, ycenter=0, a=.1, angle_deg=10, numpoints=32
     angle_deg : float
         The interior angle, in degrees, at the trailing edge.
 
-    numpoints : int
+    num_points : int
         Number of points along the boundary
 
     See Also
@@ -120,7 +120,7 @@ def karman_trefftz_foil(xcenter=-.1, ycenter=0, a=.1, angle_deg=10, numpoints=32
 
     angle_rad = angle_deg*np.pi/180
     n = 2-angle_rad/np.pi
-    t = np.linspace(0,2*np.pi,numpoints)
+    t = np.linspace(0,2*np.pi,num_points)
     ctr = xcenter + 1j*ycenter
     r = np.linalg.norm(ctr-a)
     zeta = ctr+r*np.exp(1j*t)
@@ -131,7 +131,7 @@ def karman_trefftz_foil(xcenter=-.1, ycenter=0, a=.1, angle_deg=10, numpoints=32
     return Body(np.array([x,y]).T)
 
 def van_de_vooren_foil(semichord=1.0, thickness=0.15, angle_deg=5,
-numpoints=32):
+num_points=32):
     """Return a van de Vooren foil Body.
 
     Refer to section 6.6 of [1]_
@@ -147,7 +147,7 @@ def van_de_vooren_foil(semichord=1.0, thickness=0.15, angle_deg=5,
     angle_deg : float
         interior angle, in degrees, at the trailing edge
 
-    numpoints : int
+    num_points : int
         number of points along the boundary
 
     References
@@ -158,7 +158,7 @@ def van_de_vooren_foil(semichord=1.0, thickness=0.15, angle_deg=5,
 
     k = 2-(angle_deg*np.pi/180)
     a = 2*semichord*((1+thickness)**(k-1))*2**(-k)
-    t = np.linspace(0,2*np.pi,numpoints)
+    t = np.linspace(0,2*np.pi,num_points)
     num = (a*(np.cos(t)-1)+1j*a*np.sin(t))**k
     den = (a*(np.cos(t)-thickness)+1j*a*np.sin(t))**(k-1)
     z = (num/den)+semichord
@@ -166,42 +166,152 @@ def van_de_vooren_foil(semichord=1.0, thickness=0.15, angle_deg=5,
     y = [w.imag for w in z]
     return Body(np.array([x,y]).T)
 
-def naca_airfoil(code, num_points, zero_thick_te=False, uniform=False):
-    """Return a NACA 4-digit series airfoil"""
-    # extract parameters from 4-digit code
+def naca_airfoil(code, num_points=20, te_clamp=False, uniform=False, chord=1):
+    """Returns a NACA 4-digit series foil.
+
+    Refer to [1]_ for a detailed description and formulas.
+
+    Parameters:
+    -----------
+    code_str : string
+        The 4-digit NACA code describing the foil.  If the first digit is M, 
+        the second digit is P, and the last two digits are TH, then:
+        M = maximum camber, as a percentage of the chord
+        P = distance of max camber from leading edge in tens of percent of chord
+        TH = thickenss of the airfoil as a percentage of chord
+    
+        For example, 2315 means a cambered airfoil whose max camber is 2% of the
+        chord, located 30% chordwise from the leading edge, and whose thickness
+        is 15% of the chord.
+    
+    num_points : int
+        The number of points along each surface of the foil. There will be 
+        2*num_points-1 points returned, with num_points along the edge and 
+        num_points-1 along the bottom edge.
+    
+    te_clamp : boolean
+        Trailing edge clamp.  When false, the trailing edge has a slight 
+        nonzero thickness.  When true, the foil adjusted so that the trailing 
+        edge has exactly zero thickness.
+     
+    uniform : boolean
+        Distribution of points along chord.  When false, points are uniformly 
+        distributd along the chord.  When true, points are distributed more 
+        densly near leading edge (where the curvature is large).
+    
+    chord : double
+        Length of chord
+    
+    Returns:
+    --------
+    A two-column matrix q, whose first column encodes the x coordinates and
+    whose second column encodes the y coordinates of the points on the foil.
+    The rows are ordered so that the points traverse the foil in the
+    counterclockwise fashion from the trailing edge to the leading edge and
+    back again.
+    
+    References:
+    -----------
+    .. [1] https://en.wikipedia.org/wiki/NACA_airfoil
+    
+    Example:
+    --------
+    q = naca_foil('2315',50)
+ """
     code_str = "%04d" % int(code)
     if len(code_str) != 4:
         raise ValueError("NACA designation is more than 4 digits")
-    max_camber = 0.01 * int(code_str[0])
-    p = 0.1 * int(code_str[1])  # location of max camber
-    thickness = 0.01 * int(code_str[2:])
-    if uniform:
-        x = np.linspace(0, 1, num_points)
+    m = .01*int(code_str[0])
+    p = .1*int(code_str[1])
+    thick = .01*int(code_str[2:4])
+    if (uniform is True):
+        x = np.linspace(0,chord,num_points)
     else:
-        # closer spacing near leading edge
-        theta = np.linspace(0, 0.5 * np.pi, num_points)
-        x = 1 - np.cos(theta)
-
-    # thickness
-    coefs = [-0.1015, 0.2843, -0.3516, -0.1260, 0, 0.2969]
-    if zero_thick_te:
-        coefs[0] = -0.1036
-    y_thick = 5 * thickness * (np.polyval(coefs[:5], x) +
-                               coefs[5] * np.sqrt(x))
-
-    # camber
-    front = np.where(x <= p)
-    back = np.where(x > p)
-    y_camber = np.zeros_like(x)
+        t = np.linspace(0,0.5*np.pi,num_points)
+        x = chord*(1-np.cos(t))
+
+    # Construct the thickness line, yt
+    xc = x/chord
+    coefs = [-.1015, .2843, -.3516, -.1260, 0]
+    sqrt_coef = .2969
+    if (te_clamp is True):
+        coefs[0] = -.1036
+    yt = 5*thick*chord*(sqrt_coef*np.sqrt(xc)+np.polyval(coefs,xc))
+
+    # Avoid unnecessary work (and avoid divide-by-zero) when no camber
     if p:
-        y_camber[front] = max_camber * x[front] / p**2 * (2 * p - x[front])
-        y_camber[back] = max_camber * ((1. - x[back])/(1. - p)**2 *
-                                       (1 + x[back] - 2 * p))
-    x = np.hstack([x[-1:0:-1], x])
-    y = np.hstack([y_camber[-1:0:-1] + y_thick[-1:0:-1],
-                   y_camber - y_thick])
-    return Body(np.array([x, y]).T)
+        # Construct the camber line, yc
+        front = np.where(x <= p*chord)
+        back = np.where(x > p*chord)
+        yc = np.zeros_like(x)
+        yc[front] = (m/p**2)*x[front]*(2*p-xc[front])
+        yc[back] = (m/(1-p)**2)*(chord-x[back])*(1+xc[back]-2*p)
+
+        # Combine thickness and camber line to produce the NACA airfoil.
+        # The usual formulas for the camber adjustment involve sin(arctan(.)) 
+        # and cos(arctan(.)), but to avoid branch cuts in the arctangent 
+        # function, I use instead the equivalent Cartesian formulations:
+        # sin(arctan(x)) = x/sqrt(1+x^2)
+        # cos(arctan(x)) = 1/sqrt(1+x^2)
+        dycdx = np.zeros_like(x)
+        dycdx[front] = (2*m/p**2)*(p-xc[front])
+        dycdx[back] = (2*m/(1-p)**2)*(p-xc[back])
+        z = np.sqrt(1+dycdx**2)
+        xU = x - yt*dycdx/z # X coordinates along upper surface
+        xL = x + yt*dycdx/z # X coordinates along lower surface
+        yU = yc + yt/z      # Y coordinates along upper surface
+        yL = yc - yt/z      # Y coordinates along lower surface
+        x = np.hstack([xU[-1:0:-1],xL])
+        y = np.hstack([yU[-1:0:-1],yL])
+    else:
+        x = np.hstack([x[-1:0:-1],x])
+        y = np.hstack([yt[-1:0:-1],-yt])
+    # Note: The [-1:0:-1] indexing above accomplishes two things: (1) it
+    # reorders the vectors so the foil is described counterclockwise from the
+    # trailing edge to the leading edge and back again, and (2) it omits the
+    # repeated entry at the leading edge that would otherwise be present from
+    # concatenating two copies of x.  It is important not to have two identical
+    # successive entries, or else the matrix of influence coefficients in the
+    # panel (boundary element) method becomes ill conditioned.
+
+    return Body(np.array([x,y]).T)
 
+#def naca_airfoil(code, num_points, zero_thick_te=False, uniform=False):
+#    """Return a NACA 4-digit series airfoil"""
+#    # extract parameters from 4-digit code
+#    code_str = "%04d" % int(code)
+#    if len(code_str) != 4:
+#        raise ValueError("NACA designation is more than 4 digits")
+#    max_camber = 0.01 * int(code_str[0])
+#    p = 0.1 * int(code_str[1])  # location of max camber
+#    thickness = 0.01 * int(code_str[2:])
+#    if uniform:
+#        x = np.linspace(0, 1, num_points)
+#    else:
+#        # closer spacing near leading edge
+#        theta = np.linspace(0, 0.5 * np.pi, num_points)
+#        x = 1 - np.cos(theta)
+#
+#    # thickness
+#    coefs = [-0.1015, 0.2843, -0.3516, -0.1260, 0, 0.2969]
+#    if zero_thick_te:
+#        coefs[0] = -0.1036
+#    y_thick = 5 * thickness * (np.polyval(coefs[:5], x) +
+#                               coefs[5] * np.sqrt(x))
+#
+#    # camber
+#    front = np.where(x <= p)
+#    back = np.where(x > p)
+#    y_camber = np.zeros_like(x)
+#    if p:
+#        y_camber[front] = max_camber * x[front] / p**2 * (2 * p - x[front])
+#        y_camber[back] = max_camber * ((1. - x[back])/(1. - p)**2 *
+#                                       (1 + x[back] - 2 * p))
+#    x = np.hstack([x[-1:0:-1], x])
+#    y = np.hstack([y_camber[-1:0:-1] + y_thick[-1:0:-1],
+#                   y_camber - y_thick])
+#    return Body(np.array([x, y]).T)
+#
 
 class TransformedBody(object):
     """Base class for rigid (Euclidean) transformations of existing bodies
diff --git a/pysces/tests/test_body.py b/pysces/tests/test_body.py
index 726b1d5..ef3698d 100644
--- a/pysces/tests/test_body.py
+++ b/pysces/tests/test_body.py
@@ -19,22 +19,22 @@ def test_airfoil_uniform(self):
         npoints = 20
         body = naca_airfoil("0012", npoints, uniform=True)
         self.assertEqual(len(body.get_points()), 2*npoints - 1)
-        body2 = naca_airfoil("2412", npoints, zero_thick_te=True)
+        body2 = naca_airfoil("2412", npoints, te_clamp=True)
         self.assertEqual(len(body2.get_points()), 2*npoints - 1)
 
     def test_joukowski_foil(self):
         npts = 32
-        body = joukowski_foil(numpoints=npts)
+        body = joukowski_foil(num_points=npts)
         self.assertEqual(len(body.get_points()),npts)
 
     def test_karman_trefftz_foil(self):
         npts = 32
-        body = karman_trefftz_foil(numpoints=npts)
+        body = karman_trefftz_foil(num_points=npts)
         self.assertEqual(len(body.get_points()),npts)
 
     def test_van_de_vooren_foil(self):
         npts = 32
-        body = van_de_vooren_foil(numpoints=npts)
+        body = van_de_vooren_foil(num_points=npts)
         self.assertEqual(len(body.get_points()),npts)
 
 class TestBody(unittest.TestCase):
diff --git a/scripts/animate_foil.py b/scripts/animate_foil.py
index 191e9b6..e05a57c 100644
--- a/scripts/animate_foil.py
+++ b/scripts/animate_foil.py
@@ -4,9 +4,10 @@
 import matplotlib.animation as animation
 import sys
 
-#airfoil = naca_airfoil("0006", 20) # NACA 0012 airfoil with 20 points per side
-#airfoil = naca_airfoil("2214", 20) 
-airfoil = joukowski_foil(-.1,.1,.5,100)
+#airfoil = naca_airfoil("0006", 20)
+airfoil = naca_airfoil("2214", 20) 
+airfoil = TransformedBody(airfoil, 0, [-.25,0])
+#airfoil = joukowski_foil(-.1,.1,.5,100)
 #airfoil = van_de_vooren_foil(0.5, 0.1, 3)
 #airfoil = karman_trefftz_foil(-.1,.1,1,10,32)
 #airfoil = flat_plate(20)

From 656c96c69ffacb4c4acfe27ab5595fbfbf240534 Mon Sep 17 00:00:00 2001
From: Michael Fairchild <mjfairch@gmail.com>
Date: Mon, 4 Jan 2016 14:29:57 -0500
Subject: [PATCH 2/2] Remove commented (old) version of NACA airfoil.

---
 pysces/body.py | 37 -------------------------------------
 1 file changed, 37 deletions(-)

diff --git a/pysces/body.py b/pysces/body.py
index 808d3a6..da0d185 100644
--- a/pysces/body.py
+++ b/pysces/body.py
@@ -276,43 +276,6 @@ def naca_airfoil(code, num_points=20, te_clamp=False, uniform=False, chord=1):
 
     return Body(np.array([x,y]).T)
 
-#def naca_airfoil(code, num_points, zero_thick_te=False, uniform=False):
-#    """Return a NACA 4-digit series airfoil"""
-#    # extract parameters from 4-digit code
-#    code_str = "%04d" % int(code)
-#    if len(code_str) != 4:
-#        raise ValueError("NACA designation is more than 4 digits")
-#    max_camber = 0.01 * int(code_str[0])
-#    p = 0.1 * int(code_str[1])  # location of max camber
-#    thickness = 0.01 * int(code_str[2:])
-#    if uniform:
-#        x = np.linspace(0, 1, num_points)
-#    else:
-#        # closer spacing near leading edge
-#        theta = np.linspace(0, 0.5 * np.pi, num_points)
-#        x = 1 - np.cos(theta)
-#
-#    # thickness
-#    coefs = [-0.1015, 0.2843, -0.3516, -0.1260, 0, 0.2969]
-#    if zero_thick_te:
-#        coefs[0] = -0.1036
-#    y_thick = 5 * thickness * (np.polyval(coefs[:5], x) +
-#                               coefs[5] * np.sqrt(x))
-#
-#    # camber
-#    front = np.where(x <= p)
-#    back = np.where(x > p)
-#    y_camber = np.zeros_like(x)
-#    if p:
-#        y_camber[front] = max_camber * x[front] / p**2 * (2 * p - x[front])
-#        y_camber[back] = max_camber * ((1. - x[back])/(1. - p)**2 *
-#                                       (1 + x[back] - 2 * p))
-#    x = np.hstack([x[-1:0:-1], x])
-#    y = np.hstack([y_camber[-1:0:-1] + y_thick[-1:0:-1],
-#                   y_camber - y_thick])
-#    return Body(np.array([x, y]).T)
-#
-
 class TransformedBody(object):
     """Base class for rigid (Euclidean) transformations of existing bodies
     """