import numpy as np
import matplotlib.pyplot as plt

def planck(erg: float | np.ndarray, kT: float = 1.0) -> float | np.ndarray:
    """ Compute the Planck distribution. """
    return erg * np.exp(-erg * kT) / kT


def compare_left_riemann_sum(num_pts: int, upper_bound: int = 8):
    """ Display rectangles used in Riemann sum. """
    x = np.linspace(0, upper_bound, num_pts)
    approx = planck(x)
    # add the 'continuous' version by selecting more points & using plot()
    plt.plot(
        np.linspace(0, upper_bound, 200), 
        planck(np.linspace(0, upper_bound, 200)),
        label="planck"
    )
    plt.bar(
        x,
        approx,
        width=upper_bound / num_pts,
        align="edge",
        alpha=0.2,
        edgecolor="b",
        label="riemann"
    )
    plt.show()


for num in (50, 10, 3):
    compare_left_riemann_sum(num)