;; My Polynomial fit to data from Clegg et al. (1999)
;; Returns H-alpha line emissivity in units of erg cm^3 s^-1
;; By John Wise (April, 2008)

function balmer_emission, nelec, temp

;; The fit is only good in nelec = [1e2,1e9].  Limit the data. 
local_nelec = ((nelec < 1e9) > 1e2)
log_nelec = ALOG10(local_nelec)
log_temp = ALOG10(temp)

log_nelec2 = log_nelec^2
log_temp2 = log_temp^2

a0 = -19.7929  - 1.28249*log_nelec   + 0.198715*log_nelec2
a1 = -1.34084  + 0.61503*log_nelec   - 0.097026*log_nelec2
a2 = +0.044186 - 0.0736903*log_nelec + 0.0118072*log_nelec2

; erg cm^0 s^-1 (weighted by Ne in the table)
emission = 10d0^(a0 + a1*log_temp + a2*log_temp2)

; Divide by nelec^2 to get erg cm^3 s^-1
emission /= nelec

return, 1d21*emission

end