;+
; NAME:
;  oclc_epq03_bnu
; PURPOSE: (one line)
;  calculate bnu from rbot, dtheta
; DESCRIPTION:
;  calculate bnu from rbot, dtheta
; CATEGORY:
;  Occultation lightcurve (oclc)
; CALLING SEQUENCE:
;  bnubot = oclc_epq03_bnu(rbot, rtop, thetamid, deltathetamid, ib, ch)
; INPUTS:
;  rbot = planet radius, at shell bottom border
;  rtop = planet radius, at shell top border
;   rbot[i] = r_{i + 1/2}
;   rtop[i] = r_{i - 1/2}
;  thetamid = bending angle at shell midpoint
;  deltathetamid = change in bending angle over shell
;  ib = index of top of inversion region
;  ch = 1 if adding chapman top, 0 otherwise
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  snubot = bending angle at shell borders
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  
; PROCEDURE:
;  For the bins, follow
;  Elliot, Person and Qu 2003, 
;    eq 35, 57 - similarity of integrand for Inu and Bnu
;    eq. 51, 52 - summation form for Snu approx. to Inu
;  For the topmost bin, note:
;
; eq 13 can be written
;
;            1   / inf 
; nu(r) = - ---  |     theta(r')/r' dx
;            pi  / 0
;           
; We are concerned with the integral above the topmost bin.
; Let's call the radius of the topmost bin R = rtop[0] = r_1/2.
; Then the integral starts at 
; X = sqrt(R^2 - r^2)
;
;               1   / inf 
; nu_ch(r) = - ---  |     theta(r')/r' dx
;               pi  / X
;
; Assuming that theta/r is an exponential with scale height H, this
; integral is
;
;               1  theta(R)
; nu_ch(r) = - --- -------- * H * ch(R/H, asin(r/R) )
;               pi    R
; 
; where ch is the Chapman grazing incidence integral
;  R/H = X_p is the radius/scale_height at base of integral
;  asin(r/R) = zenith angle of integration path at R
;
; MODIFICATION HISTORY:
;  Written 2006 Aug 23, Leslie Young, SwRI
;-

function oclc_epq03_bnu, rbot, rtop, thetamid, deltathetamid, ib, ch

    imax = n_elements(rbot)

    bnubot = dblarr(imax)
    
    ; get thetatop0 and h, scale height and value at 
    ; top of first bin
    thetartop0 = (thetamid[0] - deltathetamid[0]/2.)/rtop[0]
    h_theta = - (rbot[0]-rtop[0]) * thetamid[0]/deltathetamid[0]
    h = 1/(1/h_theta + 1/rtop[0])
    Xp = rtop[0]/h  ; first argment to Smith and Smith    

    for i = ib, imax - 1 do begin

        zm = rtop[0:ib-1]/rbot[i] ; minus
        zp = rbot[0:ib-1]/rbot[i] ; plus

        chi = asin(rbot[i]/rtop[0]) ; chi

        bnubot[i] = -(1/!dpi) *  $
          total( (1/(zp - zm)) * $
                 ( ( zp * acosh(zp) - sqrt(zp^2 - 1)  ) - $
                   ( zm * acosh(zm) - sqrt(zm^2 - 1)  ) ) * $
                 deltathetamid[0:ib-1] )

        if ch eq 1 then begin
            bnubot[i] =  bnubot[i] - (1/!dpi) * $
              ( thetartop0 * h * chapfun(Xp, chi,0.) )
        end
    end

    return, bnubot

end