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.
          """
  
          # Distance term
          R = np.sqrt(dists.rjb ** 2 + 11.29 ** 2)
  
          # Magnitude term
          M = rup.mag - 6
  
          # Site term only distinguishes between lava and ash;
          # since ash sites have Vs30 in the range 60-200m/s,
          # we use this upper value as class separator
          S = np.zeros(R.shape)
          S[sites.vs30 <= 200] = 1
  
          # Mean ground motion (log10)
          mean = (0.518 + 0.387*M - np.log10(R) - 0.00256*R + 0.335*S)
  
          # Converting to natural log
          mean /= np.log10(np.e)
  
          # Check for standard deviation type
          assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                     for stddev_type in stddev_types)
  
          # Constant (total) standard deviation
          stddevs = [0.237/np.log10(np.e) + np.zeros(R.shape)]
  
          return mean, stddevs