"""Original NACA 2412 normal-offset evaluation; writes white PNG in cwd."""
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt


def mean_line(x,m=.02,p=.4):
    x=np.asarray(x)
    if m==0:return np.zeros_like(x),np.zeros_like(x)
    assert 0<p<1
    y=np.where(x<=p,m/p**2*(2*p*x-x*x),m/(1-p)**2*(1-2*p+2*p*x-x*x))
    slope=np.where(x<=p,2*m/p**2*(p-x),2*m/(1-p)**2*(p-x))
    return y,slope


def coordinates(x):
    yc,slope=mean_line(x)
    yt=.6*(.2969*np.sqrt(x)-.126*x-.3516*x*x+.2843*x**3-.1015*x**4)
    theta=np.arctan(slope)
    return x-yt*np.sin(theta),yc+yt*np.cos(theta),x+yt*np.sin(theta),yc-yt*np.cos(theta),yc,yt,theta


def main():
    x=(1-np.cos(np.linspace(0,np.pi,1001)))/2
    xu,yu,xl,yl,yc,yt,theta=coordinates(x)
    np.testing.assert_allclose((xu+xl)/2,x,atol=1e-15)
    np.testing.assert_allclose((yu+yl)/2,yc,atol=1e-15)
    np.testing.assert_allclose(np.hypot(xu-xl,yu-yl),2*yt,atol=1e-15)
    # Offset is perpendicular to the mean-line tangent.
    np.testing.assert_allclose((xu-x)+(yu-yc)*np.tan(theta),0,atol=1e-15)
    a=coordinates(np.array([.4]))
    np.testing.assert_allclose([a[1][0],a[3][0]],[.07803010847647902,-.038030108476479016])
    np.testing.assert_allclose(mean_line(np.array([.4]))[1],0)
    np.testing.assert_allclose(mean_line(x,m=0,p=0)[0],0)
    assert np.isclose(.0005**2/(2*.00004),.003125)
    assert np.isclose(.001**2/(2*.00021),.002380952380952381)
    fig,axes=plt.subplots(1,2,figsize=(9.6,3.6),dpi=100,facecolor='white',gridspec_kw={'width_ratios':[1.5,1]})
    ax=axes[0]
    ax.plot(xu,yu,color='#2475ad',lw=2,label='Airfoil surfaces')
    ax.plot(xl,yl,color='#2475ad',lw=2)
    ax.plot(x,yc,color='#c47c22',ls='--',lw=1.8,label='Mean camber line')
    ax.axhline(0,color='#555555',ls=':',lw=1,label='Chord reference')
    ax.scatter([.4,.4],[a[1][0],a[3][0]],color='#b43f32',s=25,zorder=3)
    ax.set(xlim=(-.025,1.025),ylim=(-.08,.13),xlabel='Chordwise coordinate x/c',ylabel='Surface coordinate y/c',title='NACA 2412, normal-offset surfaces')
    ax.set_aspect('equal',adjustable='box')
    axes[1].plot(x,yc,color='#c47c22',lw=2)
    axes[1].axvline(.4,color='#888888',ls='--',lw=1)
    axes[1].axhline(.02,color='#888888',ls='--',lw=1)
    axes[1].scatter([.4],[.02],color='#b43f32',s=30,zorder=3)
    axes[1].set(xlim=(0,1),ylim=(0,.025),xlabel='Mean-line parameter ξ',ylabel='Mean-line ordinate y_c/c',title='Camber function (vertical scale expanded)')
    axes[1].text(.48,.021,'m = 0.02; p = 0.4',fontsize=10)
    for axis in axes:axis.grid(alpha=.2)
    handles,labels=ax.get_legend_handles_labels()
    fig.legend(handles,labels,loc='lower center',ncol=3,fontsize=9)
    fig.tight_layout(rect=(0,.1,1,1))
    fig.savefig(Path.cwd()/(Path(__file__).stem+'.png'),facecolor='white')
    plt.close(fig)


if __name__=='__main__':main()