docs: def incremental_value(self, slip_moment, mmax, mag_value, bbar, dbar):
          """
          Returns the incremental rate of earthquakes with M = mag_value
          """
          delta_m = mmax - mag_value
          dirac_term = np.zeros_like(mag_value)
          dirac_term[np.fabs(delta_m) < 1.0E-12] = 1.0
          a_1 = self._get_a1(bbar, dbar, slip_moment, mmax)
          return a_1 * (bbar * np.exp(bbar * delta_m) * (delta_m > 0.0)) +\
              a_1 * dirac_term