code: def _calc_mod_reduc(self, strains, strain_ref, x_1, x_1_mean, x_2,
                          x_2_mean, x_3, x_3_mean):
          
  
          ones = np.ones_like(strains)
          # Predictor
          x_4 = np.log(self._lab_consol_ratio) * ones
          x = np.c_[ones, x_1, x_2, x_3, x_4, (x_1 - x_1_mean) * (
              x_2 - x_2_mean), (x_1 - x_1_mean) * (x_3 - x_3_mean), (
                  x_2 - x_2_mean) * (x_3 - x_3_mean), (x_1 - x_1_mean) * (
                      x_2 - x_2_mean) * (x_3 - x_3_mean)]
          # Coefficients
          denom = np.log(1 / strain_ref + strains / strain_ref)  # TODO: is this percent or decimal?
          b = np.c_[5.11 * ones, -0.729 * ones, (1 - 0.37 * x_3_mean * (1 + ((
              np.log(strain_ref) - x_1_mean) / denom))), -0.693 * ones, 0.8 - 0.4
                    * x_3, 0.37 * x_3_mean / denom, 0.0 * ones, -0.37 * (1 + (
                        np.log(strain_ref) - x_1_mean) / denom), 0.37 / denom, ]
          ln_shear_mod = (b * x).sum(axis=1)
          shear_mod = np.exp(ln_shear_mod)
          mod_reduc = shear_mod / shear_mod[0]
          return mod_reduc
docs: Compute the shear modulus reduction using Equation (1).