docs: def get_mfd(self, slip, area, shear_modulus=30.0):
          '''
          Calculates activity rate on the fault
  
          :param float slip:
              Slip rate in mm/yr
  
          :param fault_area:
              Width of the fault (km)
  
          :param float shear_modulus:
              Shear modulus of the fault (GPa)
  
          :returns:
              * Minimum Magnitude (float)
              * Bin width (float)
              * Occurrence Rates (numpy.ndarray)
          '''
  
          # Convert shear modulus GPa -> dyne-cm, area km ** 2 -> cm ** 2 and
          # slip mm/yr -> cm/yr
          slip_moment = (shear_modulus * 1E10) * (area * 1E10) * (slip / 10.)
          dbar = D_VALUE * np.log(10.0)
          bbar = self.b_value * np.log(10.0)
  
          mags = np.arange(self.mmin - (self.bin_width / 2.),
                           self.mmax + self.bin_width,
                           self.bin_width)
          if bbar >= dbar:
              print('b-value larger than 1.5 will produce invalid results in '
                    'Anderson & Luco models')
              self.occurrence_rate = np.nan * np.ones(len(mags) - 1)
              return self.mmin, self.bin_width, self.occurrence_rate
          self.occurrence_rate = np.zeros(len(mags) - 1, dtype=float)
          for ival in range(0, len(mags) - 1):
              self.occurrence_rate[ival] = \
                  RECURRENCE_MAP[self.mfd_type].cumulative_value(
                      slip_moment, self.mmax, mags[ival], bbar, dbar) - \
                  RECURRENCE_MAP[self.mfd_type].cumulative_value(
                      slip_moment, self.mmax, mags[ival + 1], bbar, dbar)
          return self.mmin, self.bin_width, self.occurrence_rate