def debye_spectral_density(omega: float, cutoff_freq: float,
                           reorg_energy: float) -> float:

    """
    Calculates the Debye spectral density at frequency omega, with
    a given cutoff frequency omega_c and reorganisation energy lam.
    It is normalised so that if omega=omega_c and lam=1, the
    spectral density evaluates to 1. Implements an asymmetric
    spectral density, evaluating to zero for omega <= 0. It is
    given by:

    .. math::
        \\omega^2 J(\\omega_{ab})
            = f \\frac{2\\omega_c\\omega}{(\\omega_c^2
                                           + \\omega^2) \\omega^2}

    Parameters
    ----------
    omega : float
        The frequency at which the spectral density will be
        evaluated, in units of rad ps^-1. If omega <= 0, the
        spectral density evaluates to zero.
    cutoff_freq : float
        The cutoff frequency at which the spectral density
        evaluates to 1 (or the reorg_energy value if f not equal
        to 1), in units of rad ps^-1. Must be a non-negative float.
    reorg_energy : float
        The factor by which the spectral density should be scaled
        by. Should be passed in units of rad ps^-1. Must be a
        non-negative float.

    Returns
    -------
    float
        The Debye spectral density at frequency omega, in units of
        rad ps^-1.
    """

    assert cutoff_freq >= 0., (
        'The cutoff freq must be a non-negative float, in units of rad ps^-1')
    assert reorg_energy >= 0., (
        'The scaling factor must be a non-negative float, in units of rad ps^-1')

    if omega <= 0 or cutoff_freq == 0 or reorg_energy == 0:
        # Zero if omega < 0 as an asymmetric spectral density used.
        # Zero if omega = 0 or cutoff = 0 to avoid DivideByZero error.
        return 0.
    # Returned in units of rad ps^-1
    return 2 * reorg_energy * omega * cutoff_freq / (omega**2 + cutoff_freq**2)