 codex/equal-radiant-exposure-pulses.py

"""Original equal-exposure examples; write the PNG to the working directory."""
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

target_exposure = 20.0  # J / m^2
pulses = [
    (1.0, 20.0, "#0072B2", "Long pulse: 20 W/m² for 1 s"),
    (0.25, 80.0, "#D55E00", "Short pulse: 80 W/m² for 0.25 s"),
]
time = np.linspace(-0.05, 1.05, 1201)
fig, axes = plt.subplots(1, 2, figsize=(9.2, 4.1), dpi=100, facecolor="white")

for duration, irradiance, color, label in pulses:
    assert np.isclose(duration * irradiance, target_exposure)
    # Duplicate times encode vertical pulse edges without introducing a ramp.
    pulse_time = np.array([-0.05, 0.0, 0.0, duration, duration, 1.05])
    pulse_height = np.array([0.0, 0.0, irradiance, irradiance, 0.0, 0.0])
    assert np.isclose(np.trapezoid(pulse_height, pulse_time), target_exposure)
    cumulative = irradiance * np.clip(time, 0.0, duration)
    assert np.isclose(cumulative[-1], target_exposure)
    assert np.all(np.diff(cumulative) >= 0.0)
    axes[0].plot(pulse_time, pulse_height, color=color, lw=2.3, label=label)
    axes[0].fill_between(pulse_time, pulse_height, color=color, alpha=0.08)
    axes[1].plot(time, cumulative, color=color, lw=2.3)
    axes[1].plot(duration, target_exposure, "o", color=color, markersize=5)

assert np.isclose(pulses[1][1] / pulses[0][1], 4.0)
assert np.isclose(20.0 * 1e-4, 2e-3)  # 1 cm^2 patch receives 2 mJ.
energy, radius = 1e-6, 1e-3
peak = 2 * energy / (np.pi * radius**2)
disk_mean = energy * (-np.expm1(-2.0)) / (np.pi * radius**2)
assert np.isclose(peak, 0.6366197723675814)
assert np.isclose(disk_mean, 0.275231327627)
charge, responsivity, area = 0.40e-9, 0.40, 1e-6
assert np.isclose(charge / responsivity, 1e-9, rtol=1e-12, atol=0.0)
assert np.isclose(charge / (responsivity * area), 1e-3)

axes[0].set(title="Instantaneous irradiance", ylabel="Irradiance E (W/m²)", ylim=(-2, 88))
axes[1].set(
    title="Accumulated radiant exposure",
    ylabel="Radiant exposure H (J/m²)",
    ylim=(-0.5, 23),
)
axes[1].axhline(target_exposure, color="#555555", lw=0.9, ls="--", zorder=0)
axes[1].text(0.47, 21.0, "Same final value: 20 J/m²", fontsize=9)
for ax in axes:
    ax.set(xlabel="Time t (s)", xlim=(-0.05, 1.05))
    ax.set_facecolor("white")
    ax.grid(alpha=0.18)
    ax.spines[["top", "right"]].set_visible(False)
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, labels, loc="lower center", ncol=2, frameon=False, fontsize=9)
fig.suptitle("Equal radiant exposure can come from different pulse histories", fontsize=13)
fig.subplots_adjust(left=0.08, right=0.98, bottom=0.23, top=0.80, wspace=0.28)
fig.savefig(Path.cwd() / "equal-radiant-exposure-pulses.png", facecolor="white")
plt.close(fig)