!------------------------------------------------------------------------------
!AIRFOIL BUILDER SUBROUTINE by Chris Teubert
!
!This Subroutine takes an input of a NACA 4-digit airfoil mptt and the Number
!of points to split the airfoil into and builds arrays of x and y values of the
!points on the airfoil
!
!Inputs: N, m, p, tt
!Outputs: x(N), y(N)
!~Christopher Teubert
!------------------------------------------------------------------------------

SUBROUTINE AIRFOIL_BUILDER(N,m,p,tt,x,y)
IMPLICIT NONE

INTEGER                                 ::      N
REAL (KIND=8)                           ::      m, p, tt
REAL (KIND=8)                           ::      x(N), y(N)

INTEGER                                 ::      i
REAL (KIND=8)                           ::      yc(N), yt(N), c1, c2, c3, c4, c5

!set C Values
c1 = 0.296375
c2 = -0.12635
c3 = -0.35195
c4 = 0.283775
c5 = -0.10185

!adjusting...
m = m / 100
p = p / 10
tt = tt / 100

!Creates Airfoil
DO i=1, N
        IF (real(i) < N/2) THEN !Sets x values
                x(i) = 1. - (i-1.) * (1/(N/2))
        ELSE
                x(i) = (i-N/2)*(1/(N/2))
        END IF

        IF (x(i)<p) THEN !Sets yc (center cord
                yc(i) = m/(p*p)*(2*p*x(i)-x(i)*x(i))
        ELSE
                yc(i) = m / ((1-p)*(1-p)) * ((1-2*p) + 2*p*x(i) - x(i) * x(i))
        END IF

        !Thickness
        yt(i) = tt / .2d0 * (c1 * DSQRT(x(i)) + c2*x(i) + c3*x(i)*x(i) + c4*x(i)*x(i)*x(i) + c5*x(i)**4)

        IF (real(i) < N/2) THEN !Calculates y(i) from center cord and thickness
                y(i) = yc(i) - yt(i)
        ELSE
                y(i) = yc(i) + yt(i)
        ENDIF
END DO

END SUBROUTINE AIRFOIL_BUILDER