 codex/ir-diagonal-overlap-correlation.py

"""Original illustrative 2D line-shape model. Python 3.14.4, NumPy 2.3.5,
Matplotlib 3.10.7. Run from _media; writes the PNG in the current directory.
"""
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np

centres = (1600.0, 1650.0)
sigma = 20.0
frequency = np.linspace(1560, 1690, 300)
excitation, detection = np.meshgrid(frequency, frequency)

def peak(x, y, rho):
    return np.exp(-(x*x - 2*rho*x*y + y*y) / (2*sigma*sigma*(1-rho*rho)))

with plt.rc_context({'font.size': 10, 'savefig.facecolor': 'white'}):
    fig, axes = plt.subplots(1, 2, figsize=(9, 4.2), dpi=100, facecolor='white')
    fig.subplots_adjust(left=.08, right=.88, bottom=.19, top=.80, wspace=.26)
    for ax, rho in zip(axes, (0.0, .8)):
        signal = sum(peak(excitation-c, detection-c, rho) for c in centres)
        plotted = ax.contourf(excitation, detection, signal, levels=np.linspace(0, 1.05, 43), cmap='Blues', extend='max')
        ax.plot([1560, 1690], [1560, 1690], linestyle='--', linewidth=.8, color='#777777')
        ax.plot([1600], [1650], marker='x', markersize=7, color='#C43B28', markeredgewidth=1.8)
        cross_value = 2*peak(0.0, centres[1]-centres[0], rho)
        ax.set(title=f'Correlation ρ = {rho:g}', xlabel='Excitation wavenumber (cm⁻¹)', aspect='equal')
        ax.set_xticks([1575, 1625, 1675]);ax.set_yticks([1575, 1625, 1675])
        ax.text(.5, -.25, f'S(1600, 1650) = {cross_value:.5f}', transform=ax.transAxes, ha='center')
    axes[0].set_ylabel('Detection wavenumber (cm⁻¹)')
    cax = fig.add_axes([.91, .23, .018, .52]);fig.colorbar(plotted, cax=cax);cax.set_title('Amplitude', fontsize=8);cax.tick_params(labelsize=8)
    fig.suptitle('Same marginal width, different off-diagonal tails', y=.97, fontsize=13)
    fig.text(.5, .88, 'Two diagonal features only; red × marks an off-diagonal coordinate', ha='center')
    fig.savefig(Path('ir-diagonal-overlap-correlation.png'), dpi=100)
    plt.close(fig)
assert abs(2*peak(0,50,0)-.087873867)<1e-8
assert abs(2*peak(0,50,.8)-.000339713)<1e-8