pro get_molecule_data, U_input, L_input, LAMBDA=lambda, ENERGY=energy, EINSTEIN_A=ein_A21, $
    EINSTEIN_B12=ein_B12, EINSTEIN_B21=ein_B21, COLLISION=col_data, WEIGHT=molwgt, $
    ABUNDANCE=abundance
    dir = 'TOOLS/projection_with_absorption/krumholz/'
    COMMON LINELUM_MOLDATA, $
        molfile, H2opr, nlev, levenergy, levwgt, radupper, radlower, $
        radfreq, radTupper, temptable, colupper, collower, coltable, $
        einstein_a, wgt

    if 1 then begin
        molname = dir + 'c18o.dat'
        abundance = 1.4d-7
        L_input = 1
        U_input = 2
    endif else begin
        molname = dir + 'n2hplus.dat'
        abundance = 5d-10
        L_input = 0
        U_input = 1
    endelse

    if N_ELEMENTS(molfile) eq 0 then molfile = ''

    if (molfile ne molname) then begin
        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 = constant('h*c')*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

    h = constant('h') & c = constant('c') & pi = constant('pi')
    energy = levenergy[U_input] - levenergy[L_input]
    nu = energy/h
    lambda = c / nu
    ein_A21 = einstein_a[U_input, L_input]
    ein_B21 = ein_A21*c^2/(2.0*h*nu^3)
    ein_B12 = ein_B21*levwgt[U_input]/levwgt[L_input]
    i = (where((colupper eq U_input) and (collower eq L_input)))[0]
    col_data = [[coltable[*,i]], [temptable]]
    molwgt = wgt
end

function build_levpop_table, T, n_H, levels
    kB = constant('kB')
    table = DBLARR(N_ELEMENTS(levels), N_ELEMENTS(T), N_ELEMENTS(n_H))
    for iT = 0L,N_ELEMENTS(T)-1 do begin
        for iH = 0L,N_ELEMENTS(n_H)-1 do begin
            levpop = getlevpop(n_H[iH], kB * T[iT], /NO_BETA)
            table[0,iT,iH] = levpop[levels]
        endfor
    endfor
    return, TRANSPOSE(table, [1,2,0])
end

function interpolate_power_2d, V, VX, VY, X, Y
    idx=min(where(x gt vx))
    idy=min(where(y gt vy))
    if idx lt 0 then idx = N_ELEMENTS(VX)-1
    if idy lt 0 then idy = N_ELEMENTS(VY)-1
    tx=alog(x/vx[idx]) / alog(vx[idx+1]/vx[idx])
    ty=alog(x/vx[idx]) / alog(vx[idx+1]/vx[idx])
    v1 = v[idx,idy  ]^(1-tx) * v[idx+1,idy  ]^tx
    v2 = v[idx,idy+1]^(1-tx) * v[idx+1,idy+1]^tx
    return, v1^(1-ty) * v2^ty
end

; construct a projection in velocity space for the given line
; five radiative effects are considered -
; spontaneous emission, stimulated emission, line absorption,
; dust emission, and dust absorption
; procedure:
; 1. extract data along line of sight
; 2. compute species abundances n1 and n2 for emitting compound
; 3. compute dust cross section and temperature
; 4. compute line emission and absorption cross-sections
; 5. compute line and dust profiles
; 6. project
; set INTEGRATE to get units of ergs/s/cm^2/ster rather than ergs/s/cm^2/ster/Hz
function construct_ppv_projection, PROJECTION=proj, USE=use, PT=point, INTEGRATE=integrate, PLOT_RESULTS=do_plot, TABLE=table
    slice = proj.slice
    xyi = slice.xyi & xi = slice.xi & yi = slice.yi & zi = slice.zi
    grids = *proj.grids & ii_l = *proj.ii_l & ii_r = *proj.ii_r & ii_d = *proj.ii_d & delta_distance = *proj.delta_distance
    border = proj.border
    w = proj.slice.image_size

    get_molecule_data, U_input, L_input, LAMBDA=lambda, ENERGY=energy, EINSTEIN_A=einstein_A, $
        EINSTEIN_B21=einstein_B21, EINSTEIN_B12=einstein_B12, COLLISION=col_data, WEIGHT=molwgt, $
        ABUNDANCE=X_abundance
    table_T   = 10.d^(INDGEN(11)/2.)
    table_n_H = 10.d^(INDGEN(8))
    if N_ELEMENTS(table) eq 0 then $
        table = build_levpop_table(table_T, table_n_H, [L_input, U_input])

    velocity_field = STRMID('xyz',zi,1) + '-velocity'
    get_projection_data, PROJECTION=proj, FIELDS=["Temperature", "Density", velocity_field], /NOPARTS
    grid_data = *proj.grid_data
    
    ; 1. extract data along line of sight
    x = point[0] & y = point[1]
    print, x, y
    xy = CMREPLICATE([x,y], N_ELEMENTS(grids))
    i = where((TOTAL(xy ge grids.left_edge[xyi,*], 1) eq 2) and $
              (TOTAL(xy le grids.right_edge[xyi,*], 1) eq 2))
    if i[0] lt 0 then STOP
    grids = grids[i] & ii_l = ii_l[i,zi] & ii_r = ii_r[i,zi] & ii_d = ii_d[i,zi]
    dxy = (delta_distance[xyi,*])[*,i] & xy = xy[*,0:N_ELEMENTS(i)-1]
    delta_distance = delta_distance[zi,i] & grid_data = grid_data[i]

    l = ((xy - grids.left_edge[xyi,*]) / dxy) + border

    count_total = LONG([0L, TOTAL(ii_d, /CUMULATIVE)]) ; cumulative sum with 0 tacked on the front
    n_depths = count_total[N_ELEMENTS(grids)]
    int_depth_values = LON64ARR(n_depths) ; pre-allocate arrays
    grid_index = INTARR(n_depths)
    grid_values = LONARR(n_depths)
    depth_index = LINDGEN(n_depths)

    dims = grid_int_coords(INFO=grids)
    scale = max(grids.level) - [grids.level]
    dims = ISHFT(REFORM(dims[zi,*])+ii_l, scale)
    int_width = ISHFT((grids[(where([grids.level] eq 0))[0]].dim[zi]), max(grids.level))
    for i=0L,N_ELEMENTS(grids)-1 do begin
        X = l[0,i] & Y = l[1,i] & Z = DINDGEN(ii_d[i])+border
        for k=0L,grid_data[i].n_dfields-1 do $
            grid_data[i].data[k] = PTR_NEW(REFORM(INTERPOLATE(*grid_data[i].data[k], X, Y, Z, /GRID)))
        int_depth_values[count_total[i]] = dims[i] + ISHFT(ii_l[i]+L64INDGEN(ii_d[i]), scale[i])
        grid_index[count_total[i]] = INDGEN(ii_d[i])
        grid_values[count_total[i]] = LONARR(ii_d[i])+i
    endfor
    unique_depths = UNIQ(int_depth_values, SORT(int_depth_values)) ; calculate index list for unique depths
    depth_values = DOUBLE(int_depth_values) / int_width
    n_depths = N_ELEMENTS(unique_depths)
    for i=0L,n_depths-1 do $
        depth_index[where(depth_values eq depth_values[unique_depths[i]])] = i

    flat_data = DBLARR(n_depths, grid_data[0].n_dfields)
    flat_depth = depth_values[unique_depths]
    flat_thickness = DBLARR(n_depths)
    order = sort(grids.level)
    for i=0L, N_ELEMENTS(grids)-1 do begin
        j = order[i]
        depth_indices = depth_index[where(grid_values eq j)]
        for k=0L, grid_data[j].n_dfields-1 do $
            flat_data[depth_indices, k] = (*grid_data[j].data[k])[0:ii_d[j]-1]
        flat_thickness[depth_indices] = delta_distance[j] * proj.distance_unit
    endfor

    names = strtrim(grid_data[0].data_fields,2)
    Density  = (where(names eq "Density"  ))[0]
    Temperature = (where(names eq "Temperature" ))[0]
    Velocity    = (where(names eq velocity_field))[0]
    n_H = flat_data[*,Density] * proj.density_unit
    T   = flat_data[*,Temperature]
    vel = flat_data[*,Velocity] * get_unit('velocity', /FORCE_USE_UNITS) * 1e5

    vel_d = [0., .25*abs(vel[0:n_depths-3] - vel[2:n_depths-1]), 0.]

    ; 2. compute species abundances of n1 and n2 for emitting compound
    n_X1 = DBLARR(n_depths)
    n_X2 = DBLARR(n_depths)
    kB = constant('kB')
    n_X = n_H * X_abundance
    for i=0L,n_depths-1 do begin
        ;levpop = getlevpop(n_H[i], kB * T[i], /NO_BETA)
        n_X1[i] = n_X[i] * interpolate_power_2d(table[*,*,0], table_T, table_n_H, T[i], n_H[i]) ;levpop[L_input]
        n_X2[i] = n_X[i] * interpolate_power_2d(table[*,*,1], table_T, table_n_H, T[i], n_H[i]) ;levpop[U_input]
    endfor
    
    ; 3. compute dust absorption coefficient and temperature
    dust_alpha_nu = dust_absorption_coefficient(lambda, T, n_H)
    dust_j_nu     = dust_emissivity(dust_alpha_nu, lambda, T)

    ; 4. compute line emission and absorption coefficients
    ; B ~ cm^2 erg^-1 Hz ~ cross section / energy / spectral distribution
    ; h nu n B phi ~ cm-1 
    line_alpha_nu_per_phi = energy * (n_X1 * einstein_B12 - n_X2 * einstein_B21)
    line_j_nu_per_phi     = energy * einstein_A * n_X2 / (4*constant('pi'))
    ; while we're at it, compute collision rate
    collision_rate = INTERPOL(col_data[*,0], col_data[*,1], T, /LSQUADRATIC)

    ; 5. compute line and dust profiles
    v_min = -20 * 1e5 ; 20 km/s
    v_max = 20 * 1e5
    v_steps = 401
    v_step = (v_max - v_min) / (v_steps - 1)
    v = v_min + v_step*DINDGEN(v_steps)

    line_spectral_profile = DBLARR(v_steps, n_depths)
    for i=0L,n_depths-1 do begin
        line_spectral_profile[0,i] = line_profile(v, LAMBDA=lambda, TEMP=T[i], MOLAR_MASS=molwgt, $
            EINSTEIN_A=einstein_A, COLLISION_RATE=collision_rate[i], VELOCITY=vel[i], $
            VEL_DISP=vel_d[i])
    endfor

    ; now compute total absorption and emission coefficients
    alpha_nu = use[0]*CMREPLICATE(line_alpha_nu_per_phi, v_steps) * TRANSPOSE(line_spectral_profile) + use[2]*CMREPLICATE(dust_alpha_nu, v_steps)
    j_nu     = use[1]*CMREPLICATE(line_j_nu_per_phi    , v_steps) * TRANSPOSE(line_spectral_profile) + use[3]*CMREPLICATE(dust_j_nu, v_steps)

    ; compute optical depths to the front of each volume element
    ds = CMREPLICATE(flat_thickness, v_steps)
    dtau = ds*alpha_nu
    tau = [DBLARR(1,v_steps),TOTAL(dtau, 1, /CUMULATIVE)]
    ; compute emergent intensity at the front of each volume element
    ix = where((dtau gt 0.01) and (alpha_nu lt 0.01)) ; substitute Taylor polynomial where calculation is poorly conditioned
    ;efficiency = ds-alpha_nu*ds^2/2.+alpha_nu^2*ds^3/6.-alpha_nu^3*ds^4/24.
    efficiency = ds*(-.5D*dtau*(.33333D*dtau*(.25D*dtau-1.D)+1.D)+1.D)
    if ix[0] ge 0 then efficiency[ix] = (1.D - exp(-dtau[ix])) / alpha_nu[ix]
    emergent_dI_nu = efficiency * j_nu
    ; diminish emergent intensities by line-of-sight absorption
    dI_nu = emergent_dI_nu * exp(-tau[0:n_depths-1,*])
    ; sum and serve with garnish
    I_nu = TOTAL(dI_nu, 1)

    if KEYWORD_SET(integrate) then I_nu *= v_step / lambda
    if KEYWORD_SET(do_plot) then plot, v*1e-5, I_nu
    return, I_nu
end