docs: def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
          """
          See :meth:`superclass method
          <.base.GroundShakingIntensityModel.get_mean_and_stddevs>`
          for spec of input and result values.
          """
          # Check that the requested standard deviation type is available
          assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                     for stddev_type in stddev_types)
          #
          # Set parameters
          magn = rup.mag
          epi = dists.repi
          theta = dists.azimuth
          #
          # Convert Mw into Ms
          if magn < 6.58:
              mag = (magn - 0.59) / 0.86
          else:
              mag = (magn + 2.42) / 1.28
          #
          # Set coefficients
          coeff = self.COEFFS[imt]
          a1ca, a1cb, a1cc, a1cd, a1ce, a2ca, a2cb, a2cc, a2cd, a2ce = \
              gc(coeff, mag)
          #
          # Get correction coefficients. Here for each site we find the
          # the geometry of the ellipses
          ras = []
          for epi, theta in zip(dists.repi, dists.azimuth):
              res = get_ras(epi, theta, mag, coeff)
              ras.append(res)
          ras = np.array(ras)
          rbs = rbf(ras, coeff, mag)
          #
          # Compute values of ground motion for the two cases. The value of
          # 225 is hardcoded under the assumption that the hypocentral depth
          # corresponds to 15 km (i.e. 15**2)
          mean1 = (a1ca + a1cb * mag +
                   a1cc * np.log((ras**2+225)**0.5 +
                                 a1cd * np.exp(a1ce * mag)))
          mean2 = (a2ca + a2cb * mag +
                   a2cc * np.log((rbs**2+225)**0.5 +
                                 a2cd * np.exp(a2ce * mag)))
          #
          # Get distances
          x = (mean1 * np.sin(np.radians(dists.azimuth)))**2
          y = (mean2 * np.cos(np.radians(dists.azimuth)))**2
          mean = mean1 * mean2 / np.sqrt(x+y)
          if imt.name == "PGA":
              mean = np.exp(mean)/g/100
          elif imt.name == "PGV":
              mean = np.exp(mean)
          else:
              raise ValueError('Unsupported IMT')
          #
          # Get the standard deviation
          stddevs = self._compute_std(coeff, stddev_types, len(dists.repi))
          #
          # Return results
          return np.log(mean), stddevs