;+
; 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.

            if (wmid gt 2809.030 and wmid lt 2818.894) then begin
                idebug=1
                print, 'DEBUG'
            endif else begin
                idebug=0
            endelse
            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

;            if jb eq 0 then $
;              print, wav[lw1], wav[lw2], form='(2F11.3)'
; FORTRAN
; for jb=0
;    2809.030   2818.894
;    2819.108   2828.851
; for jb=1
;   3900.946   3910.565
;   3911.039   3920.634
; ...
;   4261.108   4270.938
;   4270.985   4278.759



            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)

            if idebug eq 1 then begin
                iz = 1
                print, 'a[2,*]', a[iz,*],form='(A8,3E11.3)'
                print, 'a[1,*]', a[iz-1,*],form='(A8,3E11.3)'
                print, 'da[2,*]', a[iz,*]-a[iz-1,*],form='(A8,3E11.3)'
                print, 'ncol[2]', ncol[iz],form='(A8,3E11.3)'
                print, 'ncol[1]', ncol[iz-1],form='(A8,3E11.3)'
                print, 'dncol[1]', ncol[iz]-ncol[iz-1],form='(A8,3E11.3)'
                print, 'u*w*da*dncol', $
                  xu*wu*(a[iz,*]-a[iz-1,*])/(ncol[iz]-ncol[iz-1]),$
                  form='(A8,3E11.3)'
            endif

            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
                if (idebug eq 1 and iz eq 1) then begin
                    print, solflux,delw,sum,h[iz], form='(4E11.3)'
                    form='(7E11.3)' 
                endif
            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

            if (idebug eq 1) then begin
                iz=1
                daddend = x2[iz]*(dwr+fbr*x1[iz])*h[iz]*n[iz]
                print, jb+1, wav[lw1],wav[lw2], daddend,$
                  format='(I3,2F11.3,E11.3)'
                print, x2[iz],dwr,fbr,x1[iz],h[iz],n[iz], daddend, $
                  form='(7E11.3)' 
            endif
;          -----------------------------------------
;              nu-4 source function
;          -----------------------------------------
            th = ht + rw*x2*h
              
            wlast = wnext
            
        endwhile
        
    endfor
    
    return, ht

end