; nu in s^-1, T in Kelvin, B in erg/s/cm^2/sr/Hz
function B_planck, nu, T
    return, constant('2*h/c^2') * nu^3 / (exp(constant('h/kB') * nu / T) - 1.D)
end

function line_profile, v, LAMBDA=lambda, TEMPERATURE=T, MOLAR_MASS=m, $
    EINSTEIN_A=einstein_A, COLLISION_RATE=collision_rate, VELOCITY=vel, VEL_DISP=vel_disp
    c = constant('c')
    pi = constant('pi')
    kB = constant('kB')
    delta_nu = (v-vel)/lambda

    delta_nu_D = ((2*kB*T/(m*1.6605d-24) + vel_disp^2)^.5)/lambda
    Gamma = einstein_A + 2*collision_rate
    u = delta_nu / delta_nu_D
    a = Gamma / (4 * pi * delta_nu_D)
    phi = VOIGT(a, u) / (delta_nu_D * pi^.5)
    return, phi
end

; velocities in km/s
function line_emissivity, LINE_INDEX=line_index, LINE=line_name, AVAILABLE_LINES=return_lines, NEEDED_FIELDS=return_fields, DATA=a
    n_available_lines = 6
    ; line types:
    ; 0 - log temperature quadratic
    ; 1 - collisionally de-excited
    ; 2 - LTE 
    available_lines = replicate({line,name:'',lambda:0.D,needed_fields:PTR_NEW(['']),logT_quadratic:[0.D,0.D,0.D],extra_data:PTR_NEW(0)}, n_available_lines)
    available_lines[0] = {line,'H-alpha',656.3D,PTR_NEW(["HII_Density", "Electron_Density", "Temperature"]),[-21.815977, -0.37414571 , -0.071248393], PTR_NEW(0)}
    available_lines[1] = {line,'H-beta' ,486.1D,PTR_NEW(["HII_Density", "Electron_Density", "Temperature"]),[-23.021928, -0.074316281, -0.099478146], PTR_NEW(0)}
    available_lines[2] = {line,'H-gamma',434.1D,PTR_NEW(["HII_Density", "Electron_Density", "Temperature"]),[-23.581513,  0.013914791, -0.10719574 ], PTR_NEW(0)}
    available_lines[3] = {line,'H-delta',410.2D,PTR_NEW(["HII_Density", "Electron_Density", "Temperature"]),[-23.927480,  0.048607817, -0.11034571 ], PTR_NEW(0)}
    available_lines[4] = {line,'[OIII] 500.7nm',500.7D,PTR_NEW(["HI_Density", "HII_Density", "Electron_Density", "Temperature"]),[0.,0.,0.], $
        PTR_NEW({collision,abundance:3e-4,Omega:2.5,ein_A:2.8e-2,weight:8./9.})}
    available_lines[5] = {line,'[NII] 658.4nm',658.4D,PTR_NEW(["HI_Density", "HII_Density", "Electron_Density", "Temperature"]),[0.,0.,0.], $
        PTR_NEW({collision,abundance:5e-5,Omega:2.99,ein_A:3e-3,weight:2./9.})}
    if KEYWORD_SET(return_lines) then return, available_lines
    if KEYWORD_SET(line_name) then begin
        index = (where(available_lines.name eq line_name))[0]
    endif else begin
        index = line_index
    endelse
    if index lt 0 then return, -1
    p = available_lines[index].logT_quadratic
    ; squaring 5.98802e+23 (float) yields Inf, so convert to double
    density_unit = DOUBLE(get_unit('Electron_Density', /FORCE_USE_UNITS))

    fields = *available_lines[index].needed_fields
    if KEYWORD_SET(return_fields) then return, fields

    ; find the indices corresponding to fields
    data_names = strtrim(a[0].data_fields,2)
    field_indices = INTARR(n_elements(fields))
    for i=0,n_elements(fields)-1 do field_indices[i] = (where(data_names eq fields[i]))[0]

    out = PTRARR(N_ELEMENTS(a), /ALLOCATE_HEAP)

    for j=0L,n_elements(a)-1 do begin
        if (where([0,1,2,3] eq index))[0] ge 0 then begin
            logT  = alog10(*a[j].data[field_indices[2]])
            logT2 = logT^2
            
            *out[j] = (*a[j].data[field_indices[0]]) * $
                     (*a[j].data[field_indices[1]]) * $
                     10.^(p[0] + p[1]*logT + p[2]*logT2) * $
                     density_unit^2
        endif else if (where([4,5] eq index))[0] ge 0 then begin
            c = constant('c')
            h = constant('h')
            kB = constant('kB')
            nu = c / available_lines[index].lambda
            E = h * nu
            params = *available_lines[index].extra_data

            n_HII = fn(a[j], "HII_Density") * density_unit
            n_HI = fn(a[j], "HI_Density") * density_unit
            n_e = fn(a[j], "Electron_Density") * density_unit
            T = fn(a[j], "Temperature")
            n_X = (n_HI + n_HII) * params.abundance
            ; charge transfer
            ; H+ + X -> H + X+
            n_X2 = n_X / (1+(1./params.weight)*(n_HI/n_HII))
            g = 8.63e-6 * params.Omega / T^(-.5)
            *out[j] = n_e * n_X2 * g * params.ein_A * exp(-E/(kB*T)) / $
                      (params.ein_A + g * n_e)
        endif
    endfor
    return, out
end