"""Two ideal velocity profiles with the same peak but different area means.
Python 3.14, NumPy 2.3.5, Matplotlib 3.10.7; output in current directory.
"""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
s = np.linspace(0, 1, 2001)
r = np.linspace(-1, 1, 2001)
fig, axs = plt.subplots(1, 2, figsize=(8, 3.7), dpi=100, facecolor='white')
for n, color, name in [(2, '#256a97', 'Parabolic profile'), (8, '#b45b18', 'Blunter profile')]:
    v = 1 - s**n
    area_mean = n/(n+2)
    assert abs(np.trapezoid(2*s*v, s) - area_mean) < 1e-6
    axs[0].plot(r, 1 - abs(r)**n, color=color, lw=2, label=f'{name}: mean {area_mean:g}')
    axs[1].plot(s, 2*s*v, color=color, lw=2)
    axs[1].fill_between(s, 2*s*v, color=color, alpha=.12)
axs[0].set(xlabel='Signed radial position / radius', ylabel='Velocity / peak velocity', xlim=(-1,1), ylim=(0,1.12), title='Equal peak, different shape')
axs[0].legend(fontsize=9, loc='lower center')
axs[1].set(xlabel='Radius fraction s = r/R', ylabel='Annular contribution: 2s × velocity / peak', xlim=(0,1), ylim=(0,1.5), title='Area weighting changes the mean')
axs[1].text(.05, 1.42, 'Integrated contributions: 0.5 and 0.8', fontsize=9)
for ax in axs:
    ax.grid(alpha=.16)
    ax.spines[['top','right']].set_visible(False)
fig.suptitle('Same centerline velocity does not determine volume flow', fontsize=12)
fig.text(.5, .025, 'Constructed profiles • circular opening • not fits to patient data', ha='center', fontsize=9)
fig.tight_layout(rect=(0,.075,1,.92))
fig.savefig('doppler-profile-peak-to-mean.png', facecolor='white', dpi=100)
plt.close(fig)