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
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 = '' ) : 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 ) : direction_atN = pmag . fshdev ( pole1_kappa ) 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 ) : direction_atN = pmag . fshdev ( pole2_kappa ) 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 ( ) 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 ) 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
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