 codex/pressure-tap-aliasing.py

"""Illustrate an unsampled pressure peak; write the original PNG to cwd."""

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


def pressure_difference(x):
    x = np.asarray(x)
    inside = (x >= 0.4) & (x <= 0.6)
    bump = np.where(inside, 1.2 * np.sin(np.pi * (x - 0.4) / 0.2) ** 2, 0.0)
    return 0.8 + bump


def main():
    taps = np.linspace(0, 1, 6)
    x = np.linspace(0, 1, 10001)
    baseline = np.full_like(x, 0.8)
    peaked = pressure_difference(x)
    np.testing.assert_allclose(pressure_difference(taps), np.full(6, 0.8), atol=1e-14)
    areas = np.array([np.trapezoid(baseline, x), np.trapezoid(peaked, x)])
    np.testing.assert_allclose(areas, [0.8, 0.92], atol=1e-12)
    for values in [np.full(6, 0.8), pressure_difference(taps)]:
        np.testing.assert_allclose(np.trapezoid(values, taps), 0.8, atol=1e-14)
    np.testing.assert_allclose(500.0 * 0.2 * areas, [80.0, 92.0])

    plt.rcParams.update({"font.size": 11, "figure.facecolor": "white", "axes.facecolor": "white"})
    fig, ax = plt.subplots(figsize=(8.4, 4.2), dpi=100)
    ax.plot(x, peaked, color="#C55A11", linewidth=2.5, label="Profile with an unsampled peak")
    ax.plot(x, baseline, color="#276FBF", linewidth=2, linestyle="--", label="Uniform profile")
    ax.fill_between(x, baseline, peaked, color="#D98037", alpha=0.22)
    ax.scatter(taps, pressure_difference(taps), s=40, color="#222222", zorder=5,
               label="Six identical tap readings")
    ax.annotate("Extra integral = 0.12\nExtra normal force = 12 N/m\n(q = 500 Pa, chord = 0.20 m)",
                xy=(0.52, 1.25), xytext=(0.70, 1.54),
                arrowprops={"arrowstyle": "->", "color": "#555555"}, fontsize=10)
    ax.set_xlim(-0.03, 1.03)
    ax.set_ylim(0.55, 2.30)
    ax.set_xticks(taps)
    ax.set_xlabel("Chord coordinate ξ = x/c")
    ax.set_ylabel("ΔCₚ = Cₚ,lower − Cₚ,upper")
    ax.set_title("Same tap readings, different pressure integrals", fontsize=13, pad=12)
    ax.legend(loc="upper left", fontsize=10, framealpha=0.96)
    ax.grid(alpha=0.18)
    ax.spines[["top", "right"]].set_visible(False)
    fig.subplots_adjust(left=0.105, right=0.98, top=0.87, bottom=0.16)
    output = Path(__file__).with_suffix(".png").name
    fig.savefig(output, dpi=100, facecolor="white")
    plt.close(fig)
    print(f"{output}: 840×420; integrals {areas.tolist()}; identical six samples")


if __name__ == "__main__":
    main()