"""Conditioning and truncation bias in a sampled leading-edge radius estimate.

Run from the desired media output directory. NumPy and Matplotlib dependencies
are supplied by pyproject.toml. The analytic local nose is stipulated, not a NACA
profile or measured data; the PNG is written to the working directory.
"""
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np

TRUE_RADIUS_M = 0.003
QUARTIC_COEFFICIENT = 1.0e7  # m^-3
HORIZONTAL_ERROR_M = 1.0e-5


def estimates(height_m):
    x = height_m**2 / (2*TRUE_RADIUS_M) + QUARTIC_COEFFICIENT*height_m**4
    centre = height_m**2 / (2*x)
    lower = height_m**2 / (2*(x+HORIZONTAL_ERROR_M))
    upper = height_m**2 / (2*(x-HORIZONTAL_ERROR_M))
    return x, centre, lower, upper


def main():
    height = np.linspace(0.35e-3, 2.0e-3, 900)
    x, centre, lower, upper = estimates(height)
    assert np.all(x > HORIZONTAL_ERROR_M)
    fig, ax = plt.subplots(figsize=(7.2, 4.6), dpi=100, facecolor='white')
    ax.set_facecolor('white')
    ax.fill_between(height*1000, lower/TRUE_RADIUS_M, upper/TRUE_RADIUS_M,
                    color='#98c4e8', alpha=.65, label='Bounded horizontal error: ±10 μm')
    ax.plot(height*1000, centre/TRUE_RADIUS_M, color='#145d98', lw=2.3,
            label='Quadratic quotient, exact coordinates')
    ax.axhline(1.0, color='black', lw=1.3, ls='--', label='True nose radius: 3 mm')
    ax.set(xlim=(.35, 2.0), ylim=(.6, 2.05), xlabel='Sample height from nose centreline (mm)',
           ylabel='Estimated radius / true radius')
    ax.set_title('Leading-edge radius: coordinate error versus model bias', fontsize=12)
    ax.text(.98, .94, r'$x=y^2/(2r_0)+b y^4$'+'\n'+r'$b=10^7\ \mathrm{m}^{-3}$',
            transform=ax.transAxes, ha='right', va='top', fontsize=11)
    ax.legend(loc='upper center', bbox_to_anchor=(.52,.74), fontsize=8.7, framealpha=1)
    ax.annotate('Small offsets amplify\nfixed coordinate errors', xy=(.4,1.48),
                xytext=(.72,1.85), fontsize=9,
                arrowprops={'arrowstyle':'->','color':'#555555'}, color='#444444')
    ax.annotate('Larger offsets retain\nmore quartic-model bias', xy=(1.85,.83),
                xytext=(1.15,.67), fontsize=9,
                arrowprops={'arrowstyle':'->','color':'#555555'}, color='#444444')
    ax.grid(alpha=.16)
    fig.text(.5,.017,'Original local geometry; exact vertical coordinates; not measured airfoil data.',
             ha='center',fontsize=8.5,color='#555555')
    fig.tight_layout(rect=(0,.035,1,1))
    fig.savefig(Path('leading-edge-radius-window.png'),dpi=100,facecolor='white',transparent=False)
    plt.close(fig)


if __name__ == '__main__':
    main()