docs: def C_Reader_Harris_Gallagher(D, Do, rho, mu, m, taps='corner'):
      r'''Calculates the coefficient of discharge of the orifice based on the 
      geometry of the plate, measured pressures of the orifice, mass flow rate
      through the orifice, and the density and viscosity of the fluid.
      
      .. math::
          C = 0.5961 + 0.0261\beta^2 - 0.216\beta^8 + 0.000521\left(\frac{
          10^6\beta}{Re_D}\right)^{0.7}\\
          + (0.0188 + 0.0063A)\beta^{3.5} \left(\frac{10^6}{Re_D}\right)^{0.3} \\
          +(0.043 + 0.080\exp(-10L_1) -0.123\exp(-7L_1))(1-0.11A)\frac{\beta^4}
          {1-\beta^4} \\
          -  0.031(M_2' - 0.8M_2'^{1.1})\beta^{1.3}
          
      .. math::
          M_2' = \frac{2L_2'}{1-\beta}
          
      .. math::
          A = \left(\frac{19000\beta}{Re_{D}}\right)^{0.8}
          
      .. math::
          Re_D = \frac{\rho v D}{\mu}
          
          
      If D < 71.12 mm (2.8 in.) (Note this is a continuous addition; there is no
      discontinuity):
          
      .. math::
          C += 0.11(0.75-\beta)\left(2.8-\frac{D}{0.0254}\right)
          
      If the orifice has corner taps:
          
      .. math::
          L_1 = L_2' = 0
          
      If the orifice has D and D/2 taps:
          
      .. math::
          L_1 = 1
          
      .. math::
          L_2' = 0.47
      
      If the orifice has Flange taps:
      
      .. math::
          L_1 = L_2' = \frac{0.0254}{D}
          
      Parameters
      ----------
      D : float
          Upstream internal pipe diameter, [m]
      Do : float
          Diameter of orifice at flow conditions, [m]
      rho : float
          Density of fluid at `P1`, [kg/m^3]
      mu : float
          Viscosity of fluid at `P1`, [Pa*s]
      m : float
          Mass flow rate of fluid through the orifice, [kg/s]
      taps : str
          The orientation of the taps; one of 'corner', 'flange', 'D', or 'D/2',
          [-]
          
      Returns
      -------
      C : float
          Coefficient of discharge of the orifice, [-]
  
      Notes
      -----
      The following limits apply to the orifice plate standard [1]_:
          
      The measured pressure difference for the orifice plate should be under
      250 kPa.
      
      There are roughness limits as well; the roughness should be under 6 
      micrometers, although there are many more conditions to that given in [1]_.
          
      For orifice plates with D and D/2 or corner pressure taps:
          
      * Orifice bore diameter muse be larger than 12.5 mm (0.5 inches)
      * Pipe diameter between 50 mm and 1 m (2 to 40 inches)
      * Beta between 0.1 and 0.75 inclusive
      * Reynolds number larger than 5000 (for :math:`0.10 \le \beta \le 0.56`)
        or for :math:`\beta \ge 0.56, Re_D \ge 16000\beta^2`
        
      For orifice plates with flange pressure taps:
  
      * Orifice bore diameter muse be larger than 12.5 mm (0.5 inches)
      * Pipe diameter between 50 mm and 1 m (2 to 40 inches)
      * Beta between 0.1 and 0.75 inclusive
      * Reynolds number larger than 5000 and also larger than
        :math:`170000\beta^2 D`.
        
      This is also presented in Crane's TP410 (2009) publication, whereas the
      1999 and 1982 editions showed only a graph for discharge coefficients.
      
      Examples
      --------
      >>> C_Reader_Harris_Gallagher(D=0.07391, Do=0.0222, rho=1.165, mu=1.85E-5, 
      ... m=0.12, taps='flange')
      0.5990326277163659
      
      References
      ----------
      .. [1] American Society of Mechanical Engineers. Mfc-3M-2004 Measurement 
         Of Fluid Flow In Pipes Using Orifice, Nozzle, And Venturi. ASME, 2001.
      .. [2] ISO 5167-2:2003 - Measurement of Fluid Flow by Means of Pressure 
         Differential Devices Inserted in Circular Cross-Section Conduits Running
         Full -- Part 2: Orifice Plates.  
      .. [3] Reader-Harris, M. J., "The Equation for the Expansibility Factor for
         Orifice Plates," Proceedings of FLOMEKO 1998, Lund, Sweden, 1998: 
         209-214.
      .. [4] Reader-Harris, Michael. Orifice Plates and Venturi Tubes. Springer, 
         2015.
      '''
      A_pipe = 0.25*pi*D*D
      v = m/(A_pipe*rho)
      Re_D = rho*v*D/mu
      Re_D_inv = 1.0/Re_D
      
      beta = Do/D
      if taps == 'corner':
          L1, L2_prime = 0.0, 0.0
      elif taps == 'D' or taps == 'D/2':
          L1 = 1.0
          L2_prime = 0.47
      elif taps == 'flange':
          L1 = L2_prime = 0.0254/D
      else:
          raise Exception('Unsupported tap location')
          
      beta2 = beta*beta
      beta4 = beta2*beta2
      beta8 = beta4*beta4
      
      A = (19000.0*beta*Re_D_inv)**0.8
      M2_prime = 2.0*L2_prime/(1.0 - beta)
      
      delta_C_upstream = ((0.043 + 0.080*exp(-1E1*L1) - 0.123*exp(-7.0*L1))
              *(1.0 - 0.11*A)*beta4/(1.0 - beta4))
      
      # The max part is not in the ISO standard
      t1 = log10(3700.*Re_D_inv)
      if t1 < 0.0:
          t1 = 0.0
      delta_C_downstream = (-0.031*(M2_prime - 0.8*M2_prime**1.1)*beta**1.3
                            *(1.0 + 8.0*t1))
      
      # C_inf is discharge coefficient with corner taps for infinite Re
      # Cs, slope term, provides increase in discharge coefficient for lower
      # Reynolds numbers.
      x1 = (1E6*Re_D_inv)**0.3
      x2 = 22.7 - 0.0047*Re_D
      t2 = x1 if x1 > x2 else x2
      # max term is not in the ISO standard
      C_inf_C_s = (0.5961 + 0.0261*beta2 - 0.216*beta8 
                   + 0.000521*(1E6*beta*Re_D_inv)**0.7
                   + (0.0188 + 0.0063*A)*beta**3.5*(
                   t2))
      
      C = (C_inf_C_s + delta_C_upstream + delta_C_downstream)
      if D < 0.07112:
          # Limit is 2.8 inches, .1 inches smaller than the internal diameter of
          # a sched. 80 pipe.
          # Suggested to be required not becausue of any effect of small
          # diameters themselves, but because of edge radius differences.
          # max term is given in [4]_ Reader-Harris, Michael book
          # There is a check for t3 being negative and setting it to zero if so
          # in some sources but that only occurs when t3 is exactly the limit
          # (0.07112) so it is not needed
          t3 = (2.8 - D/0.0254)
          delta_C_diameter = 0.011*(0.75 - beta)*t3
          C += delta_C_diameter
      
      return C