import math


def naca4(number, x, finite_TE = False):
    """
    Gives boundaries of the NACA 4-digit object.
    
    Parameters
    ----------
    number : str
        4-digit number specifying profile
        options: ['0009', '0012', '2412', 2414', '2415', '6409' , '0006', '0008', '0010', '0012', '0015']
    x : 1D-ndarray
        all x coordinates of the NACA profile
    finite_TE : Boolean
        True for finite thick TE
        False for zero thick TE
    
    Returns
    -------
    X : 1D ndarray
        x-coordinates of boundary NACA profile
    Y : 1D ndarray
        y-coordinates of boundary NACA profile
        
          
    """
    
    m = float(number[0])/100.0
    p = float(number[1])/10.0
    t = float(number[2:])/100.0

    a0 = +0.2969
    a1 = -0.1260
    a2 = -0.3516
    a3 = +0.2843

    if finite_TE:
        a4 = -0.1015 # For finite thick TE
    else:
        a4 = -0.1036 # For zero thick TE


    yt = [5*t*(a0*math.sqrt(xx)+a1*xx+a2*math.pow(xx,2)+a3*math.pow(xx,3)+a4*math.pow(xx,4)) for xx in x]

    xc1 = [xx for xx in x if xx <= p]
    xc2 = [xx for xx in x if xx > p]

    if p == 0:
        xu = x
        yu = yt

        xl = x
        yl = [-xx for xx in yt]

        xc = xc1 + xc2
        zc = [0]*len(xc)
    else:
        yc1 = [m/math.pow(p,2)*xx*(2*p-xx) for xx in xc1]
        yc2 = [m/math.pow(1-p,2)*(1-2*p+xx)*(1-xx) for xx in xc2]

        zc = yc1 + yc2

        dyc1_dx = [m/math.pow(p,2)*(2*p-2*xx) for xx in xc1]
        dyc2_dx = [m/math.pow(1-p,2)*(2*p-2*xx) for xx in xc2]
        dyc_dx = dyc1_dx + dyc2_dx

        theta = [math.atan(xx) for xx in dyc_dx]

        xu = [xx - yy * math.sin(zz) for xx,yy,zz in zip(x,yt,theta)]
        yu = [xx + yy * math.cos(zz) for xx,yy,zz in zip(zc,yt,theta)]

        xl = [xx + yy * math.sin(zz) for xx,yy,zz in zip(x,yt,theta)]
        yl = [xx - yy * math.cos(zz) for xx,yy,zz in zip(zc,yt,theta)]

    X = xu[::-1] + xl[1:]
    Z = yu[::-1] + yl[1:]

    return X,Z