docs: def plot_cumulative_moment(year, mag, figure_size=(8, 6),
                             filename=None, filetype='png', dpi=300, ax=None):
      '''Calculation of Mmax using aCumulative Moment approach, adapted from
      the cumulative strain energy method of Makropoulos & Burton (1983)
      :param year: Year of Earthquake
      :type year: numpy.ndarray
      :param mag: Magnitude of Earthquake
      :type mag: numpy.ndarray
      :keyword iplot: Include cumulative moment plot
      :type iplot: Boolean
      :return mmax: Returns Maximum Magnitude
      :rtype mmax: Float
      '''
      # Calculate seismic moment
      m_o = 10. ** (9.05 + 1.5 * mag)
      year_range = np.arange(np.min(year), np.max(year) + 1, 1)
      nyr = np.int(np.shape(year_range)[0])
      morate = np.zeros(nyr, dtype=float)
      # Get moment release per year
      for loc, tyr in enumerate(year_range):
          idx = np.abs(year - tyr) < 1E-5
          if np.sum(idx) > 0:
              # Some moment release in that year
              morate[loc] = np.sum(m_o[idx])
      ave_morate = np.sum(morate) / float(nyr)
  
      # Average moment rate vector
      exp_morate = np.cumsum(ave_morate * np.ones(nyr))
  
      if ax is None:
          fig, ax = plt.subplots(figsize=figure_size)
      else:
          fig = ax.get_figure()
  
      ax.step(year_range, np.cumsum(morate), 'b-', linewidth=2)
      ax.plot(year_range, exp_morate, 'r-', linewidth=2)
      # Get offsets
      upper_morate = exp_morate + (np.max(np.cumsum(morate) - exp_morate))
      lower_morate = exp_morate + (np.min(np.cumsum(morate) - exp_morate))
      ax.plot(year_range, upper_morate, 'r--', linewidth=1)
      ax.plot(year_range, lower_morate, 'r--', linewidth=1)
      ax.axis([np.min(year), np.max(year), 0.0, np.sum(morate)])
      _save_image(fig, filename, filetype, dpi)