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 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)
          '''
          # Working in Nm so convert:  shear_modulus - GPa -> Nm
          # area - km ** 2. -> m ** 2.
          # slip - mm/yr -> m/yr
          moment_rate = (shear_modulus * 1.E9) * (area * 1.E6) * (slip / 1000.)
          moment_mag = _scale_moment(self.mmax, in_nm=True)
          beta = self.b_value * log(10.)
          mag = np.arange(self.mmin - (self.bin_width / 2.),
                          self.mmax + self.bin_width,
                          self.bin_width)
          if self.b_value > 1.5:
              print('b-value larger than 1.5 will produce invalid results in '
                    'Anderson & Luco models')
              self.occurrence_rate = np.nan * np.ones(len(mag) - 1)
              return self.mmin, self.bin_width, self.occurrence_rate
  
          self.occurrence_rate = np.zeros(len(mag) - 1, dtype=float)
          for ival in range(0, len(mag) - 1):
              self.occurrence_rate[ival] = (
                  self.cumulative_value(mag[ival], moment_rate, beta, moment_mag)
                  - self.cumulative_value(
                      mag[ival + 1], moment_rate, beta, moment_mag))
  
          return self.mmin, self.bin_width, self.occurrence_rate