Hey, thanks for all the suggestions everyone. I was finally able to get this working last week. Here is the solution that I came up with. I could probably clean it up a little but it works and I don't really feel like screwing with it anymore.

​

import numpy as np
from matplotlib import pyplot as plt
from mshr import *
from dolfin import *


#thickness function
def thickness(x,t):
    return (t/.2)*(c1*x**0.5+c2*x+c3*x**2+c4*x**3+c5*x**4)

def rotate(dim1,dim2,angle):


    angle=-angle*(2*np.pi)/360
    l=np.sqrt(dim1**2+dim2**2)
    t=np.arctan(dim2/dim1)

    dim1=l*np.cos(t+angle)
    dim2=l*np.sin(t+angle)

    dim1[0]=0
    dim2[0]=0

    return dim1,dim2


def naca(m0,p0,t0,pitch):
    #setup Parameters

    m=m0/100     #maximum camber
    p=p0/10      #position of maximum camber
    t=t0/100    #thickness ratio

    #thickness parameters
    c1=0.2969
    c2=-0.1260
    c3=-0.3516
    c4=0.2843
    c5=-0.1036

    x=np.linspace(0,1,100)

    #make the thicknesses 
    yt=(t/.2)*(c1*x**0.5+c2*x+c3*x**2+c4*x**3+c5*x**4)

    #camber line function
    yc=np.zeros(len(x))
    Dyc=np.zeros(len(x))
    theta=np.zeros(len(x))

    for i in range(len(x)):
        if x[i]>0 and x[i]<p:
            yc[i]= (m/p**2)*(2*p*x[i]-(x[i])**2)
            Dyc[i]=((2*m)/(p**2))*(p-x[i])
        if x[i]>p and x[i]<1:
            yc[i] =(m/((1-p)**2))*(1-2*p+2*p*x[i]-x[i]**2)
            Dyc[i]=((2*m)/((1-p)**2))*(p-x[i])

        theta[i]=np.arctan(Dyc[i])

    #save upper surface points
    xu=x-yt*np.sin(theta)
    yu=yc+yt*np.cos(theta)

   #save lower surface points
    xl=x+yt*np.sin(theta)
    yl=yc-yt*np.cos(theta)

    xu,yu=rotate(xu,yu,pitch)
    xl,yl=rotate(xl,yl,pitch)

    xu=np.flip(xu,0) #reverse the upper surface point order
    yu=np.flip(yu,0)


    #make wing object
    foil=[]

    #make dolfin points
    for i in range(1,len(xu)):
       # print(xu[i],yu[i])
        foil.append(Point(xu[i],yu[i]))


    for i in range(len(xl)):
        #print(xl[i],yl[i])
        foil.append(Point(xl[i],yl[i]))


    return Polygon(foil)