docs: def get_mean(self, C, rup, sites, dists):
          """
          Returns the mean ground motion in terms of log10 m/s/s, implementing
          equation 2 (page 502)
          """
          # W2 needs to be a 1 by 5 matrix (not a vector
          w_2 = np.array([
              [C["W_21"], C["W_22"], C["W_23"], C["W_24"], C["W_25"]]
              ])
          # Gets the style of faulting dummy variable
          sof = self._get_sof_dummy_variable(rup.rake)
          # Get the normalised coefficients
          p_n = self.get_pn(rup, sites, dists, sof)
          mean = np.zeros_like(dists.rjb)
          # Need to loop over sites - maybe this can be improved in future?
          # ndenumerate is used to allow for application to 2-D arrays
          for idx, rval in np.ndenumerate(p_n[0]):
              # Place normalised coefficients into a single array
              p_n_i = np.array([rval, p_n[1], p_n[2][idx], p_n[3], p_n[4]])
              # Executes the main ANN model
              mean_i = np.dot(w_2, np.tanh(np.dot(self.W_1, p_n_i) + self.B_1))
              mean[idx] = (0.5 * (mean_i + C["B_2"] + 1.0) *
                           (C["tmax"] - C["tmin"])) + C["tmin"]
          return mean