"""Exact numerical frequency and accumulated phase of an undamped contact mode.

Python 3.14, NumPy and Matplotlib; see pyproject.toml. Writes its PNG to cwd.
"""
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt


def main():
    z = np.linspace(0.001, 1.9, 800)
    central = 2 * np.arcsin(z / 2) / z - 1
    trapezoidal = 2 * np.arctan(z / 2) / z - 1
    step = np.geomspace(1e-6, 1e-3, 600)
    omega, duration = 1000.0, 0.1
    phase_c = duration * (2 * np.arcsin(step * omega / 2) / step - omega)
    phase_t = duration * (omega - 2 * np.arctan(step * omega / 2) / step)
    fig, axes = plt.subplots(1, 2, figsize=(9, 4.25), dpi=120, layout='constrained')
    fig.patch.set_facecolor('white')
    axes[0].plot(z, 100 * central, label='Central difference', color='#176da5')
    axes[0].plot(z, 100 * trapezoidal, label='Trapezoidal', color='#bf5726')
    axes[0].axhline(0, color='#777777', linewidth=0.7)
    axes[0].set(xlabel=r'Dimensionless step $h\omega$', ylabel='Numerical frequency error (%)',
                title='Both curves are within their stable ranges')
    axes[1].loglog(step * 1e6, phase_c, label='Central difference', color='#176da5')
    axes[1].loglog(step * 1e6, phase_t, label='Trapezoidal', color='#bf5726')
    axes[1].axhline(0.01, color='#555555', linestyle='--', label='0.01 rad phase target')
    axes[1].set(xlabel='Step h (microseconds)', ylabel='Absolute accumulated phase error (rad)',
                title='1000 rad/s mode over 0.1 s')
    for axis in axes:
        axis.grid(alpha=0.2, which='both')
        axis.legend(fontsize=8)
    fig.suptitle('Stable contact-mode integration can still drift in phase', fontsize=13)
    fig.savefig(Path('contact-mode-phase-accumulation.png'), facecolor='white', transparent=False)
    plt.close(fig)


if __name__ == '__main__':
    main()