;+
; NAME:
;  rt_cascade
; PURPOSE: (one line)
;  calculate methane cascade heating
; DESCRIPTION:
;
; CATEGORY:
;  RT
; CALLING SEQUENCE:
;  hcsc = rt_cascade(efc, rmass, rh, $
;                 wav0, wdel, xu, wu, nu, atm, ep1, ep2, $
;                 ssolar, th)
; INPUTS:
;
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  cascade heating rate in erg/cm^3/s
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
; 
; 
; MODIFICATION HISTORY:
;  Written 2004 May 8, Leslie Young SwRI
;  Modified from cascade.f, Roger Yelle
;-

function rt_cascade,efc, rmass, rh, $
                   wav0, wdel, xu, wu, nu, atm, ep1, ep2, $
                   ssolar, th

    tmp0 = 296.d0 
    tmax = 7.d0 
    
    nb = ssolar.nb
    wav = ssolar.wav
    
    nz = atm.nz
    x2 = ep2/(1.d + ep2)
    x1 = ep1/(1.d + ep1)
    im = 1 ; index of CH4 in atm structure arrays
    ncol = sqrt(atm.ncol[0:nz-2,im]*atm.ncol[1:nz-1,im])
    ncol=[ncol,ncol[nz-2]^2/ncol[nz-3]]
    t = 0.5*(atm.t[0:nz-2]+atm.t[1:nz-1])
    t = [t,t[nz-2]]
    p = sqrt(atm.p[0:nz-2]*atm.p[1:nz-1])
    p=[p,p[nz-2]^2/p[nz-3]]
    n = atm.n[*,im]

;       -------------------------------------------------------
;                      Loop through near IR bands
;       -------------------------------------------------------
    th = dblarr(nz)
    ha = dblarr(nz)
    ht = dblarr(nz)
    h = dblarr(nz)
    for jb = 0, nb-1 do begin
;    for jb = nb-1, nb-1 do begin  ; DEBUG DEBUG DEBUG
        i0 = ssolar.nl[jb]
        ni = ssolar.nlines[jb]
        i1 = ssolar.nl[jb]+ni-1
        wend = wav[i1]
        wlast = wav[i0]

        while (wlast lt wend) do begin

            wnext = min([wlast+wdel,wend])
            wmid =(wlast+wnext)/2.
            delw = wnext - wlast
            solflux = rt_solirflx(wmid)/rh^2.

            wavi = wav[i0:i1]
            lw1 = min( [ rt_locate(wavi,wlast)+1, ni-1] )+i0
            lw2 = max( [ rt_locate(wavi,wnext),   0 ] )+i0
            nlw = lw2 - lw1 + 1

            a = dblarr(nz,nu)
            str=ssolar.str[lw1:lw2]
            wid=ssolar.wid[lw1:lw2]
            eng = ssolar.eng[lw1:lw2]
            for iu = 0, nu-1 do $
              a[*,iu] = rt_random(xu[iu],p,t,ncol,$
                  str,wid,eng, efc,rmass,wmid,delw)

            for iz = 0, nz-1 do begin
                sum = 0.
                for i = 0, nu-1 do begin
                    if (iz gt 0) then begin
                        sum=sum+xu[i]*wu[i]*(a[iz,i]-a[iz-1,i]) $
                          /(ncol[iz]-ncol[iz-1])
                    endif else begin
                        sum=sum+xu[i]*wu[i]*a[0,i]/ncol[0]
                    endelse
                endfor
                h[iz] = 0.5d*solflux*delw*sum
            endfor
            
            ha = ha + h*n

;          -----------------------------------------
;              Cascade heating rate (ergs/cm3/s)
;          -----------------------------------------
            fb = float(jb)+1.
            dwr = (wmid-(fb+1.)*wav0)/wmid
            rw = wav0/wmid
            fbr = fb*rw
            ht = ht + x2*(dwr+fbr*x1)*h*n

;          -----------------------------------------
;              nu-4 source function
;          -----------------------------------------
            th = ht + rw*x2*h
              
            wlast = wnext
            
        endwhile
        
    endfor
    
    return, ht

end