code: def _Kw(rho, T):
      
      # Check input parameters
      if rho < 0 or rho > 1250 or T < 273.15 or T > 1073.15:
          raise NotImplementedError("Incoming out of bound")
  
      # The internal method of calculation use rho in g/cm³
      d = rho/1000.
  
      # Water molecular weight different
      Mw = 18.015268
  
      gamma = [6.1415e-1, 4.825133e4, -6.770793e4, 1.01021e7]
      pKg = 0
      for i, g in enumerate(gamma):
          pKg += g/T**i
  
      Q = d*exp(-0.864671+8659.19/T-22786.2/T**2*d**(2./3))
      pKw = -12*(log10(1+Q)-Q/(Q+1)*d*(0.642044-56.8534/T-0.375754*d)) + \
          pKg+2*log10(Mw/1000)
      return pKw
docs: Equation for the ionization constant of ordinary water
  
      Parameters
      ----------
      rho : float
          Density, [kg/m³]
      T : float
          Temperature, [K]
  
      Returns
      -------
      pKw : float
          Ionization constant in -log10(kw), [-]
  
      Notes
      ------
      Raise :class:`NotImplementedError` if input isn't in limit:
  
          * 0 ≤ ρ ≤ 1250
          * 273.15 ≤ T ≤ 1073.15
  
      Examples
      --------
      >>> _Kw(1000, 300)
      13.906565
  
      References
      ----------
      IAPWS, Release on the Ionization Constant of H2O,
      http://www.iapws.org/relguide/Ionization.pdf