function dust_extinction_data
COMMON dust_data, dust_data
    if N_ELEMENTS(dust_data) eq 0 then begin
        fn = 'milky_way_extinction.dat'
        header_lines = 13
        data_lines = FILE_LINES(fn) - header_lines
        openr, lun, fn, /GET_LUN
        header = STRARR(header_lines)
        readf, lun, header
        ;  lambda    albedo   <cos>  C_ext/H    K_abs   <cos^2>
        ; (micron)                   (cm^2/H)  (cm^2/g)          comment
        ;----------- ------  ------ --------- --------- ------- --------
        ;1.00000E+04 0.0000  0.0000 1.225E-28 6.553E-03 0.66939
        data = DBLARR(6, data_lines)
        readf, lun, data
        free_lun, lun
        data[0,0] = data[0,*] * 1d-4 ; convert microns to cm
        dust_data = data
    endif
    return, dust_data
end

function dust_sublimation_fraction, T
    ; this is entirely ad-hoc, but extends Krumholz et al.
    ; and provides a plausible temperature-continuous value.
    ; the grains' icy mantles sublimate at 110 K
    frac = (2. - 1. / (1. + exp(-(T-110.)/5.)))
    ; cross section falls by a factor of 8 above a few thousand K
    frac *= (1. - 0.875 * (1. / (1. + exp(-(T-1200.) / 100.))))
    return, frac
end

; lambda in cm, T in K
function dust_absorption_coefficient, lambda, T, n_H
    dust_data = dust_extinction_data()
    cs = INTERPOL(dust_data[3,*], dust_data[0,*], lambda, /LSQUADRATIC)
    return, n_H * dust_sublimation_fraction(T) * cs[0]
end

function dust_emissivity, alpha_nu, lambda, T
    dust_data = dust_extinction_data()
    albedo = INTERPOL(dust_data[1,*], dust_data[0,*], lambda, /LSQUADRATIC)
    ; Kirchhoff's law
    return, alpha_nu * (1-albedo[0]) * (constant('2*h*c') / lambda^3) / (exp(constant('h*c/kB') / (lambda * T)) - 1)
end

function dust_cross_section, LAMBDA=lambda, NEEDED_FIELDS=return_fields, DATA=a
    density_unit = get_unit('Electron_Density', /FORCE_USE_UNITS)

    fields = ["HII_Density", "HI_Density", "Temperature"]
    if KEYWORD_SET(return_fields) then return, fields

    out = PTRARR(N_ELEMENTS(a))

    for j=0L,n_elements(a)-1 do begin
        out[j] = PTR_NEW(/ALLOCATE_HEAP)

        n_H = (fn(a[j], "HI_Density") + fn(a[j], "HII_Density")) * density_unit
        T = fn(a[j], "Temperature")
        *out[j] = dust_absorption_coefficient(lambda, T, n_H)
        ;T_sub = 1200; K (estimated sublimation temperature)
        ;sigma_per_H = (1.76459d * lambda^2 - 3.17236d * lambda + 1.69854d) * 1d-21
        ;; now calculate how much dust would survive at this temperature
        ;; using the (vaguely) empirical relation
        ;*out[j] = sigma_per_H * n_H * (1. - 0.875 * (1. / (1. + exp(-(T-T_sub) / 100))))
    endfor
    return, out
end