;+
FUNCTION ABSORPTIVITY, Nbar, T, U, L, Molname, $ 
	               SAME_MOL=same_mol, SAME_LEVPOP=same_levpop

; Block for molecule data
COMMON LINELUM_MOLDATA, $
	molfile, H2opr, nlev, levenergy, levwgt, radupper, radlower, $
	radfreq, radTupper, temptable, colupper, collower, coltable, $
	einstein_a, wgt

; Block to pass around input data
COMMON LINELUM_INPUT, nbar_input, kT_input, tau10_input, mach_input, $
	U_input, L_input

; Return to caller on error
On_Error, 2

; Constants
c = 3.0d10
kb = 1.38d-16
h = 6.63d-27
hc = h*c

; Store input data
kT = kb*T
nbar_input = nbar
kT_input = kT
U_input = U
L_input = L

; ----------------------------------------
; STEP 1: Get molecular data
; ----------------------------------------

if (not keyword_set(same_mol)) and (not keyword_set(same_levpop)) then begin

    ; Read file
    readmol, molname, levtable, radtable, temptable, coltable, molwgt, $
	npart = npart

    ; Handle ortho- and para-H_2 issues
    if npart gt 1 then begin
        H2opr = 0.25
	para = 1.0/(1.0+H2opr)
	ortho = H2opr/(1.0+H2opr)
	temptable1 = *temptable[0]
	coltable1 = para*(*coltable[0])+ortho*(*coltable[1])
	ptr_free, temptable
	ptr_free, coltable
	temptable = temptable1
	coltable = coltable1
    endif

    ; Store the molecular information we want
    molfile = molname
    levenergy = hc*reform(levtable[0,*])
    levwgt = reform(levtable[1,*])
    nlev = (size(levtable))[2]
    radupper = fix(reform(radtable[0,*]))
    radlower = fix(reform(radtable[1,*]))
    einstein_a = dblarr(nlev, nlev)
    einstein_a[radupper, radlower] = radtable[2,*]
    radfreq = reform(radtable[3,*])*1d9
    radTupper = reform(radtable[4,*])
    ntemp = n_elements(temptable)
    colupper = fix(reform(coltable[0,*]))
    collower = fix(reform(coltable[1,*]))
    coltable = reform(coltable[2:ntemp+1,*])
    wgt = molwgt

endif else begin

    ; check that we do have stored molecule information if we've been
    ; told to re-use it
    if n_elements(molfile) eq 0 then begin
	message, 'Error: no molecular data stored, please call without setting SAME_MOL'
    endif

endelse

; ----------------------------------------
; STEP 2: Compute level populations and opacities
; ----------------------------------------

; compute level population if we aren't using a stored one
if not keyword_set(same_levpop) then $
	levpop = getlevpop(nbar, kT_input, /no_beta)

; compute luminosity at this density
transfreq=(levenergy[U_input] - levenergy[L_input])/h
einstein_b = einstein_a[U_input, L_input]*c^2/(2.0*h*transfreq^3) * $
             levwgt[U_input]/levwgt[L_input]
kappa = einstein_b * levpop[L_input]*h*transfreq/(4.0*!pi)

return, kappa

end