;+
; 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)
;  col[ntau,5] - vertical column densities of N2, CH4, C2H2, C2H6, HCN (cm^-2)
;  den[ntau,5] - number densities of N2, CH4, C2H2, C2H6, HCN (cm^-3)
;  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
;-
function rt_uvheat, rh, col, den, wave, dwave, flux, crs, ht4

physconstants

ntau = ( size(col) )[1]
nwav = n_elements(wave)

wave0 = [1265.d0, 1870.d0, 3000.d0, 2000.d0]

ht4 = dblarr(ntau,4)

; 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, 3 do begin
  eff = ( (1.d - wave/wave0[js]) > 0)
  fac = 0.5d * ergs * eff * dwave
  for jt = 0, ntau-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.)) )
  endfor
endfor
ht = total(ht4,2)

return, ht

end

pro rt_uvheatTEST

; example for p=10 ubar, fCH4=0.1
physconstants
psurf=10
tsurf=100.
nsurf=psurf/(!phys.k * tsurf)
h = 50 ; scale height in km
ntau = 50
dz = 20
z = (findgen(ntau)+0.5)*dz
den = fltarr(ntau,5)
n = nsurf*exp(-z/h)
den[*,0] = n*0.9
den[*,1] = n*0.1
col = den * (h*1e5) ; multiply by H in cm

rt_solread, wave, dwave, flux, crs, nsol, cycle=1

rh = 30. ; Pluto-ish

ht = rt_uvheat(rh, col, den, wave, dwave, flux, crs, ht4)

plot, z, ht
oplot, z, ht4[*,0]>1e-30, /li  ; clip at 1e-30 to avoid "floating underflow"
oplot, z, ht4[*,1]>1e-30, li=2

; some diagnostic printing

ergs = (1.d8/wave)*(!phys.h * !phys.c) * flux/(rh^2) 
n2 = where(wave lt 1000.)
ch4 = where(wave gt 1000. and wave lt 1500.)

print, 'total N2 heating (erg/cm^2/s)= ', total(ht4[*,0]) * dz*1e5
print, 'total solar flux lam<1000', total(ergs[n2]*dwave[n2])/4.
print, 'total CH4 heating (erg/cm^2/s)= ', total(ht4[*,1]) * dz*1e5
print, 'total solar flux 1000<lam<1500', total(ergs[ch4]*dwave[ch4])/4.



end