import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt

@np.vectorize
def chapman(x, theta):
    f = lambda l, x, theta: np.exp(-(np.sqrt(x**2 + l**2 + 2*l*x*np.cos(theta)) - x))
    result, err = scipy.integrate.quad(f, 0, np.inf, args=(x, theta))
    return result

th = np.linspace(0, 90, 200)

fig, ax = plt.subplots(dpi=200, figsize=(3,3))
ax.set(yscale="log", ylim=(1, 3e2), xlim=(0, 90), xlabel="$z$", ylabel=r"$\mathrm{ch}(x,z)$")
ax.xaxis.set_major_formatter('{x:.0f}°')

for x in np.logspace(4, 0, 5):
    ch = chapman(x, np.deg2rad(th))
    ax.plot(th, ch, label=f"$x = 10^{{{np.log10(x):.0f}}}$", lw=2)

ax.legend(loc="upper left", fontsize="small")

fig.tight_layout()
fig.savefig("chapman_function.svg")