NAME:
  bidr2
 PURPOSE: (one line)
  Compute the bi-directional reflectance (newer Hapke formula).
 DESCRIPTION:
    This function is coded from equation 12.55 on page 346 in Hapke's book,
    "Theory of Reflectance and Emittance Spectroscopy".
 CATEGORY:
  Miscellaneous
 CALLING SEQUENCE:
  ans = bidr2(w,emu,imu,g,holes,p,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).
  p     - Parameters of the single particle phase function
          (default "function" is a constant of value p).
          There are four legal input forms for p:
            1. a scalar
            2. an array of dimensionality same as w, emu, imu, etc.
            3. an array of dimensionality Pparms
            4. an array of dimensionality (n_elements(w),Pparms)
  b0    - Backscatter value.
  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.
  Pfn   - Specify function to use for P(g) instead of the default constant.
            The function must be a procedure taking arguments g,a,F,/radians
              g  phase angle in radians (with keyword /radians set)
              a  an array of Pparms parameters
              F  phase function evaluated at phase angles g
            Use "fn_hg3.pro" as a model.
  Pparms- Specify number of parameters to be passed to P(g).
             Really only needed if the number of parameters isn't one.
  pedantic - Flag, if set, returns NaN for w > 1.
 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
  1997/09/18, MWB, added surface roughness parameter
  2001/10/29, WMG, added "Pfn", "Pparms", and "pedantic" parameters
  2003/04/03, MWB, changes loop counters to long
  2004/06/24, MWB, fixed bug in imue equation, error in calculated
                       value happened to be small but was a real
                       computational error.  Thanks to Frank Seelos,
                       EPRSL, Washington University in St. Louis