def susceptibility_Losc(w, w0, f, w_p, y0, eps_b=1.0):
    """Dielectric Susceptibility for the Lorentzian oscillator model
    f - oscillator strength (unitless)
    The following should all use the same units i.e Hz or meV, natural frequency
    or real frequency ...
    w_p - plasma frequency 
    w - frequency axis
    w0 - transition frequency
    y0 - transition broadening
    
    eps_b is the background dielectric constant, it is included because w_p usually
    includes eps_b in its definition and so it is necessary to include this factor
    for this formula to calculate the susceptibility. Only really interesting when using
    a frequency dependent background dielectric constant, as it slightly affects 
    how we calculate the 'total' dielectric constant. It's the difference between 
    having to using the result from this function as 
    eps_b = eps_b0.(1 + Xi + Xi2...) 
    or (by including eps_b!=1.0 here) as
    eps_b = epsb0 + Xi + Xi2 +...  
    """
    chi = eps_b * w_p ** 2 * f / (w0 ** 2 - w ** 2 - 2.0j * y0 * w)
    return chi
