docs: def continuum_seismicity(self, threshold_moment, e1h, e2h, err,
                               region_params):
          '''
          Function to implement the continuum seismicity calculation given
          vectors of input rates e1h, e2h [np.ndarray] and a dictionary of
          the corresponding regionalisation params
          returns a vector of the corresponding seismicity rates
          Python implementation of the CONTINUUM_SEISMICITY subroutine of
          SHIFT_GSRM.f90
  
          :param float threshold_moment:
              Target moment for calculation of activity rate
  
          :param np.ndarray e1h:
              First principal strain rate
  
          :param np.ndarray e1h:
              Second principal strain rate
  
          :param np.ndarray err:
              Vertical strain rate
  
          :param dict region_params:
              Activity rate parameters specific to the tectonic region under
              consideration
  
          :returns:
              Cumulative seismicity rate greater than or equal to the
              threshold magnitude
          '''
  
          strain_values = np.column_stack([e1h, e2h, err])
          e1_rate = np.amin(strain_values, axis=1)
          e3_rate = np.amax(strain_values, axis=1)
          e2_rate = 0. - e1_rate - e3_rate
          # Pre-allocate seismicity rate with zeros
          seismicity_rate = np.zeros(
              [np.shape(strain_values)[0], len(threshold_moment)],
              dtype=float)
          # Calculate moment rate per unit area
          temp_e_rate = 2.0 * (-e1_rate)
          id0 = np.where(e2_rate < 0.0)[0]
          temp_e_rate[id0] = 2.0 * e3_rate[id0]
          M_persec_per_m2 = (
              region_params['assumed_mu'] * temp_e_rate *
              region_params['coupled_thickness'])
  
          # Calculate seismicity rate at the threshold moment of the CMT
          # catalogue - Eq 6 in Bird et al (2010)
          seismicity_at_cmt_threshold = region_params['CMT_pure_event_rate'] * \
              (M_persec_per_m2 / region_params['tGR_moment_rate'])
          # Adjust forecast rate to desired rate using tapered G-R model
          # Taken from Eq 7 (Bird et al. 2010) and Eq 9 (Bird & Kagan, 2004)
          for iloc, moment_thresh in enumerate(threshold_moment):
              g_function = calculate_taper_function(
                  region_params['CMT_moment'],
                  moment_thresh,
                  region_params['corner_moment'],
                  region_params['beta'])
              seismicity_rate[:, iloc] = g_function * seismicity_at_cmt_threshold
          return seismicity_rate