;+
; NAME:
;  bidr2
; PURPOSE: (one line)
;  Compute the bi-directional reflectance (newer Hapke formula).  This
;    is coded from equation 12.55 on page 346 in Hapke's book, "Theory of
;    Reflectance and Emittance Spectroscopy".
; DESCRIPTION:
; CATEGORY:
;  Miscellaneous
; CALLING SEQUENCE:
;  ans = bidr2(w,emu,imu,g,holes,pzero,b0,theta)
; INPUTS:
;  w     - Single scattering albedo.
;  emu   - Cosine of the emission angle.
;  imu   - Cosine of the incidence angle.
;  g     - Phase angle, in radians.
;  holes - Compaction parameter value (1986 formalism).
;  pzero - Value of the single particle phase function.
;  theta - Surface roughness value.  (radians)
; OPTIONAL INPUT PARAMETERS:
;  None.
; KEYWORD PARAMETERS:
;  H93   - Flag, if set, uses the 1993 version of Hapke's approximation to the
;            Chandresekar H function.  The 1993 version is more accurate but
;            considerably slower to compute.
; OUTPUTS:
;  Return value is the bi-directional reflectance.
; COMMON BLOCKS:
;  None.
; SIDE EFFECTS:
;  None.
; RESTRICTIONS:
;  Any input may be a vector.  If more than one is a vector then the
;     lengths must match.  The return will have the same dimensions as
;     the input.
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written by Marc W. Buie, Lowell Observatory, 1997/08/21
;  97/09/18, MWB, added surface roughness parameter
;-
function bidr2,in_w,in_emu,in_imu,in_g,in_holes,in_pzero,in_b0,in_theta,H93=h93

   check=[n_elements(in_w),n_elements(in_emu),n_elements(in_imu), $
          n_elements(in_g),n_elements(in_holes),n_elements(in_pzero), $
          n_elements(in_b0),n_elements(in_theta)]

   h93 = keyword_set(h93)

   tlen = max(check)
   z=where(check ne 1 and check ne tlen,count)
   IF count ne 0 THEN BEGIN
      print,'BIDR2: Error, lengths of inputs must match or be 1.'
      return,0.0
   ENDIF

   ; Promote all inputs to same length
   IF n_elements(in_w)     eq 1 THEN w     = replicate(in_w,    tlen) ELSE w     = in_w
   IF n_elements(in_emu)   eq 1 THEN emu   = replicate(in_emu,  tlen) ELSE emu   = in_emu
   IF n_elements(in_imu)   eq 1 THEN imu   = replicate(in_imu,  tlen) ELSE imu   = in_imu
   IF n_elements(in_g)     eq 1 THEN g     = replicate(in_g,    tlen) ELSE g     = in_g
   IF n_elements(in_holes) eq 1 THEN holes = replicate(in_holes,tlen) ELSE holes = in_holes
   IF n_elements(in_pzero) eq 1 THEN pzero = replicate(in_pzero,tlen) ELSE pzero = in_pzero
   IF n_elements(in_b0)    eq 1 THEN b0    = replicate(in_b0,   tlen) ELSE b0    = in_b0
   IF n_elements(in_theta) eq 1 THEN theta = replicate(in_theta,tlen) ELSE theta = in_theta

   gamma = sqrt(1-w)

   bscat=replicate(0.,tlen)
   m = where( g lt !pi/2 and b0 ne 0.0, count)
   if count ne 0 then begin
      bscat[m] = 1.0/(1.0 + 1.0/holes[m]*tan(g[m]/2.0))
   end
   bscat = b0*bscat

   bidr = fltarr(tlen)

   ; This is the case where the surface roughness term is turned off.
   z=where(theta eq 0.0,count)
   IF count ne 0 THEN BEGIN

      twoemu = emu[z] + emu[z]
      twoimu = imu[z] + imu[z]

      IF h93 THEN BEGIN
         r0 = 2.0/(1.0+gamma[z])-1.0
         lnt = replicate(0.0,count)
         zn0 = where(emu[z] ne 0.0)
         IF zn0[0] ne -1 THEN lnt[zn0]=alog((1.0+emu[z[zn0]])/emu[z[zn0]])
         hobs = 1.0/(1.0 - (1.0-gamma[z])*emu[z]*(r0+(1.0-0.5*r0-r0*emu[z])*lnt))
         lnt = replicate(0.0,count)
         zn0 = where(imu[z] ne 0.0)
         IF zn0[0] ne -1 THEN lnt[zn0]=alog((1.0+imu[z[zn0]])/imu[z[zn0]])
         hsun = 1.0/(1.0 - (1.0-gamma[z])*imu[z]*(r0+(1.0-0.5*r0-r0*imu[z])*lnt))
      ENDIF ELSE BEGIN
         hobs = (1+twoemu) / (1 + gamma[z] * twoemu)
         hsun = (1+twoimu) / (1 + gamma[z] * twoimu)
      ENDELSE

      bidr[z] = w[z]*imu[z]/(imu[z]+emu[z])*((1+bscat[z])*pzero[z]+hobs*hsun-1)/(4*!pi)

   ENDIF

   ; With surface roughness
   z=where(theta ne 0.0,count)
   IF count ne 0 THEN BEGIN

      emue  = emu[z]
      imue  = imu[z]
      emue0 = emu[z]
      imue0 = imu[z]
      sfun  = emu[z]

      ; ancillary values
      tanthe  = tan(theta[z])
      costhe  = cos(theta[z])
      cotthe  = 1.0/tanthe
      cotthe2 = cotthe^2

      i       = acos(imu[z])
      sini    = sin(i)
      e       = acos(emu[z])
      sine    = sin(e)

      cosg    = cos(g[z])
      cosphi  = replicate(1.0,count)
      zz = where(i*e ne 0.0,countzz)
      if countzz ne 0 then cosphi[zz]=(cosg - imu[z[zz]]*emu[z[zz]])/(sini[zz]*sine[zz])
      zz = where(cosphi gt 1.0,countzz)
      if countzz ne 0 then cosphi[zz] = 1.0
      zz = where(cosphi lt -1.0,countzz)
      if countzz ne 0 then cosphi[zz] = -1.0
      phi     = acos(cosphi)
      sinphi2_2=sin(phi/2.0)^2

      gold=1.0e-7
      z0=where(abs(sini) lt gold,count0)
      IF count0 ne 0 THEN sini[z0]=gold
      z0=where(abs(sine) lt gold,count0)
      IF count0 ne 0 THEN sine[z0]=gold

      coti    = imu[z]/sini
      coti2   = coti^2
      cote    = emu[z]/sine
      cote2   = cote^2
   ;az=where(finite(phi) eq 0,countz)
   ;print,countz

   ;print,minmax(-2.0/!pi*cotthe*coti)
      e1i = exp( -2.0/!pi*cotthe*coti   )    ; eqn. 12.45b, p. 344
      e2i = exp( -1.0/!pi*cotthe2*coti2 )    ; eqn. 12.45c, p. 344
      e1e = exp( -2.0/!pi*cotthe*cote   )
      e2e = exp( -1.0/!pi*cotthe2*cote2 )

      chi = 1.0/sqrt(1.0+!pi*tanthe^2)       ; eqn. 12.45a, p. 344
      fg = exp( -2.0 * tan(phi/2.0) )        ; eqn. 12.51, p. 345

      emue0 = chi * ( emu[z] + sine * tanthe * e2e / ( 2.0 - e1e ) )
      imue0 = chi * ( imu[z] + sini * tanthe * e2i / ( 2.0 - e1i ) )

      ; e >= i
      zz = where(e ge i,countee)
      IF countee ne 0 THEN BEGIN
         denom = 2.0 - e1e[zz] - (phi[zz]/!pi)*e1i[zz]
         imue[zz] = chi[zz] * ( imu[z[zz]] + sini[zz] * tanthe * $
                              ( cosphi[zz]*e2e[zz] + sinphi2_2[zz]*e2i[zz] ) / denom )
         emue[zz] = chi[zz] * ( emu[z[zz]] + sine[zz] * tanthe[zz] * $
                              ( e2e[zz] - sinphi2_2[zz]*e2i[zz] ) / denom )
         sfun[zz] = emue[zz]/emue0[zz] * imu[z[zz]]/imue0[zz] * chi[zz] / $
                     ( 1.0 - fg[zz] + fg[zz]*chi[zz]*imu[z[zz]]/imue0[zz] )
      ENDIF

      ; e < i
      zz = where(e lt i,countee)
      IF countee ne 0 THEN BEGIN
         denom = 2.0 - e1i[zz] - (phi[zz]/!pi)*e1e[zz]
         imue[zz] = chi[zz] * ( imu[z[zz]] + sini[zz] * tanthe[zz] * $
                              ( e2i[zz] - sinphi2_2[zz]*e2e[zz] ) / denom )
         emue[zz] = chi[zz] * ( emu[z[zz]] + sine[zz] * tanthe[zz] * $
                              ( cosphi[zz]*e2i[zz] + sinphi2_2[zz]*e2e[zz] ) / denom )
         sfun[zz] = emue[zz]/emue0[zz] * imu[z[zz]]/imue0[zz] * chi[zz] / $
                     ( 1.0 - fg[zz] + fg[zz]*chi[zz]*emu[z[zz]]/emue0[zz] )
      ENDIF

      twoemu = emue + emue
      twoimu = imue + imue

      IF h93 THEN BEGIN
         r0 = 2.0/(1.0+gamma[z])-1.0
         lnt = replicate(0.0,count)
         zn0 = where(emue ne 0.0)
         IF zn0[0] ne -1 THEN lnt[zn0]=alog((1.0+emue[zn0])/emue[zn0])
         hobs = 1.0/(1.0 - (1.0-gamma[z])*emue*(r0+(1.0-0.5*r0-r0*emue)*lnt))
         lnt = replicate(0.0,count)
         zn0 = where(imue ne 0.0)
         IF zn0[0] ne -1 THEN lnt[zn0]=alog((1.0+imue[zn0])/imue[zn0])
         hsun = 1.0/(1.0 - (1.0-gamma[z])*imue*(r0+(1.0-0.5*r0-r0*imue)*lnt))
      ENDIF ELSE BEGIN
         hobs = (1+twoemu) / (1 + gamma[z] * twoemu)
         hsun = (1+twoimu) / (1 + gamma[z] * twoimu)
      ENDELSE

      bidr[z] = w[z]*imue/(imue+emue)*((1+bscat[z])*pzero[z]+hobs*hsun-1)/(4*!pi)*sfun

   ENDIF

   err=check_math()
   IF (err AND '1B'X) ne 0 THEN $
      print,'BIDR: Math error detected, code ',err

   return,bidr

end