import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import gaussian_filter1d
from scipy.fft import fft, fftfreq

n = 100
n_pad = 1024
x = np.arange(0,n)
y = np.zeros(n)
y[n//2:] = 1
esf = gaussian_filter1d(y, sigma=10)
lsf = np.gradient(esf)
lsf_padded = np.pad(lsf, n_pad, mode='constant', constant_values=0)
mtf = np.abs(fft(lsf_padded))
mtf /= mtf[0]  # Normalize
freq = fftfreq(len(lsf_padded))
half = len(freq) // 2
freq = freq[:half]
mtf = mtf[:half]

fig, axs = plt.subplots(1, 3, constrained_layout=True, figsize=(6, 2.5))

ax = axs[0]
ax.plot(x,esf)
ax.set_title("ESF")
ax.set_xlabel("Pixel #")
ax.set_ylabel("Pixel count")
ax.set_yticks([])

ax = axs[1]
ax.plot(x,lsf)
ax.set_title("LSF")
ax.set_xlabel("Pixel #")
ax.set_ylabel('$\mathit{dESF/dx}$')
ax.set_yticks([])

ax = axs[2]
ax.plot(freq,mtf)
ax.set_title("MTF")
ax.set_xlabel("Frequency [cycles/pixel]")
ax.set_ylabel(r"$\mathcal{F}[LSF]$")
ax.set_xlim(0,0.07)

plt.show()