docs: def plate_rate_mc(pole1_plon, pole1_plat, pole1_kappa, pole1_N, pole1_age, pole1_age_error,
                    pole2_plon, pole2_plat, pole2_kappa, pole2_N, pole2_age, pole2_age_error,
                    ref_loc_lon, ref_loc_lat, samplesize=10000, random_seed=None, plot=True,
                    savefig=True, save_directory='./', figure_name=''):
      """
      Determine the latitudinal motion implied by a pair of poles and utilize
      the Monte Carlo sampling method of Swanson-Hysell (2014) to determine the
      associated uncertainty.
  
      Parameters:
      ------------
      plon : longitude of pole
      plat : latitude of pole
      kappa : Fisher precision parameter for VPGs in pole
      N : number of VGPs in pole
      age : age assigned to pole in Ma
      age_error : 1 sigma age uncertainty in million years
      ref_loc_lon : longitude of reference location
      ref_loc_lat : latitude of reference location
      samplesize : number of draws from pole and age distributions (default set to 10000)
      random_seed : set random seed for reproducible number generation (default is None)
      plot : whether to make figures (default is True, optional)
      savefig : whether to save figures (default is True, optional)
      save_directory = default is local directory (optional)
      figure_name = prefix for file names (optional)
  
      Returns
      --------
      rate : rate of latitudinal motion in cm/yr along with estimated 2.5 and 97.5
      percentile rate estimates
      """
      ref_loc = [ref_loc_lon, ref_loc_lat]
      pole1 = (pole1_plon, pole1_plat)
      pole1_paleolat = 90 - pmag.angle(pole1, ref_loc)
      pole2 = (pole2_plon, pole2_plat)
      pole2_paleolat = 90 - pmag.angle(pole2, ref_loc)
      print("The paleolatitude for ref_loc resulting from pole 1 is:" +
            str(pole1_paleolat))
      print("The paleolatitude for ref_loc resulting from pole 2 is:" +
            str(pole2_paleolat))
      rate = old_div(((pole1_paleolat - pole2_paleolat) * 111 *
                      100000), ((pole1_age - pole2_age) * 1000000))
      print("The rate of paleolatitudinal change implied by the poles pairs in cm/yr is:" + str(rate))
  
      if random_seed != None:
          np.random.seed(random_seed)
      pole1_MCages = np.random.normal(pole1_age, pole1_age_error, samplesize)
      pole2_MCages = np.random.normal(pole2_age, pole2_age_error, samplesize)
  
      plt.hist(pole1_MCages, 100, histtype='stepfilled',
               color='darkred', label='Pole 1 ages')
      plt.hist(pole2_MCages, 100, histtype='stepfilled',
               color='darkblue', label='Pole 2 ages')
      plt.xlabel('Age (Ma)')
      plt.ylabel('n')
      plt.legend(loc=3)
      if savefig == True:
          plot_extension = '_1.svg'
          plt.savefig(save_directory + figure_name + plot_extension)
      plt.show()
  
      pole1_MCpoles = []
      pole1_MCpole_lat = []
      pole1_MCpole_long = []
      pole1_MCpaleolat = []
      for n in range(samplesize):
          vgp_samples = []
          for vgp in range(pole1_N):
              # pmag.dev returns a direction from a fisher distribution with
              # specified kappa
              direction_atN = pmag.fshdev(pole1_kappa)
              # this direction is centered at latitude of 90 degrees and needs to be rotated
              # to be centered on the mean pole position
              tilt_direction = pole1_plon
              tilt_amount = 90 - pole1_plat
              direction = pmag.dotilt(
                  direction_atN[0], direction_atN[1], tilt_direction, tilt_amount)
              vgp_samples.append([direction[0], direction[1], 1.])
          mean = pmag.fisher_mean(vgp_samples)
          mean_pole_position = (mean['dec'], mean['inc'])
          pole1_MCpoles.append([mean['dec'], mean['inc'], 1.])
          pole1_MCpole_lat.append(mean['inc'])
          pole1_MCpole_long.append(mean['dec'])
          paleolat = 90 - pmag.angle(mean_pole_position, ref_loc)
          pole1_MCpaleolat.append(paleolat[0])
  
      pole2_MCpoles = []
      pole2_MCpole_lat = []
      pole2_MCpole_long = []
      pole2_MCpaleolat = []
      for n in range(samplesize):
          vgp_samples = []
          for vgp in range(pole2_N):
              # pmag.dev returns a direction from a fisher distribution with
              # specified kappa
              direction_atN = pmag.fshdev(pole2_kappa)
              # this direction is centered at latitude of 90 degrees and needs to be rotated
              # to be centered on the mean pole position
              tilt_direction = pole2_plon
              tilt_amount = 90 - pole2_plat
              direction = pmag.dotilt(
                  direction_atN[0], direction_atN[1], tilt_direction, tilt_amount)
              vgp_samples.append([direction[0], direction[1], 1.])
          mean = pmag.fisher_mean(vgp_samples)
          mean_pole_position = (mean['dec'], mean['inc'])
          pole2_MCpoles.append([mean['dec'], mean['inc'], 1.])
          pole2_MCpole_lat.append(mean['inc'])
          pole2_MCpole_long.append(mean['dec'])
          paleolat = 90 - pmag.angle(mean_pole_position, ref_loc)
          pole2_MCpaleolat.append(paleolat[0])
  
      if plot is True:
          plt.figure(figsize=(5, 5))
          map_axis = make_mollweide_map()
  
          plot_vgp(map_axis, pole1_MCpole_long, pole1_MCpole_lat, color='b')
          plot_vgp(map_axis, pole2_MCpole_long, pole2_MCpole_lat, color='g')
          if savefig == True:
              plot_extension = '_2.svg'
              plt.savefig(save_directory + figure_name + plot_extension)
          plt.show()
  
      # calculating the change in paleolatitude between the Monte Carlo pairs
      pole1_pole2_Delta_degrees = []
      pole1_pole2_Delta_kilometers = []
      pole1_pole2_Delta_myr = []
      pole1_pole2_degrees_per_myr = []
      pole1_pole2_cm_per_yr = []
  
      for n in range(samplesize):
          Delta_degrees = pole1_MCpaleolat[n] - pole2_MCpaleolat[n]
          Delta_Myr = pole1_MCages[n] - pole2_MCages[n]
          pole1_pole2_Delta_degrees.append(Delta_degrees)
          degrees_per_myr = old_div(Delta_degrees, Delta_Myr)
          cm_per_yr = old_div(((Delta_degrees * 111) * 100000),
                              (Delta_Myr * 1000000))
          pole1_pole2_degrees_per_myr.append(degrees_per_myr)
          pole1_pole2_cm_per_yr.append(cm_per_yr)
  
      if plot is True:
          plotnumber = 100
          plt.figure(num=None, figsize=(10, 4))
          plt.subplot(1, 2, 1)
          for n in range(plotnumber):
              plt.plot([pole1_MCpaleolat[n], pole2_MCpaleolat[n]],
                       [pole1_MCages[n], pole2_MCages[n]], 'k-', linewidth=0.1, alpha=0.3)
          plt.scatter(pole1_MCpaleolat[:plotnumber],
                      pole1_MCages[:plotnumber], color='b', s=3)
          plt.scatter(pole1_paleolat, pole1_age, color='lightblue',
                      s=100, edgecolor='w', zorder=10000)
          plt.scatter(pole2_MCpaleolat[:plotnumber],
                      pole2_MCages[:plotnumber], color='g', s=3)
          plt.scatter(pole2_paleolat, pole2_age, color='lightgreen',
                      s=100, edgecolor='w', zorder=10000)
          plt.plot([pole1_paleolat, pole2_paleolat], [
                   pole1_age, pole2_age], 'w-', linewidth=2)
          plt.gca().invert_yaxis()
          plt.xlabel('paleolatitude (degrees)', size=14)
          plt.ylabel('time (Ma)', size=14)
  
          plt.subplot(1, 2, 2)
          plt.hist(pole1_pole2_cm_per_yr, bins=600)
          plt.ylabel('n', size=14)
          plt.xlabel('latitudinal drift rate (cm/yr)', size=14)
          # plt.xlim([0,90])
          if savefig == True:
              plot_extension = '_3.svg'
              plt.savefig(save_directory + figure_name + plot_extension)
          plt.show()
  
      twopointfive_percentile = stats.scoreatpercentile(
          pole1_pole2_cm_per_yr, 2.5)
      fifty_percentile = stats.scoreatpercentile(pole1_pole2_cm_per_yr, 50)
      ninetysevenpointfive_percentile = stats.scoreatpercentile(
          pole1_pole2_cm_per_yr, 97.5)
      print("2.5th percentile is: " +
            str(round(twopointfive_percentile, 2)) + " cm/yr")
      print("50th percentile is: " + str(round(fifty_percentile, 2)) + " cm/yr")
      print("97.5th percentile is: " +
            str(round(ninetysevenpointfive_percentile, 2)) + " cm/yr")
      return rate[0], twopointfive_percentile, ninetysevenpointfive_percentile