NAME:
  bidr4
 PURPOSE: (one line)
  Compute the bi-directional reflectance (Hapke 2012+enhancements).
 DESCRIPTION:
  This function is based on the rough reflectance treatment
    from equation 12.55 on page 323 in Hapke's book,
    "Theory of Reflectance and Emittance Spectroscopy", second edition.
    All aspects of the Hapke formalism that is built for rough surfaces
    is used.  Also included is an additional modification and the
    relationship between K and h_s taken from Helfenstein & Shepard,
    Icarus, v215, p83, (2011).
  This function is intended to be as close as possible to bidr1 and bidr2
     but the change in the theory does require additional inputs.
 CATEGORY:
  Miscellaneous
 CALLING SEQUENCE:
  ans = bidr4(ssalb,emu,imu,phang,hs,hc,ppp,b0s,b0c,theta)
 INPUTS:
  all input arguments are protected against modification by this program
  ssalb - Single scattering albedo.
  emu   - Cosine of the emission angle.
  imu   - Cosine of the incidence angle.
  phang - Phase angle, in radians.
  hs    - Shadow hiding (SHOE) compaction parameter value (2012 formalism).
  hc    - Coherent backscatter (CBOE) compaction parameter (2012 formalism).
  ppp   - Parameters of the single particle phase function
          (default "function" is a constant of value ppp
          There are four legal input forms for ppp
            1. a scalar
            2. an array of dimensionality same as ssalb, emu, imu, etc.
            3. an array of dimensionality Pparms
            4. an array of dimensionality (n_elements(ssalb),Pparms)
  b0s   - SHOE backscatter amplitude.
  b0c   - CBOE backscatter amplitude.
  theta - Surface roughness value.  (radians)
 OPTIONAL INPUT PARAMETERS:
  None.
 KEYWORD PARAMETERS:
  Pfn   - Specify procedure to use for P(phang) 
            The routine must be a procedure taking arguments phang,a,F,/radians
              phang - phase angle in radians (with keyword /radians set)
              a     - an array of Pparms parameters
              F     - phase function evaluated at phase angles phang
            Use "fn_hg3.pro" as a model.
  Pparms- Specify number of parameters to be passed to P(phang).
             (default=1)
  pedantic - Flag, if set, returns NaN for ssalb > 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, 2023/08/01
    and was cloned from bidr2