;+
; NAME:
;  rt_uvheat
; PURPOSE: (one line)
;  calculate UV heating rate
; DESCRIPTION:
;  calculate global average of UV heating rate given data from solread
; CATEGORY:
;  RT
; CALLING SEQUENCE:
;  ht = rt_uvheat(rh,col,den,wave,flux,crs)
; INPUTS:
;  rh - heliocentric distance (AU)
;  atm structure with
;   ncol[nz,5] - vertical column densities of N2, CH4, C2H2, C2H6, HCN (cm^-2)
;   n[nz,5] - number densities of N2, CH4, C2H2, C2H6, HCN (cm^-3)
;  sol structure with
;   wave[nsol] - wavelength in A
;   dwave[nsol] - wavelength interval in A
;   flux[nsol] - flux in photon/cm^2/s/A
;   crs[nsol,4] - cross section of N2, CH4, C2H2, C2H6 (cm^2)  
; OPTIONAL INPUT PARAMETERS:
; KEYWORD INPUT PARAMETERS:
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  ht[ntau] - heating rate in erg/cm^3/s
;  ht4[ntau, 4] - heating rate in erg/cm^2/s for N2, CH4, C2H2, C2H6
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2004 Mar 20, Leslie Young SwRI
;  based on uvheat.f, written by Roger Yelle, 1990
;  with more array arithmetic
;  2004 May 9, LAY SwRI - take uv structure as input
;-
function rt_uvheat, rh, atm, sol, h
    
    physconstants
    
    nz = atm.nz
    col = atm.ncol
    den = atm.n
    
    nwav = sol.nsol
    ncrs = sol.ncrs
    wave = sol.wave
    dwave = sol.dwave
    flux = sol.flux * dwave
    crs = sol.crs
    
    
    wave0 = [1265.d0, 1870.d0, 3000.d0, 2000.d0]

    h = dblarr(nz,ncrs, nwav) ; debugging purposes    
    ht4 = dblarr(nz,ncrs)
    
; convert to flux in erg/cm^2/s/dlam at distance rh
    ergs = (1.d8/wave)*(!phys.h * !phys.c) * flux/(rh^2)

; the call to expint is to perform the global average:
; expint(2,tau) = integral( exp(-tau/mu),mu=0 to 1)
;
; the efficiency factor accounts for how much energy
; goes into local heating and how much into ionospheric chemistry.
; This rough formula is directly from RVY's 1990 version of uvheat.f
    
    for js = 0, ncrs-1 do begin
        eff = ( (1.d - wave/wave0[js]) > 0)
        fac = 0.5d * ergs * eff
        for jt = 0, nz-1 do begin
            opp = den[jt, js] * crs[*, js]
            tau = col[jt,js] * crs[*,js]
            ht4[jt,js] = total( fac * opp * expint(2, double(tau<100.)) )
            h[jt,js,*] = fac * opp * expint(2, double(tau<100.)) 
;            if js eq 2 and jt eq nz-1  then begin
;                jw = nwav-1
;                print, 'uvheat: js = 1, lowest atm, last wavelength'
;                print, 'wave = ', wave[jw]
;                print, 'den = ', den[jt,js]
;                print, 'crs = ', crs[jw,js]
;                print, 'opp = ', opp[jw]
;                print, 'tau = ', tau[jw]
;                print, 'en2(tau)= ', expint(2, double(tau[jw]<100.))
;                print, 'fac = ', fac[jw]
;                print, 'h = ', h[jt,js,jw]
;            end
        endfor
        
    endfor
    ht = total(ht4,2)

    return, ht

end