# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# +
###Made by <NAME>, Geoscience Dept, UW-Madison
###To process Iolite baseline-subtracted NP-II files for U-Th-Pb LA(SS)-ICPMS.

# +
### Makes Jupyter Notebook full width ###

from IPython.core.display import display, HTML
display(HTML("<style>.container { width:80% !important; }</style>"))
# -

# Mod 11-37-2021:
#
# -Make sure that your sample names for MKED-1 and/or MudTank match the code!
#
# -Make sure that your analysis order is correct
#
# To Do:
#
# -Write function to chop lines up into smaller intervals, to test how long of a line we really need.
#

# +
import pandas as pd
import os
import re
import copy
import numpy as np



from sklearn.linear_model import LinearRegression
import math
from math import log10, floor
from scipy import stats
from scipy import optimize
from scipy.stats import linregress
import webbrowser

#Graphing stuff
from matplotlib import pyplot as plt
import matplotlib.patches as mpatches
# %matplotlib inline

# #%pip install PyPDF2
from PyPDF2 import PdfFileMerger, PdfFileReader

# #%pip install pdfkit
import pdfkit

# #%pip install xlsxwriter
import xlsxwriter

#pd.set_option("display.precision", 8)

# +
#Constants and published values

# Initial coefficient

l238U = 1.55125 * 10 ** (-10)  # lambda_238U Jaffey et al., (1971); Half-life = 4.4683 ± 0.0024 [10^9 years]
l235U = 9.8485 * 10 ** (-10)  # lambda_235U Jaffey et al., (1971); Half-life = 7.0381 ± 0.0048 [10^8 years]
l232Th = 4.9475 * 10 ** (-11)  # lambda_232Th 
#l231Pa = 2.10 * 10 ** (-5)  # Rempfer2017epsl
#l230Th = 9.17 * 10 ** (-6)  # Cheng2013epsl
U85r = 137.818  # 238U/235U; Hiess et al., (????) Default in IsoplotR, but the Isoplot default is usually 137.88


#Published values dictionary MAY NEED TO ADD 207/235 -KD 10/18/21

# These ratios need to be checked for accuracy, have reference added, and clarified for whether ratio is corrected or uncorrected for initial Pb. -KD 11/3/2021


published_ratios_610 = {
    '206/238': 0.258174418887103,
    '208/232': 0.546910763260703,
    '207/206': 0.909778846717897,
    '207/235': (0.258174418887103) * (0.909778846717897) * U85r, #Calculated
    '208/206': 2.16900334369684,
    '206/204': 17.047,
    '207/204': 15.509,
    '208/204': 36.975,

}
published_ratios_612 = {
    '206/238': 0.289090155580635,
    '208/232': 0.598868250924868,
    '207/206': 0.907335907335907,
    '207/235': (0.289090155580635) * (0.907335907335907) * U85r, #Calculated
    '208/206': 2.16450216450216,
    '206/204': 17.094,
    '207/204': 15.51,
    '208/204': 37,      
}
published_ratios_BHVO = {
    '206/238': 1.24298582359808,
    '208/232': 0.811596334965975,
    '207/206': 0.831206960977953,
    '207/235': (1.24298582359808) * (0.831206960977953) * U85r, #Calculated
    '208/206': 2.042918913147926,
    '206/204': 18.733,
    '207/204': 15.571,
    '208/204': 38.27,      
}   
published_ratios_BCR_2G = {
    '206/238': 1.90697056349179,
    '208/232': 1.09082521297421,
    '207/206': 0.832720490274447,
    '207/235': (1.90697056349179) * (0.832720490274447) * U85r, #Calculated
    '208/206': 2.06394884092726,
    '206/204': 18.765,
    '207/204': 15.626,
    '208/204': 38.73,      
}
published_ratios_GSD_1G = {
    '206/238': 0.367307298285922,
    '208/232': 0.706245210573853,
    '207/206': 0.80417794575821,
    '207/235': (0.367307298285922) * (0.80417794575821) * U85r, #Calculated
    '208/206': 1.98723121712038,
    '206/204': 19.579,
    '207/204': 15.745,
    '208/204': 38.908,      
}

#Mod 9-1-2021: Changed GSE-1G values and MKED values
published_ratios_GSE_1G = { 
    '206/238': 0.367307298285922, ###These are just the GSD-1G ratios
    '208/232': 0.706245210573853, ###These are just the GSD-1G ratios
    '207/206': 0.79232,
    '208/206': 1.96460,
    '206/204': 19.9250,
    '207/204': 15.7870,
    '208/204': 39.1450,      
}
#Mod 9-13-2021: Change these to measured values.
published_ratios_MKED = {
    '206/238': 0.265752478, #Total Sample from Spandler et al., 2015
    #'206/238': 0.265752478 * 1.02, # Testing value -KD
    
    
    #'206/238': 0.2634133, #ID-ICPMS (unpublished)
    '208/232': 0.079191124, #ID-ICPMS (unpublished)
    '207/206': 0.096, #Total Sample from Spandler et al., 2015
    #'207/206': 0.095909435, #ID-ICPMS (unpublished)
     '207/235': (0.265752478) * (0.096) * U85r, #Calculated with Spandler et al., 2015 values
    '208/206': 0.852420222, #ID-ICPMS (unpublished)
    #'206/204': 9090.909, #Total Sample from Spandler et al., 2015
    '206/204': 10474.84592, #ID-ICPMS (unpublished)
    '207/204': 1004.661031, #ID-ICPMS (unpublished)
    '208/204': 8928.948295, #ID-ICPMS  (unpublished)   
}    
published_ratios_BLR_1 = {
    '206/238': 0.176,
    '208/232': 0.0533134203278787,
    '207/206': 0.0743078,
    '207/235': (0.176364) * (0.0743078) * U85r, #Calculated
    '208/206': 0.17982,
    '206/204': 975.36,
    '207/204': 72.47686,
    '208/204': 175.38924, 
}    

published_ratios_BLR_1a = { #Arizona Laserchron values for 'BLS'
    '206/238': 0.1769, #Arizona Laserchron values for 'BLS'
    '208/232': 0.05332, #Arizona Laserchron values for 'BLS'
    '207/206': 0.07431629, #Arizona Laserchron values for 'BLS'
    '207/235': (0.1769) * (0.07431629) * U85r, #Calculated #Arizona Laserchron values for 'BLS'
    '208/206': 0.015720513, # Calculated using U/Th ratio from Arizona... ?
    '206/204': 975.36, # Aleinikoff et al., (2007) ?
    '207/204': 72.47686, # Aleinikoff et al., (2007) ?
    '208/204': 175.38924, # Aleinikoff et al., (2007) ?     
    
    
}
published_ratios_BLR_1b = { #KD MOD
    '206/238': 0.1794688, 
    '208/232': 0.05332, #Arizona Laserchron values for 'BL
    '207/206': 0.08714508,#Arizona Laserchron values for 'BLS'
    '207/235': (0.1794688) * (0.08714508) * U85r,
    '208/206': 0.015720513, # Calculated using U/Th ratio from Arizona... ?
    '206/204': 975.36, # Aleinikoff et al., (2007) ?
    '207/204': 72.47686, # Aleinikoff et al., (2007) ?
    '208/204': 175.38924, # Aleinikoff et al., (2007) ?     
    
    
}
published_ratios_OLT_1 = {
    '206/238': 0.170681666666667,
    '208/232': 0.0516589274888366,
    '207/206': 0.0731801314551979,
    '207/235': (0.170681666666667) * (0.0731801314551979) * U85r, #Calculated
    '208/206': 0.58025,
    '206/204': 1847.16667,
    '207/204': 135.17590,
    '208/204': 1071.81688,      
}
published_ratios_OLT_2 = {
    '206/238': 0.16693,
    '208/232': 0.05064,
    '207/206': 0.07248,
    '207/235': (0.170681666666667) * (0.0731801314551979) * U85r, #Calculated
    '208/206': 0.63815,
    '206/204': 1453.00000,
    '207/204': 105.31206,
    '208/204': 927.22713,      
}
published_ratios_TCB = {
    '206/238': 0.17086,
    '208/232': 0.05178,
    '207/206': 0.07327,
    '207/235': (0.17086) * (0.07327) * U85r, #Calculated
    '208/206': 0.63348,
    '206/204': 2016.75000,
    '207/204': 147.76034,
    '208/204': 1277.56162,      
}
published_ratios_MudTank = {
    '206/238': 0.052269449,
    '208/232': 0.049381079,
    '207/206': 0.071923689,
    '207/235': (0.052269449) * (0.071923689) * U85r, #Calculated
    '208/206': 0.068164043,
    '206/204': 766.0145618,
    '207/204': 55.09459,
    '208/204': 52.21465,      
                                
}
                                
###Needed to add stuff for apatite
#WRONG DO NOT USE. These are not correct at all.
published_ratios_MT_apa = {
    '206/238': 0.052269449,
    '208/232': 0.049381079,
    '207/206': 0.071923689,
    '208/206': 0.068164043,
    '206/204': 766.0145618,
    '207/204': 55.09459,
    '208/204': 52.21465,      
}
published_ratios_MT_apa = {
    '206/238': 0.052269449,
    '208/232': 0.049381079,
    '207/206': 0.071923689,
    '208/206': 0.068164043,
    '206/204': 766.0145618,
    '207/204': 55.09459,
    '208/204': 52.21465,      
}

published_ratios_dict = {
    'SRM NIST 610': published_ratios_610,
    'SRM NIST 612': published_ratios_612,
    'BHVO-2G': published_ratios_BHVO,
    'BCR-2G': published_ratios_BCR_2G,
    'GSD-1G': published_ratios_GSD_1G,
    'GSE-1G': published_ratios_GSE_1G,
    'Bear Lake': published_ratios_BLR_1,
    'MKED-1': published_ratios_MKED,
    'Mud Tank': published_ratios_MudTank,
    'OLT-1': published_ratios_OLT_1,
    'OLT-2': published_ratios_OLT_2,
    'TCB': published_ratios_TCB,
    'BLR-1': published_ratios_BLR_1b
}

published_df = pd.DataFrame(published_ratios_dict)
published_df = pd.DataFrame.transpose(published_df)

def published(std_name, ratio):
    '''Returns the published ratio for a given reference material from the published values dataframe. Be sure that these are the correct values if you are using this function.'''
    return published_df.loc[std_name][ratio]

def published_full(std_name):
    '''Returns dictionary of published ratios for the input reference material. Be sure that these are the correct values if you are using this function.'''
    return published_ratios_dict[std_name]


# + code_folding=[7, 12, 17, 22, 33, 38, 44]
#Functions for Age Calculation: calc_t75(), calc_t68(), funct_t76(), calc_t76(), calc_t82(), calc_t76_2()

##### Credit should be given to: https://github.com/anoda/UPbplot.py/blob/master/UPbplot.py code #####

# ------------------------------------------------
# Age calculation for isotopic ratios (no correction)

def calc_t75(r):
    """ Calculation for 207Pb/235U age."""
    t = 1 / l235U * np.log(r + 1)
    return t

def calc_t68(r):         #PASSED
    """ Calculation for 207Pb/235U age."""
    t = 1 / l238U * np.log(r + 1)
    return t

def func_t76(t, r):
    """ Function for calc_t76."""
    res = abs(U85r * r - (np.exp(l235U * t) - 1) / (np.exp(l238U * t) - 1))
    return res

def calc_t76(r76):
    '''This takes a list in as input. -KD'''
    T = []
    for i, r in enumerate(r76):
        t = 1 / l238U * np.log(r + 1)  # initial time for calculation
        T76 = optimize.leastsq(func_t76, t, args=(r))[0][0]
        T.append(T76)
    return T

### KD modification###

def calc_t82(r):         
    """ Calculation for 208Pb/232U age. -KD"""
    t = 1 / l232Th * np.log(r + 1)
    return t

def calc_t76_2(r):
    '''Calculation of 207/206 Age from single ratio (float) inout. -KD'''
    t = 1 / l238U * np.log(r + 1)  # initial time for calculation
    T76 = optimize.leastsq(func_t76, t, args=(r))[0][0]
    return T76

def calc_concordant_ratios(mill_age):
    ''' Use the decay equations for U238, U235, and Th232 to calculate radiogenic (*) ratios for a given age [Ma].'''
    
    ratio_dict = {}
    age = mill_age * 1000000

    ratio_dict['207*/206*'] = abs ((1/U85r) * (np.exp(l235U * age) - 1) / (np.exp(l238U * age) - 1))
    ratio_dict['206*/238'] = np.exp(l238U * age) - 1
    ratio_dict['207*/235'] = np.exp(l235U * age) - 1
    ratio_dict['208*/232'] = np.exp(l232Th * age) - 1
    ratio_dict['238/206*'] = 1 / ratio_dict['206*/238']
    
    return ratio_dict


# IN PROGRESS 04Nov2021
# def funct_t86(t, r86, r68, r82):
#      """ Function for calc_t86."""
#     res = abs((r68*r86**1/r82) * r86 - (np.exp(l232Th * t) - 1) / (np.exp(l238U * t) - 1))
#     return res

# def calc_t86(r):
#     '''Calculation of 208/206 Age from single ratio (float) inout. -KD'''
#     t = 1 / l232Th * np.log(r + 1)  # initial time for calculation
#     T86 = optimize.leastsq(func_t86, t, args=(r))[0][0]
#     return T86

    

### END KD mod####

# + code_folding=[8]
### Functions for rounding significant figures and (unfinished) Stacey-Kramers corrections
   
def round_sig(x, sig=5):
    '''Rounds value to specified number of significant figures.'''
    return round(x, sig-int(floor(log10(abs(x))))-1)

#Stacey-Kramer

def SK_values(mill_age):
    '''Input age [Ma] younger than 3.7 ga and get the Stacey-Kramer values.'''
   
    age = mill_age * 1000000

    SK_dict = {}

    SK_dict['SK_206_204'] = round_sig(11.152 + 9.7357*( np.exp(l238U * 3700000000)- np.exp(l238U * age)))

    SK_dict['SK_207_204'] = round_sig(12.998 + (9.7357/137.88)* (np.exp(l235U * 3700000000)- np.exp(l235U * age)))

    SK_dict['SK_208_204'] = round_sig(31.23 + 36.84 * (np.exp(l232Th * 3700000000) - np.exp(l232Th * age)))

    SK_dict['SK_207_206'] = round_sig(SK_dict['SK_207_204'] / SK_dict['SK_206_204'])

    SK_dict['SK_208_206'] = round_sig(SK_dict['SK_208_204'] / SK_dict['SK_206_204'])


    return SK_dict

# ex./ SK_values(960)

# + code_folding=[18]
# INPUTS
# measured_76 = 0.0863356  # Measured 207/206

# initial_76 = 0.9161 # Initial 207/206; Stacey-Kramers, intercept, or independent...
# expected_76 = 0.0756 # The (radiogenic?) expected 207/206

# measured_68 = 0.183344 # Measured 206/238
# measured_82 = 5  # Measured 208/232
# measured_86 = 5  # Measured 208/206
# measured_75 = measured_68 * measured_76 * U85r  #'Measured' 207/235

# ratio_238_to_232 = 5 # Measured 238/232

# initial_86 = 1 # Initial 208/206; Stacey-Kramers, intercept, or independent...

####


def correction_207Pb(measured_68, measured_76,measured_82, measured_86, iterations = 5):
    '''Performs 5 iterations of 207-Pb correction using Stacey-Kramers.'''
    correction_dict = {}
    correction_dict['206/238_uncorr'] = measured_68
    for idx in range (iterations):
        if idx == 0:
            correction_dict['206/238_Age'] = calc_t68(correction_dict['206/238_uncorr']) / 1E6
        else:
            correction_dict['206/238_Age'] = calc_t68(correction_dict['206/238_corr']) / 1E6
            
        #correction_dict['206/238_Age'] = calc_t68(correction_dict['206/238_uncorr']) / 1E6
        correction_dict['207/206_uncorr'] = measured_76
        correction_dict['207/206_initial_SK'] = SK_values(correction_dict['206/238_Age'])['SK_207_206']
        correction_dict['207/206_rad_expected'] = calc_concordant_ratios(correction_dict['206/238_Age'])['207*/206*']
        correction_dict['f_206'] = (correction_dict['207/206_uncorr'] - correction_dict['207/206_rad_expected']) / (correction_dict['207/206_initial_SK'] - correction_dict['207/206_rad_expected'])
        correction_dict['206/238_corr'] = (1 - correction_dict['f_206']) * measured_68
      
    
    correction_dict['206/238_Age'] = calc_t68(correction_dict['206/238_corr']) / 1E6
    
    correction_dict['207/206_corr'] = (measured_76 - initial_76 * correction_dict['f_206']) / (1 - correction_dict['f_206']) 
    correction_dict['208/206_corr'] = (measured_86 - initial_86 * correction_dict['f_206']) / (1 - correction_dict['f_206']) 
    correction_dict['208/232_corr'] =  measured_82 - (measured_68 * initial_86 * ratio_238_to_232 * correction_dict['f_206'])
    correction_dict['207/235_corr'] = measured_75 - (measured_68 * U85r * initial_76 *  correction_dict['f_206'])
       
    return correction_dict


# + code_folding=[1, 33, 57]
### Hard-code function for grouping samples by sampleID, and order dataframe based on chronologic order: group_samples(), comb_groups(), chronologic()
def group_samples(data):
    '''Function for grouping up samples by Sample ID. Hard-coded for names with "1.n" in title.'''
    samples = data.index.values.tolist()
    group = None
    group_names = []
    start = [0]
    end = []
    for idx in range(len(samples)):
        if group == None:
            group = samples[idx].split(' 1')[0]     #Hard-coding
            group_names.append(group)
        #print(samples[idx])
        if group in samples[idx]:
            pass
        else:
            end.append(idx)

            start.append(idx)
            name = samples[idx].split(' 1')[0]      #Hard-coding
            group = name
            group_names.append(name)
    end.append(len(samples))

#     print('start:', start)
#     print('end:', end)
#     print('group names:', group_names)
    
    start_index = start
    end_index = end
    
    return start_index, end_index, group_names

def comb_groups(group_list, data):
    """ Function to combine a list of group names from dataset and output a list of rows to plot as well as the number of labels that will be necessary."""
    start, end, names = group_samples(data)
    #test_group = ['Bancroft', 'BLR-1']

    #len(test_group)
    indx_list = []
    for val in group_list:
        #print(val)
        indx_list.append(names.index(val))
    #print(indx_list)
    y_list = []
    for ind in indx_list:
        for i in range(end[ind]-start[ind]):
            y1 = data.iloc[start[ind] + i]
            y_list.append(y1)


    df3 = pd.DataFrame(y_list)

    labels = len(df3.index.tolist())

    return y_list, labels

def chronologic(result_df, order):
    '''This function will reorder the rows of a dataframe by the index order that you designate.
    In this case it is used to re-order dataframes by the Sample ID in chronologic order.'''
    
    new_result = result_df.copy(deep=True)
    return new_result.reindex(order)


# + code_folding=[]
### Functions for read_np2_timeseries(), calc_CPS(), and statistics_NP2()

def read_np2_timeseries(excel_file):
    ''' Excel input file is your baseline corrected time series export from Iolite for the NP-II.'''
    df = pd.read_excel(excel_file, sheet_name = None)
    keys = df.keys()
    header_row = 0
    new_dict = {}
    for key in keys:
        if 'time' in key: #Kind of hard-coded right now, so if names get weird may need to change
            
            df_test = df[key]

            df_test.columns = df_test.iloc[header_row]
            df_test = df_test.drop(header_row)
            df_test = df_test.reset_index(drop=True)
            
            
            new_string = key.split('time')[0].rstrip()
            new_dict[new_string] = df_test #test1_new
    return new_dict

def calc_CPS(np2_dict):
    """ Function for calculating CPS from baseline-corrected signal, OPZ (On-Peak Zero), SNR (Signal-to-Noise Ratio), and isotope ratios."""
    columns = ['Absolute Time',
     'Elapsed Time',
     'm238_CPS',
     'm232_CPS',
     'm208_CPS',
     'm207_CPS',
     'm206_CPS',
     'm204_CPS',
     'm202_CPS']

    new_col = ['Absolute Time',
     'Elapsed Time',
     '238_CPS',
     '232_CPS',
     '208_CPS',
     '207_CPS',
     '206_CPS',
     '204_CPS',
     '202_CPS']

    cut_col = ['238_CPS',
     '232_CPS',
     '208_CPS',
     '207_CPS',
     '206_CPS',
     '204_CPS',
     '202_CPS']

    calc_dict = {}
    for key in np2_dict:
        #print(key)
        test_df1 = np2_dict[key]

        for col in columns:
                test_df2 = test_df1.apply(lambda x: x * 62500000 if 'CPS' in x.name else x)
                test_df2 = test_df2[['Absolute Time',
             'Elapsed Time',
                 'm238_CPS',
                 'm232_CPS',
                 'm208_CPS',
                 'm207_CPS',
                 'm206_CPS',
                 'm204_CPS',
                 'm202_CPS',]]
        test_df2.columns = new_col
        test_df2 = test_df2[cut_col]
        result = pd.concat([test_df1, test_df2], axis=1)
        
         #Calculating OPZ
        result['OPZ_238'] = result.apply(lambda x: x['m238'] - x['m238_CPS'], axis=1)
        result['OPZ_232'] = result.apply(lambda x: x['m232'] - x['m232_CPS'], axis=1)
        result['OPZ_208'] = result.apply(lambda x: x['m208'] - x['m208_CPS'], axis=1)
        result['OPZ_207'] = result.apply(lambda x: x['m207'] - x['m207_CPS'], axis=1)
        result['OPZ_206'] = result.apply(lambda x: x['m206'] - x['m206_CPS'], axis=1)
        result['OPZ_204'] = result.apply(lambda x: x['m204'] - x['m204_CPS'], axis=1)
        result['OPZ_202'] = result.apply(lambda x: x['m202'] - x['m202_CPS'], axis=1)
        result['OPZ_208/206'] = result.apply(lambda x: x['OPZ_208'] / x['OPZ_206'], axis=1)
        result['OPZ_207/206'] = result.apply(lambda x: x['OPZ_207'] / x['OPZ_206'], axis=1)
        result['OPZ_206/204_Hg-corrected'] = result.apply(lambda x: x['OPZ_206'] / (x['OPZ_204'] - (x['OPZ_202'] * 6.87/29.86)) , axis=1)
       
        #Calculating Signal-to-Noise Ratios [V]/[V]
        result['SNR_238'] = result.apply(lambda x: x['m238_CPS']/ x['OPZ_238'], axis=1)
        result['SNR_232'] = result.apply(lambda x: x['m232_CPS']/ x['OPZ_232'], axis=1)
        result['SNR_208'] = result.apply(lambda x: x['m208_CPS']/ x['OPZ_208'], axis=1)
        result['SNR_207'] = result.apply(lambda x: x['m207_CPS']/ x['OPZ_207'], axis=1)
        result['SNR_206'] = result.apply(lambda x: x['m206_CPS']/ x['OPZ_206'], axis=1)
        result['SNR_204'] = result.apply(lambda x: x['m204_CPS']/ x['OPZ_204'], axis=1)
        result['SNR_202'] = result.apply(lambda x: x['m202_CPS']/ x['OPZ_202'], axis=1)
        
        #Calculating Ratios
        result['206/238'] = result.apply(lambda x: x['206_CPS']/x['238_CPS'], axis=1)
        result['208/232'] = result.apply(lambda x: x['208_CPS']/x['232_CPS'], axis=1)
        result['207/206'] = result.apply(lambda x: x['207_CPS']/x['206_CPS'], axis=1)
        result['208/206'] = result.apply(lambda x: x['208_CPS']/x['206_CPS'], axis=1)
        result['206/204'] = result.apply(lambda x: x['206_CPS']/x['204_CPS'], axis=1)
        result['208/204'] = result.apply(lambda x: x['208_CPS']/x['204_CPS'], axis=1)
        result['207/204'] = result.apply(lambda x: x['207_CPS']/x['204_CPS'], axis=1)
        ### KD added 10/18/2021
        result['207/235'] = result.apply(lambda x: x['207/206'] * x['206/238'] * U85r, axis=1)
        
        #Calculating OPZ CPS
        result['OPZ_238_CPS'] = result.apply(lambda x: x['OPZ_238'] * 62500000, axis=1)
        result['OPZ_232_CPS'] = result.apply(lambda x: x['OPZ_232'] * 62500000, axis=1)
        result['OPZ_208_CPS'] = result.apply(lambda x: x['OPZ_208'] * 62500000, axis=1)
        result['OPZ_207_CPS'] = result.apply(lambda x: x['OPZ_207'] * 62500000, axis=1)
        result['OPZ_206_CPS'] = result.apply(lambda x: x['OPZ_206'] * 62500000, axis=1)
        result['OPZ_204_CPS'] = result.apply(lambda x: x['OPZ_204'] * 62500000, axis=1)
        result['OPZ_202_CPS'] = result.apply(lambda x: x['OPZ_202'] * 62500000, axis=1)
        
        calc_dict[key] = result
    
    return calc_dict

def statistics_NP2(calc_dict):
    """Function for calculating statistics of CPS and isotope ratios, prior to error minimization and correction."""
    
    calc_list = ['238_CPS', '232_CPS',
           '208_CPS', '207_CPS', '206_CPS', '204_CPS', '202_CPS', '206/238',
           '208/232', '207/206', '208/206', '206/204','208/204','207/204', '207/235'  ]
    mega_dict = {}

    for sheet in calc_dict:
        tester = calc_dict[sheet]
        stats_dict = {}
        for col in tester:

            if col in calc_list:
                #print(col)
                if '/' in col:
                    key = col + '_before rejection'
                else:
                    key = col + '_mean'
                df_mean = tester[col].mean()
                stats_dict[key] = df_mean
                #This is not 2SE%, should probably fix labels... KD 14 June 2021
                df_precision = (2 * tester[col].sem()) / df_mean * 100
                
                stats_dict[col + '_2se%'] = df_precision
            if 'OPZ' in col:
                 stats_dict[col + '_mean'] = tester[col].mean()
            if 'SNR' in col:
                 stats_dict[col + '_mean'] = tester[col].mean() 
                    
        stats_dict['Time (s)'] = tester['Elapsed Time'].max()
        
        #new_string = sheet.replace('time series data', '')
        new_string = sheet.split('time')[0].rstrip()
        mega_dict[new_string] = stats_dict

    df_1 = pd.DataFrame(mega_dict)
    df_flip = pd.DataFrame.transpose(df_1)
    return df_flip


# + code_folding=[2, 25, 96]
### Functions for group_samples(), regress_group(), regress_indiv(), and regress_compare()

def group_samples2(tester):
    '''Takes in the output from calc_CPS() and outputs a dictionary with individual analyses grouped by sample.'''
    samples = tester.keys()
    group_list = []
    group_dict = {}
    for sample in samples:
        group = sample.split(' 1.')[0]

        if len(group_list) == 0:
            group_list.append(group)
        elif group == group_list[-1]:
            pass
        else:
            group_list.append(group)

    for group in group_list:
        group_dict[group] = []
        for sample in samples:
            if sample.split(' 1.')[0] == group:
                group_dict[group].append(sample)

    return group_dict 

def regress_group(tester, group_name,choice = False, start_clip = 0):
    '''Does a linear regression for each sample and outputs a dictionary with mean, intercept, slope, coefficient of determination, and (intercept/average).'''
    mked = group_samples2(tester)[group_name]

    tester_list = []
    for entry in mked:
        tester_list.append(tester[entry])

    mked_full = pd.concat(tester_list)
    regress_data = mked_full[['Elapsed Time', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]

    ### Regression Math ###
    #start_clip = 1.0
    mked_full = regress_data[regress_data['Elapsed Time'] > start_clip]

    x = np.array(mked_full['Elapsed Time']).reshape((-1, 1))
    y = mked_full['206/238']

    model = LinearRegression().fit(x, y)

    model_dict = {}
    model_dict['average'] = y.mean()
    model_dict['intercept'] = model.intercept_
    model_dict['slope'] = float(model.coef_)

    model_dict['coefficient of determination'] = model.score(x, y)
    
    new = model_dict['intercept']
    model_dict['intercept / average'] = new/model_dict['average']
   
    if choice:
        print('\nRegression for', group_name)
        print('intercept:', model.intercept_)
        print('slope:', float(model.coef_))
        print('coefficient of determination:', model.score(x, y))
        print('intercept / average', model_dict['intercept / average'] )

    return model_dict

def regress_indiv(tester, sample_name, choice = False, start_clip = 1, end_clip = 1):
    '''Does a linear regression for each individual analysis and outputs a dictionary with mean, intercept, slope, coefficient of determination, and (intercept/average).'''
    sample = tester[sample_name]
    regress_data = sample[['Elapsed Time', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]
    mked_full1 = regress_data[regress_data['Elapsed Time'] > start_clip]
    mked_full = mked_full1[mked_full1['Elapsed Time'] > end_clip]
    
    x = np.array(mked_full['Elapsed Time']).reshape((-1, 1))
    y = mked_full['206/238']

    model = LinearRegression().fit(x, y)

    model_dict = {}
    model_dict['average'] = y.mean()
    model_dict['intercept'] = model.intercept_
    model_dict['slope'] = float(model.coef_)

    model_dict['coefficient of determination'] = model.score(x, y)
    
    new = model_dict['intercept']
    model_dict['intercept / average'] = new / model_dict['average']
    
    
    if choice:
        print('\nRegression for', sample_name)
        print('average:', model_dict['average'])
        print('intercept:', new)
        print('slope:', float(model.coef_))
       #print('coefficient of determination:', model.score(x, y))
       #print('intercept / average', new/orig)

    return model_dict

def regress_compare(tester, std1, std2, clip1 = 0, clip2  = 0, endclip = 35, choice = False):    
    '''This function takes in the output from calc_CPS(), and two selected sample groups, plus clip positions for start and end. It will plot things.'''
    ###MKED###
    #start_clip = 2
    mked = group_samples2(tester)[std1]

    tester_list = []
    for entry in mked:
        tester_list.append(tester[entry])

    mked_full = pd.concat(tester_list)
    regress_data = mked_full[['Elapsed Time', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]

    x_orig = np.array(mked_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    mked_full = regress_data[regress_data['Elapsed Time'] > clip1]

    mked_full = mked_full[mked_full['Elapsed Time'] < endclip]
    
    x = np.array(mked_full['Elapsed Time']).reshape((-1, 1))
    y = mked_full['206/238']
    model = LinearRegression().fit(x, y)
    model_dict = {}
    model_dict['average'] = y.mean()
    model_dict['intercept'] = model.intercept_
    model_dict['slope'] = float(model.coef_)
    model_dict['coefficient of determination'] = model.score(x, y)

    fit = model.fit(x,y)
    y_pred = model.predict(x_orig)

    ### MKED ###

    fig, ax = plt.subplots(1, 3, sharex = True, figsize = (24, 6))
    ax[0].scatter(x,y, color="black")
    ax[0].plot(x_orig, y_pred, color='blue', linewidth = 3)
    
    # these are matplotlib.patch.Patch properties
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

    # place a text box in upper left in axes coords
    ax[0].text(0.97, 0.97, std1, transform=ax[0].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation = 'y = ' + str(round_sig(model_dict['slope'], 3)) + 'x + ' + str(round_sig(model_dict['intercept']))
    ax[0].text(0.97, 0.88, equation, transform=ax[0].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    

    ### BLR-1 ###
    #start_clip = 8
    blr = group_samples2(tester)[std2]

    tester_list = []
    for entry in blr:
        tester_list.append(tester[entry])

    blr_full = pd.concat(tester_list)
    regress_data = blr_full[['Elapsed Time', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]


    x_orig2 = np.array(blr_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    blr_full = regress_data[regress_data['Elapsed Time'] > clip2]
    blr_full = blr_full[blr_full['Elapsed Time'] < endclip]
    x2 = np.array(blr_full['Elapsed Time']).reshape((-1, 1))
    y2 = blr_full['206/238']
    model = LinearRegression().fit(x2, y2)

    model_dict2 = {}
    model_dict2['average'] = y2.mean()
    model_dict2['intercept'] = model.intercept_
    model_dict2['slope'] = float(model.coef_)
    model_dict2['coefficient of determination'] = model.score(x, y)


    fit2 = model.fit(x2,y2)
    y_pred2 = model.predict(x_orig2)

    ### BLR-1 ###
    ax[1].scatter(x2,y2, color="black")
    ax[1].plot(x_orig2, y_pred2, color='blue', linewidth = 3)

    # place a text box in upper left in axes coords
    ax[1].text(0.97, 0.97, std2, transform=ax[1].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation2 = 'y = ' + str(round_sig(model_dict2['slope'], 3)) + 'x + ' + str(round_sig(model_dict2['intercept']))
    ax[1].text(0.97, 0.88, equation2, transform=ax[1].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    
    y3 = y_pred2 / (y_pred[:len(y_pred2)])
    y4 = model_dict2['intercept'] / (y_pred[:len(y_pred2)])

    ax[2].plot(y3[0:50], color='green', linewidth = 3)
    ax[2].plot(y4[0:50], color='blue', linewidth = 3)
    
    
#     #Still in Development
#     ax[0].spines['left'].set_position(('data', 0))
#     ax[1].spines['left'].set_position(('data', 0))
#     ax[2].spines['left'].set_position(('data', 0))
    
    ###
    #print(equation)
    #print(model_dict)
    #print(model_dict2)
    print('int2/int1: ',model_dict2['intercept']/model_dict['intercept'])
    print('int2/avg1 ', model_dict2['intercept']/model_dict['average'])
    print(std1, ' chopped at: ',clip1, '   ', model_dict2['intercept']/y_pred[clip1])
    print()
    if choice:
        return model_dict, model_dict2
    


# + code_folding=[]
### Functions for error minimization (Could use some revision...): ranked_minimization2() and statistics_ranktest2()

### Getting rid of negatives

def ranked_minimization2(sheet, ratio, reject_percentage = 20):
    ''' Function for excluding negative isotope ratios and then performing a ranked error minimization by excluding a specified percentage of points that lie outside of 1SD of mean. Excludeds negative values BEFORE error minimization.'''
    
    mytest = tester[sheet].copy(deep=True)

    df_mean_before = mytest[ratio].mean()
    df_1std_before = mytest[ratio].std()
    df_count_before = mytest[ratio].count()
    df_2se_perc_before = (2 * mytest[ratio].sem()) / df_mean_before * 100

    dif_mean = ratio + '_dif_from_mean'
    dif_1SD = ratio + '_dif_from_1SD'
    mytest[dif_mean] = mytest.apply(lambda x: abs(x[ratio] - df_mean_before), axis=1)
    mytest[dif_1SD] = mytest.apply(lambda x: x[dif_mean] - df_1std_before, axis=1)

    mytest_noNeg = mytest[mytest[ratio] > 0]
    
    mytest2 = mytest_noNeg.sort_values(by = dif_1SD, ascending = False)
    #mytest2.head()

    ratios_to_reject = int(mytest_noNeg[ratio].count() * reject_percentage / 100)
    #print(ratios_to_reject)

    after_rejection = mytest2[ratios_to_reject:]
    
    

    df_mean_after = after_rejection[ratio].mean()
    df_1std_after = after_rejection[ratio].std()
    df_count_after = after_rejection[ratio].count()
   
    df_2se_perc_after = (2 * after_rejection[ratio].sem()) / df_mean_after * 100
    df_2sd_perc_after = (2 * after_rejection[ratio].std()) / df_mean_after * 100
    

    # print(df_mean_after)
    # print(df_1std_after)
    # print(df_2se_perc_after)

    results_dict = {}
    
    results_dict['avg_before'] = df_mean_before
    results_dict['1sd_before'] = df_1std_before
    results_dict['2se%_before'] = df_2se_perc_before
    results_dict['avg_after'] = df_mean_after
    results_dict['1sd_after'] = df_1std_after
    results_dict['2se%_after'] = df_2se_perc_after
    results_dict['2σ%_after'] = df_2sd_perc_after
    
    
    
    
    return results_dict

def statistics_ranktest2(calc_dict, reject_percentage = 20):
    '''Calculates statistics while Incorporating the ranked error minimization.'''
    calc_list = ['238_CPS', '232_CPS',
           '208_CPS', '207_CPS', '206_CPS', '204_CPS', '202_CPS', '206/238',
           '208/232', '207/206', '208/206', '206/204','208/204','207/204', '207/235'  ]
    mega_dict = {}

    for sheet in calc_dict:
        tester = calc_dict[sheet]
        stats_dict = {}
        
        
        ###Regression tests###
        
        regress_dict = regress_indiv(calc_dict, sheet, False, 1, 1)
        
        
        
        
        ###End Regression tests###
        
        for col in tester:

            if col in calc_list:
                #print(col)
                if '/' in col:
                    key_bf = col + '_before rejection'
                    key_af = col + '_after rejection'
                    key_bf_se = col + '_before rejection 2se%'
                    key_af_se = col + '_after rejection 2se%'
                    key_af_2sd = col + '_after rejection 2σ%'
                    
                    ranked_dict = ranked_minimization2(sheet, col, reject_percentage)
                    
                    stats_dict[key_bf] = ranked_dict['avg_before']
                    stats_dict[key_bf_se] = ranked_dict['2se%_before']
                    stats_dict[key_af] = ranked_dict['avg_after']
                    
                    stats_dict[key_af_se] = ranked_dict['2se%_after']
                    stats_dict[key_af_2sd] = ranked_dict['2σ%_after']
                    
                    ###Regression tests### 
                    if col == '206/238':
                        key_int = col + '_intercept'
                        stats_dict[key_int] = regress_dict['intercept']

                    
                    ###End Regression tests###
                    
                else:
                    key = col + '_mean'
                    df_mean = tester[col].mean()
                    stats_dict[key] = df_mean
                    #This is not 2SE%, should probably fix labels... KD 14 June 2021
                    df_precision = (2 * tester[col].std()) / df_mean * 100
                    stats_dict[col + '_2σ%'] = df_precision
            if 'OPZ' in col:
                 stats_dict[col + '_mean'] = tester[col].mean()
            if 'SNR' in col:
                 stats_dict[col + '_mean'] = tester[col].mean()
                    
                    
                    
                    
                    
                
        stats_dict['Time (s)'] = tester['Elapsed Time'].max()
        
        #new_string = sheet.replace('time series data', '')
        new_string = sheet.split('time')[0].rstrip()
        mega_dict[new_string] = stats_dict

    df_1 = pd.DataFrame(mega_dict)
    df_flip = pd.DataFrame.transpose(df_1)
   
    return df_flip


# + code_folding=[2, 41, 123]
### Functions for fractionation output .xlsx: frac_np2_timeseries(), frac_calc_CPS(), fractionation_output()

def frac_np2_timeseries(excel_file):
    ''' Excel input file is your baseline corrected time series export from Iolite for the NP-II.'''
    df = pd.read_excel(excel_file, sheet_name = None)
    keys = df.keys()
    header_row = 0
    new_dict = {}
    frac_dict = []
    for key in keys:
        if 'time' in key: #Kind of hard-coded right now, so if names get weird may need to change
            
            df_test = df[key]

            df_test.columns = df_test.iloc[header_row]
            df_test = df_test.drop(header_row)
            df_test = df_test.reset_index(drop=True)
            
            
            new_string = key.split('time')[0].rstrip()
            
            #Adding column with SAMPLE ID
            cols = list(df_test.columns.values)
            df_test['SAMPLE ID'] = new_string
            new_cols = ['SAMPLE ID'] + cols
            df_test = df_test[new_cols]
            
            new_dict[new_string] = df_test #test1_new
            frac_dict.append(new_string)

    combo = new_dict[frac_dict[0]]
    for idx in enumerate(frac_dict):
        #print(idx[0])
        if idx[0] == 0:
            continue
        else:
            combo = pd.concat([combo, new_dict[idx[1]]])
    #print(frac_dict)
    #print(combo)
    return new_dict, combo

def frac_calc_CPS(combo):
    columns = ['Absolute Time',
     'Elapsed Time',
     'm238_CPS',
     'm232_CPS',
     'm208_CPS',
     'm207_CPS',
     'm206_CPS',
     'm204_CPS',
     'm202_CPS']

    new_col = ['Absolute Time',
     'Elapsed Time',
     '238_CPS',
     '232_CPS',
     '208_CPS',
     '207_CPS',
     '206_CPS',
     '204_CPS',
     '202_CPS']

    cut_col = ['238_CPS',
     '232_CPS',
     '208_CPS',
     '207_CPS',
     '206_CPS',
     '204_CPS',
     '202_CPS']

#     calc_dict = {}
#     for key in np2_dict:
#         #print(key)
#         test_df1 = np2_dict[key]
    test_df1 = combo
    for col in columns:
            test_df2 = test_df1.apply(lambda x: x * 62500000 if 'CPS' in x.name else x)
            test_df2 = test_df2[['Absolute Time',
         'Elapsed Time',
             'm238_CPS',
             'm232_CPS',
             'm208_CPS',
             'm207_CPS',
             'm206_CPS',
             'm204_CPS',
             'm202_CPS',]]
    test_df2.columns = new_col
    test_df2 = test_df2[cut_col]
    result = pd.concat([test_df1, test_df2], axis=1)

     #Calculating OPZ
    result['OPZ_238'] = result.apply(lambda x: x['m238'] - x['m238_CPS'], axis=1)
    result['OPZ_232'] = result.apply(lambda x: x['m232'] - x['m232_CPS'], axis=1)
    result['OPZ_208'] = result.apply(lambda x: x['m208'] - x['m208_CPS'], axis=1)
    result['OPZ_207'] = result.apply(lambda x: x['m207'] - x['m207_CPS'], axis=1)
    result['OPZ_206'] = result.apply(lambda x: x['m206'] - x['m206_CPS'], axis=1)
    result['OPZ_204'] = result.apply(lambda x: x['m204'] - x['m204_CPS'], axis=1)
    result['OPZ_202'] = result.apply(lambda x: x['m202'] - x['m202_CPS'], axis=1)
    result['OPZ_208/206'] = result.apply(lambda x: x['OPZ_208'] / x['OPZ_206'], axis=1)
    result['OPZ_207/206'] = result.apply(lambda x: x['OPZ_207'] / x['OPZ_206'], axis=1)
    result['OPZ_206/204_Hg-corrected'] = result.apply(lambda x: x['OPZ_206'] / (x['OPZ_204'] - (x['OPZ_202'] * 6.87/29.86)) , axis=1)

    #Calculating Signal-to-Noise Ratios [V]/[V]
    result['SNR_238'] = result.apply(lambda x: x['m238_CPS']/ x['OPZ_238'], axis=1)
    result['SNR_232'] = result.apply(lambda x: x['m232_CPS']/ x['OPZ_232'], axis=1)
    result['SNR_208'] = result.apply(lambda x: x['m208_CPS']/ x['OPZ_208'], axis=1)
    result['SNR_207'] = result.apply(lambda x: x['m207_CPS']/ x['OPZ_207'], axis=1)
    result['SNR_206'] = result.apply(lambda x: x['m206_CPS']/ x['OPZ_206'], axis=1)
    result['SNR_204'] = result.apply(lambda x: x['m204_CPS']/ x['OPZ_204'], axis=1)
    result['SNR_202'] = result.apply(lambda x: x['m202_CPS']/ x['OPZ_202'], axis=1)

    #Calculating Ratios
    result['206/238'] = result.apply(lambda x: x['206_CPS']/x['238_CPS'], axis=1)
    result['208/232'] = result.apply(lambda x: x['208_CPS']/x['232_CPS'], axis=1)
    result['207/206'] = result.apply(lambda x: x['207_CPS']/x['206_CPS'], axis=1)
    result['208/206'] = result.apply(lambda x: x['208_CPS']/x['206_CPS'], axis=1)
    result['206/204'] = result.apply(lambda x: x['206_CPS']/x['204_CPS'], axis=1)
    result['208/204'] = result.apply(lambda x: x['208_CPS']/x['204_CPS'], axis=1)
    result['207/204'] = result.apply(lambda x: x['207_CPS']/x['204_CPS'], axis=1)

    
    return result

def fractionation_output(excel_file):
    frac_tester, frac_combo = frac_np2_timeseries(excel_file)
    combo_tester = frac_calc_CPS(frac_combo)
    
    new_filename = str(excel_file.split('.')[0]) + '_fractionation.xlsx'   
    with pd.ExcelWriter(new_filename) as writer:
        combo_tester.to_excel(writer, sheet_name = 'Fractionation_ALL', index = False)
       
        worksheet = writer.sheets['Fractionation_ALL']  # pull worksheet object
        worksheet.set_column(0, 0, 24)
        worksheet.freeze_panes(1,1)   
        
    #return combo_tester
       


# + code_folding=[2, 22]
### Functions for exporting the time-series processed data to Excel: files_process_toEXCEL(), file_process_combine()

def files_process_toEXCEL(calc_dict, excel_name, order):
    ''' Takes the dataframe from calcCPS() and outputs the time-series processed .xlsx with group statistics in group order and chronologic order. 
    If chronologic order is not updated, it might break...'''
    
    with pd.ExcelWriter(excel_name) as writer:
        for sheet in calc_dict:
            calc_dict[sheet].to_excel(writer, sheet_name = sheet, index = False)
        
        statistics_NP2(calc_dict).to_excel(writer, sheet_name = 'Statistics', index = True)
        
        chronologic(statistics_NP2(calc_dict), order).to_excel(writer, sheet_name = 'Chrono_Statistics', index = True)
        
        worksheet = writer.sheets['Statistics']  # pull worksheet object
        worksheet.set_column(0, 0, 24)
        worksheet.freeze_panes(1,1)
        
        worksheet = writer.sheets['Chrono_Statistics']  # pull worksheet object
        worksheet.set_column(0, 0, 24)
        worksheet.freeze_panes(1,1)
        
def file_process_combine(filename):
    ''' CURRENTLY NOT IN USE. This function will take the filename and output the time-series data processed in .xlsx format. Useful for quick data check.'''
    calc_dict = calc_CPS(read_np2_timeseries(filename))
    new_filename = str(filename.split('.')[0]) + '_processed.xlsx'
    files_process_toEXCEL(calc_dict, new_filename)


# + code_folding=[1, 85]
### Functons for graphing and report generation: U_Pb_plots(), U_Pb_report()
def U_Pb_plots(calc_dict, sample, choice = True):
    ''' Makes plots for CPS and isotope ratios plotted against time. The input is the dataframe from calc_CPS() and a sample.'''
    key_list = ['238_CPS', '232_CPS',
       '208_CPS', '207_CPS', '206_CPS', '204_CPS', '202_CPS', '206/238',
       '208/232', '207/206', '208/206', '206/204']
    
    zet = calc_dict[sample]
    new_string = sample.split('time')[0].rstrip()
    y_list = []
    for key in key_list:
        y_list.append(zet[key])
    
    x = zet['Elapsed Time']
    
    #mod 12 September 2021
    y_238over232 = zet['238_CPS'] / zet['232_CPS']
    new_list = ['238_CPS', '232_CPS',
       '208_CPS', '207_CPS', '206_CPS', '204_CPS', '238/232', '206/238',
       '208/232', '207/206', '208/206', '206/204']
    
    y_list[6] = y_238over232
    
    
    fig, axs = plt.subplots(4, 3, sharex = True, figsize = (12, 12))
    fig.suptitle(new_string, fontsize=24)
    
    ax_list = [
        axs[0, 0], 
        axs[0, 1],   
        axs[0, 2], 
        axs[1, 0], 
        axs[1, 1],
        axs[1, 2],
        axs[2, 0], 
        axs[2, 1], 
        axs[2, 2],
        axs[3, 0], 
        axs[3, 1], 
        axs[3, 2]   
        ]

    axs[0, 0].plot(x, y_list[0])
    axs[0, 1].plot(x, y_list[1])
    axs[0, 2].plot(x, y_list[2])
    axs[1, 0].plot(x, y_list[3])
    axs[1, 1].plot(x, y_list[4])
    axs[1, 2].plot(x, y_list[5])
    axs[2, 0].plot(x,  y_list[6])
    axs[2, 1].plot(x, y_list[7])
    axs[2, 2].plot(x, y_list[8])
    axs[3, 0].plot(x, y_list[9])
    axs[3, 0].set(xlabel = 'Time (s)')
    axs[3, 1].plot(x, y_list[10])
    axs[3, 1].set(xlabel = 'Time (s)')
    axs[3, 2].plot(x, y_list[11])
    axs[3, 2].set(xlabel = 'Time (s)')
    for idx in range(len(ax_list)):
        ax_list[idx].ticklabel_format(axis='y', style='sci', scilimits=(0,0))
        ax_list[idx].set_title(new_list[idx])
        y_mean = [np.mean(y_list[idx])]*len(x)
        # Plot the average line
        mean_line = ax_list[idx].plot(x,y_mean, label='Mean', linestyle='--', color = "black")
        # Make a legend
        legend = ax_list[idx].legend(loc='upper right')
    
    MYDIR = ("Figures")
    CHECK_FOLDER = os.path.isdir(MYDIR)

    # If folder doesn't exist, then create it.
    if not CHECK_FOLDER:
        os.makedirs(MYDIR)
        #print("created folder : ", MYDIR)
    
    #new_string = sample.replace('time series data', '').rstrip()
    
    filename = os.path.join(MYDIR, new_string + '.pdf')
    plt.savefig(filename)
    print('Plot for ', new_string, " is complete.")
    
    if choice == False:
        plt.close()
    #else:
        #plt.close()
        
def U_Pb_report(calc_dict, intro_filename, intro = False, output_name = 'U-Pb_output.pdf'):
    """ Makes the U-Pb plots into a super long report. The input is the dataframe from calcCPS(). Could definitely be revised..."""
    MYDIR = ("Figures")
    mergedObject = PdfFileMerger()
    
    if intro:
        mergedObject.append(PdfFileReader(intro_filename, 'rb'))
        print(f'Succesfully incorporated {intro_filename} into PDF.')
    pd.set_option('precision', 2)
    stats = statistics_NP2(calc_dict)
    stat_dict = {}
    stat_dict['stat1'] = stats.iloc[:, 14:]
    stat_dict['stat2'] = stats.iloc[:, :8]
    stat_dict['stat3'] = stats.iloc[:, 8:14]
    html_list = []

    for key in stat_dict:
        name = key + ".pdf"
        stats_html = stat_dict[key].to_html()
        pdfkit.from_string(stats_html, name)
        mergedObject.append(PdfFileReader(name, 'rb'))

    file_list = []
    keys = calc_dict.keys()

    for key in keys:
        #print(key)
        U_Pb_plots(calc_dict, key, False)
        new_string = key.split('time')[0].rstrip()
        filename = os.path.join(MYDIR, new_string + '.pdf')

        mergedObject.append(PdfFileReader(filename, 'rb'))

    if '.pdf' in output_name:
        pass
    else:
        output_name = output_name + '.pdf'
    
    #output_name = "U-Pb_output.pdf"  
    mergedObject.write(output_name)

    print(f'PDF file named: {output_name} is complete.')        
        


# + code_folding=[2, 15, 67, 87, 101, 112]
### Functions for calculating excess variance (and weighted mean and MSWD, if needed): excV_prep(), MSWD_stats(), excess_var(), all_excVariance(), group_excVariance(), stat_to_excV()

def excV_prep(stat_test1):
    '''CURRENTLY NOT IN USE. Reduces and changes order of ratios and uncertainties required to calculate excess variance. Might be sueful later for error propogation or Excel export.'''
    
    
    new_key_order = ['206/238_after rejection','206/238_after rejection 2se%', '208/232_after rejection', '208/232_after rejection 2se%',
        '207/206_after rejection','207/206_after rejection 2se%', '208/206_after rejection', '208/206_after rejection 2se%',
        '206/204_after rejection','206/204_after rejection 2se%', '208/204_after rejection', '208/204_after rejection 2se%',
       '207/204_after rejection','207/204_after rejection 2se%', 
            ]
    stat_test2 = stat_test1[new_key_order]
    
    return stat_test2

def MSWD_stats(ratio_input, excess_variance, mked_df, choice = 'MSWD'):
    '''Given an input ratio, excess variance, and dataframe with the correct columns, this function will output the MSWD for that ratio. If choice is not "MSWD" this function will output all steps of the MSWD calculation.'''
    ratio = ratio_input +'_after rejection'

    mked_stat_dict = {}

    #print(ratio)
    error_string = str(ratio + ' 2se%')

    new_string1 = ratio + ' abs err'

    mod_string1 = ratio + ' mod err'
    #print(error_string)
    new_string2 = ratio + ' x/(σ^2)'
    new_string3 = ratio + ' 1/(σ^2)'

    mked_df[new_string1] = mked_df.apply(lambda x: (x[error_string] / 100 ) * (x[ratio]), axis=1)
    mked_df[mod_string1] = mked_df.apply(lambda x: (math.sqrt( x[error_string]**2 + excess_variance**2 ) / 100 ) * (x[ratio]), axis=1)

    mked_df[new_string2] = mked_df.apply(lambda x: x[ratio] / ((x[mod_string1] / 2)**2), axis=1)
    mked_df[new_string3] = mked_df.apply(lambda x: (x[new_string2] / x[ratio]), axis=1)

    wtd_mean_string = str(ratio.split()[0]) + '_wtdMean'
    one_sig2_string = str(ratio.split()[0]) + '_sig^2'
    one_sig_string = str(ratio.split()[0]) + '_sig'

    new_string2 = ratio + ' x/(σ^2)'
    new_string3 = ratio + ' 1/(σ^2)'


    mked_stat_dict[wtd_mean_string] = mked_df[new_string2].sum() / (mked_df[new_string3].sum())
    mked_stat_dict[one_sig2_string] = 1 / (mked_df[new_string3].sum())

    mked_stat_dict[one_sig_string] = math.sqrt(mked_stat_dict[one_sig2_string])

    z_string = str(ratio.split()[0]) + '_Z'
    wtd_mean_string = str(ratio.split()[0]) + '_wtdMean'
    new_string3 = ratio + ' 1/(σ^2)'
    S_string = ratio + '_S'
    MSWD_string = ratio + '_MSWD'


    mked_df[z_string] = mked_df.apply(lambda x: ( ((x[ratio]  - mked_stat_dict[wtd_mean_string])**2) * x[new_string3]), axis=1)
    mked_stat_dict[S_string]=  mked_df[z_string].sum()
    mked_stat_dict[ratio + "_MSWD"] = mked_stat_dict[S_string] / ((mked_df[z_string].count()) - 1)

    if choice == 'MSWD':
        return mked_stat_dict[MSWD_string]  

    else:
        return MKED_stat_dict
    
def excess_var(ratio, mked_df):
    ''' Given a specific ratio and the primary standard dataframe with correct columns, 
    this function will calculate the excess variance (rounded to less than 6 significant figures) needed to correct for overdispersion.'''
    eVar = 0
    #pick = '207/206'
    MSWD = MSWD_stats(ratio, eVar, mked_df)
    print('Initial MSWD=', MSWD)
    
    for i in range(5):
        while MSWD > 1.000:
            eVar = eVar + 1/(10**i)
            MSWD = MSWD_stats(ratio, eVar, mked_df)

        eVar = eVar - 1/(10**i)
        MSWD = MSWD_stats(ratio, eVar, mked_df)

    print('excess variance=', round_sig(eVar, 6), 'MSWD=', round_sig(MSWD, 5), '\n')
    
    return  round_sig(eVar, 6) 

def all_excVariance(mked_df):
    ''' Uses the excess_var() function to calculate the excess variance for all relevant ratios.'''
    ratios = ['206/238','208/232','207/206','208/206',
            '206/204','207/204','208/204'   ]

    eVar_dict = {}

    for rat in ratios:
        print(rat)
        new_string = str(rat + '_excess variance')
        eVar_dict[new_string] = excess_var(rat, mked_df)
    
    return eVar_dict

def group_excVariance(stat_test2, group_name):
    ''' Takes an input from statistics_ranktest2() and a chosen primary standard and calculates the excess variance for each isotope ratio.'''
    AB_err_tester = stat_test2.copy()

    y_list, labels = comb_groups([group_name],AB_err_tester)
    mked_df = pd.DataFrame(y_list)
    print('Calculating excess variance for: ', group_name, '\n')
    return all_excVariance(mked_df)   

#mked_excV = group_excVariance(stat_test2, 'MKED-1')

def stat_to_excV(tester, primary_std = 'MKED-1', reject_percentage = 0):
    ''' Function to take the resulting df from calc_CPS() and calculate the excess variance for a chosen primary standard and error minimization.'''
    stattie = statistics_ranktest2(tester, reject_percentage)
    new_excV = group_excVariance(stattie, primary_std)
    
    return new_excV



# + code_folding=[]
### Functions for doing correction math and error minimization: stat_correcter2(), stat_rank_and_correct()

excv_dict = {}
def stat_correcter2(stats_df, std = 'glass', excV_dict = excv_dict):
    ''' Input dataframe that was output from statistics_ranktest2
    function. Performs correction using specified primary standard. Slightly hard-coded for published standard values.'''
    headers_to_keep = [
        '206/238_after rejection',
        '206/238_after rejection 2se%',
        '206/238_intercept',       #Regression test
        '208/232_after rejection',
        '208/232_after rejection 2se%',
        '207/206_after rejection',
        '207/206_after rejection 2se%',
        '208/206_after rejection',
        '208/206_after rejection 2se%',
        '206/204_after rejection',
        '206/204_after rejection 2se%',
        '208/204_after rejection',
        '208/204_after rejection 2se%',
        '207/204_after rejection',
        '207/204_after rejection 2se%',
        '207/235_after rejection',
        '207/235_after rejection 2se%',
        
    ]

    headers_to_math = [
        '206/238_after rejection',
        '206/238_intercept',       #Regression test
        '208/232_after rejection',  
        '207/206_after rejection',    
        '208/206_after rejection',    
        '206/204_after rejection',    
        '208/204_after rejection',   
        '207/204_after rejection',
        '207/235_after rejection',
    ]
    
### Correction Math
    
    #Calling relevant standard information
    
    if std == 'glass':
        standard = 'SRM NIST 610'   
    elif std == '612':
        standard = 'SRM NIST 612'        
    elif std == 'BHVO':
        standard = 'BHVO-2G' 
    elif std == 'BLR-1':
        standard = 'BLR-1'
    else:
        standard = 'MKED-1'
       
    published_ratios = published_full(standard) 
    
    short_tester = stats_df[headers_to_keep] #Reduces columns. Might break because 2sigma has been excluded.
    result = short_tester.copy(deep=True)
    
    #Getting list of group names and corresponding index positons
    
    start, end, group_names = group_samples(result)

    #print('start:', start)
    #print('end:', end)
    #print('group names:', group_names)

 
    stats_dict = {}
    corrected_stats_dict = {}
    for idx in range(len(group_names)):
        group_dict = {}
        #print(idx)
        group_df = short_tester.iloc[start[idx]:end[idx]]
        #print(group_df.index.values.tolist())
        for col in headers_to_math:
            
            ###Regression Test###
            if col == '206/238_intercept':
                df_mean = group_df[col].mean()
                group_dict[col + '_mean'] = df_mean
            ###End Regression Test###
            else:
                #print(group_df[col])
                name = col.split('_')[0]
                df_mean = group_df[col].mean()
                #print(df_mean)
                group_dict[name + '_mean'] = df_mean
                #External error
                df_2sd_perc = (2 * group_df[col].std())/df_mean * 100

                group_dict[name + '_2sd%'] = df_2sd_perc

                df_2se_perc = (2 * group_df[col].sem())/df_mean * 100

                group_dict[name + '_2se%'] = df_2se_perc        #Not sure if the group 2se% is really a relevant statistic...

        stats_dict[group_names[idx]] = group_dict
    
    grouper = pd.DataFrame(stats_dict)
       
    primary_std = grouper[standard]
        
    
        #Applying correction Hard coding the primary standard
    result['206/238 corrected'] = short_tester.apply(lambda x: x['206/238_after rejection']/ (primary_std['206/238_mean'] / published_ratios['206/238']), axis=1)
    result['208/232 corrected'] = short_tester.apply(lambda x: x['208/232_after rejection']/ (primary_std['208/232_mean'] / published_ratios['208/232']), axis=1)
    result['207/206 corrected'] = short_tester.apply(lambda x: x['207/206_after rejection']/ (primary_std['207/206_mean'] / published_ratios['207/206']), axis=1)
    result['208/206 corrected'] = short_tester.apply(lambda x: x['208/206_after rejection']/(primary_std['208/206_mean'] / published_ratios['208/206']), axis=1)
    result['206/204 corrected'] = short_tester.apply(lambda x: x['206/204_after rejection']/(primary_std['206/204_mean'] / published_ratios['206/204']), axis=1)
    result['207/204 corrected'] = short_tester.apply(lambda x: x['207/204_after rejection']/(primary_std['207/204_mean'] / published_ratios['207/204']), axis=1)
    result['208/204 corrected'] = short_tester.apply(lambda x: x['208/204_after rejection']/(primary_std['208/204_mean'] / published_ratios['208/204']), axis=1)
    result['207/235 corrected'] = short_tester.apply(lambda x: x['207/235_after rejection']/(primary_std['207/235_mean'] / published_ratios['207/235']), axis=1)
    
    result['206/238_int corrected'] = short_tester.apply(lambda x: x['206/238_intercept']/ (primary_std['206/238_intercept_mean'] / published_ratios['206/238']), axis=1)
    
#     #Making 2se for plotting in EXCEL
    
    result['206/238_after rejection 2se'] = short_tester.apply(lambda x: x['206/238_after rejection'] * x['206/238_after rejection 2se%'] / 100, axis=1)
    result['208/232_after rejection 2se'] = short_tester.apply(lambda x: x['208/232_after rejection'] * x['208/232_after rejection 2se%'] / 100, axis=1)
    result['207/206_after rejection 2se'] = short_tester.apply(lambda x: x['207/206_after rejection'] * x['207/206_after rejection 2se%'] / 100, axis=1)
    result['208/206_after rejection 2se'] = short_tester.apply(lambda x: x['208/206_after rejection'] * x['208/206_after rejection 2se%'] / 100, axis=1)
    result['206/204_after rejection 2se'] = short_tester.apply(lambda x: x['206/204_after rejection'] * x['206/204_after rejection 2se%'] / 100, axis=1)
    result['207/204_after rejection 2se'] = short_tester.apply(lambda x: x['207/204_after rejection'] * x['207/204_after rejection 2se%'] / 100, axis=1)
    result['208/204_after rejection 2se'] = short_tester.apply(lambda x: x['208/204_after rejection'] * x['208/204_after rejection 2se%'] / 100, axis=1)
    result['207/235_after rejection 2se'] = short_tester.apply(lambda x: x['207/235_after rejection'] * x['207/235_after rejection 2se%'] / 100, axis=1)  
    
    result_headers = [
        '206/238 corrected',
        '208/232 corrected',
        '207/206 corrected',
        '208/206 corrected',
        '206/204 corrected',
        '207/204 corrected',
        '208/204 corrected',
        '207/235 corrected',
    ]      
    
    for idx in range(len(group_names)):
        group_dict = {}
        #print(idx)
        group_df = result.iloc[start[idx]:end[idx]]
        #print(group_df.index.values.tolist())
        #print(group_df)
        for col in result_headers:
            #print(group_df[col])
            name = col.split('_')[0]
            df_mean = group_df[col].mean()
            #print(df_mean)
            group_dict[name + '_mean'] = df_mean
            df_2sd_perc = (2 * group_df[col].std())/df_mean * 100
            group_dict[name + '_2sd%'] = df_2sd_perc 
            df_2sd = 2 * group_df[col].std()
            group_dict[name + '_2sd'] = df_2sd
            
            df_2se_perc = (2 * group_df[col].sem())/df_mean * 100
            group_dict[name + '_2se%'] = df_2se_perc
            df_2se = 2 * group_df[col].sem()
            group_dict[name + '_2se'] = df_2se
        
        corrected_stats_dict[group_names[idx]] = group_dict
   
    corrected_grouper = pd.DataFrame(corrected_stats_dict)
     
    grouper_flip = pd.DataFrame.transpose(grouper)
    corrected_flip = pd.DataFrame.transpose(corrected_grouper)
    grouper_list = [grouper_flip, corrected_flip]
    
    grouper_comb = pd.concat(grouper_list, axis=1)
    
        
    #print('Primary Standard= ', standard)
    
    #Calculate excess variance
    #excV_dict = group_excVariance(result, standard)
    
    #print(excV_dict)
    ####
    
    external_err_std = grouper_comb.loc[standard]  
     
    result['206/238 corrected BB error%'] = short_tester.apply(lambda x: ((x['206/238_after rejection 2se%'])**2 + (external_err_std['206/238 corrected_2sd%'])**2)**0.5, axis=1)
    result['208/232 corrected BB error%'] = short_tester.apply(lambda x: ((x['208/232_after rejection 2se%'])**2 + (external_err_std['208/232 corrected_2sd%'])**2)**0.5, axis=1)
    result['207/206 corrected BB error%'] = short_tester.apply(lambda x: ((x['207/206_after rejection 2se%'])**2 + (external_err_std['207/206 corrected_2sd%'])**2)**0.5, axis=1)
    result['208/206 corrected BB error%'] = short_tester.apply(lambda x: ((x['208/206_after rejection 2se%'])**2 + (external_err_std['208/206 corrected_2sd%'])**2)**0.5, axis=1)
    result['206/204 corrected BB error%'] = short_tester.apply(lambda x: ((x['206/204_after rejection 2se%'])**2 + (external_err_std['206/204 corrected_2sd%'])**2)**0.5, axis=1)
    result['207/204 corrected BB error%'] = short_tester.apply(lambda x: ((x['207/204_after rejection 2se%'])**2 + (external_err_std['207/204 corrected_2sd%'])**2)**0.5, axis=1)
    result['208/204 corrected BB error%'] = short_tester.apply(lambda x: ((x['208/204_after rejection 2se%'])**2 + (external_err_std['208/204 corrected_2sd%'])**2)**0.5, axis=1)
    result['207/235 corrected BB error%'] = short_tester.apply(lambda x: ((x['207/235_after rejection 2se%'])**2 + (external_err_std['207/235 corrected_2sd%'])**2)**0.5, axis=1)
       
    result['206/238 corrected BB error'] = result.apply(lambda x: (x['206/238 corrected BB error%'])* x['206/238 corrected']/100, axis=1)
    result['208/232 corrected BB error'] = result.apply(lambda x: (x['208/232 corrected BB error%'])* x['208/232 corrected']/100, axis=1)
    result['207/206 corrected BB error'] = result.apply(lambda x: (x['207/206 corrected BB error%'])* x['207/206 corrected']/100, axis=1) 
    result['208/206 corrected BB error'] = result.apply(lambda x: (x['208/206 corrected BB error%'])* x['208/206 corrected']/100, axis=1)
    result['206/204 corrected BB error'] = result.apply(lambda x: (x['206/204 corrected BB error%'])* x['206/204 corrected']/100, axis=1)   
    result['207/204 corrected BB error'] = result.apply(lambda x: (x['207/204 corrected BB error%'])* x['207/204 corrected']/100, axis=1)
    result['208/204 corrected BB error'] = result.apply(lambda x: (x['208/204 corrected BB error%'])* x['208/204 corrected']/100, axis=1)
    result['207/235 corrected BB error'] = result.apply(lambda x: (x['207/235 corrected BB error%'])* x['207/235 corrected']/100, axis=1)
    
    #Propagating the measurement error (internal precision; 2se%) with 
    #the primary standard excess variance(% to add to the primary standard internal precision to correct for overdispersion). 
    #This is following the direction of Horstwood et al., (2016) on how to propogate the population uncertainty (Sx).
    
    
    
    result['206/238 corrected Sx%'] = result.apply(lambda x: ((x['206/238_after rejection 2se%'])**2 + (excV_dict['206/238_excess variance'])**2)**0.5, axis=1)
    result['208/232 corrected Sx%'] = result.apply(lambda x: ((x['208/232_after rejection 2se%'])**2 + (excV_dict['208/232_excess variance'])**2)**0.5, axis=1)
    result['207/206 corrected Sx%'] = result.apply(lambda x: ((x['207/206_after rejection 2se%'])**2 + (excV_dict['207/206_excess variance'])**2)**0.5, axis=1)
    result['208/206 corrected Sx%'] = result.apply(lambda x: ((x['208/206_after rejection 2se%'])**2 + (excV_dict['208/206_excess variance'])**2)**0.5, axis=1)
    result['206/204 corrected Sx%'] = result.apply(lambda x: ((x['206/204_after rejection 2se%'])**2 + (excV_dict['206/204_excess variance'])**2)**0.5, axis=1)
    result['207/204 corrected Sx%'] = result.apply(lambda x: ((x['207/204_after rejection 2se%'])**2 + (excV_dict['207/204_excess variance'])**2)**0.5, axis=1)
    result['208/204 corrected Sx%'] = result.apply(lambda x: ((x['208/204_after rejection 2se%'])**2 + (excV_dict['208/204_excess variance'])**2)**0.5, axis=1)
    
    result['207/235 corrected Sx%'] = result.apply(lambda x: ((x['206/238 corrected Sx%'])**2 + (x['207/206 corrected Sx%'])**2)**0.5, axis=1) #Just quadratic addition of 206/238 error and 207/206 error
      
     
    new_headers = [
        '206/238 corrected BB error%',
        '208/232 corrected BB error%',
        '207/206 corrected BB error%',
        '208/206 corrected BB error%',
        '206/204 corrected BB error%',
        '207/204 corrected BB error%',
        '208/204 corrected BB error%',
        '207/235 corrected BB error%',
        '206/238 corrected BB error',
        '208/232 corrected BB error',
        '207/206 corrected BB error',
        '208/206 corrected BB error',
        '206/204 corrected BB error',
        '207/204 corrected BB error',
        '208/204 corrected BB error',
        '207/235 corrected BB error',
        '206/238 corrected Sx%',
        '208/232 corrected Sx%',
        '207/206 corrected Sx%',
        '208/206 corrected Sx%',
        '206/204 corrected Sx%',
        '207/204 corrected Sx%',
        '208/204 corrected Sx%',
        '207/235 corrected Sx%',
        
    ]
    
    new_stats_dict = {}
    for idx in range(len(group_names)):
        group_dict = {}
        #print(idx)
        group_df = result.iloc[start[idx]:end[idx]]
        #print(group_df)
        for col in new_headers:
            df_mean = group_df[col].mean()
            group_dict[col + '_mean'] = df_mean
        
        new_stats_dict[group_names[idx]] = group_dict 
            
    extra_grouper = pd.DataFrame(new_stats_dict)
     
    corrected_flip2 = pd.DataFrame.transpose(extra_grouper)
    grouper_list2 = [grouper_comb, corrected_flip2]
    
    grouper_comb_final = pd.concat(grouper_list2, axis=1)
    #print('GROUPER',grouper_comb_final.keys())
    #print('RESULT', result.keys())
    ############

    
    grouper_order = ['206/238_mean', '206/238_2sd%','206/238_2se%', '208/232_mean', '208/232_2sd%','208/232_2se%',
       '207/206_mean', '207/206_2sd%','207/206_2se%', '208/206_mean', '208/206_2sd%', '208/206_2se%',
       '206/204_mean', '206/204_2sd%','206/204_2se%', '208/204_mean', '208/204_2sd%','208/204_2se%',
       '207/204_mean', '207/204_2sd%', '207/204_2se%','207/235_mean', '207/235_2sd%', '207/235_2se%',
                
        '206/238 corrected_mean','206/238 corrected_2sd','206/238 corrected_2sd%','206/238 corrected_2se%',
       '206/238 corrected_2se', '206/238 corrected BB error_mean', '206/238 corrected BB error%_mean', '206/238 corrected Sx%_mean',
                     
       '208/232 corrected_mean','208/232 corrected_2sd', '208/232 corrected_2sd%', '208/232 corrected_2se%', 
        '208/232 corrected_2se','208/232 corrected BB error_mean', '208/232 corrected BB error%_mean', '208/232 corrected Sx%_mean',
                     
       '207/206 corrected_mean', '207/206 corrected_2sd', '207/206 corrected_2sd%','207/206 corrected_2se%',
       '207/206 corrected_2se', '207/206 corrected BB error_mean','207/206 corrected BB error%_mean', '207/206 corrected Sx%_mean',
                     
        '208/206 corrected_mean', '208/206 corrected_2sd','208/206 corrected_2sd%', '208/206 corrected_2se%', 
        '208/206 corrected_2se','208/206 corrected BB error_mean', '208/206 corrected BB error%_mean', '208/206 corrected Sx%_mean',
                     
        '206/204 corrected_mean','206/204 corrected_2sd','206/204 corrected_2sd%', '206/204 corrected_2se%',
       '206/204 corrected_2se', '206/204 corrected BB error_mean', '206/204 corrected BB error%_mean', '206/204 corrected Sx%_mean',
                     
      '207/204 corrected_mean','207/204 corrected_2sd', '207/204 corrected_2sd%', '207/204 corrected_2se%', 
        '207/204 corrected_2se', '207/204 corrected BB error_mean', '207/204 corrected BB error%_mean', '207/204 corrected Sx%_mean',
                     
      '208/204 corrected_mean', '208/204 corrected_2sd','208/204 corrected_2sd%','208/204 corrected_2se%',
       '208/204 corrected_2se', '208/204 corrected BB error_mean', '208/204 corrected BB error%_mean', '208/204 corrected Sx%_mean',
                     
       '207/235 corrected_mean', '207/235 corrected_2sd', '207/235 corrected_2sd%','207/235 corrected_2se%',
       '207/235 corrected_2se', '207/235 corrected BB error_mean','207/235 corrected BB error%_mean', '207/235 corrected Sx%_mean',
            ]
    
    result_order = ['206/238_after rejection', '206/238_after rejection 2se%','206/238_after rejection 2se',
                   '208/232_after rejection','208/232_after rejection 2se%','208/232_after rejection 2se',
                   '207/206_after rejection',  '207/206_after rejection 2se%','207/206_after rejection 2se',
                   '208/206_after rejection', '208/206_after rejection 2se%','208/206_after rejection 2se',
                    '206/204_after rejection', '206/204_after rejection 2se%','206/204_after rejection 2se',
                   '208/204_after rejection', '208/204_after rejection 2se%','208/204_after rejection 2se',
                   '207/204_after rejection',  '207/204_after rejection 2se%','207/204_after rejection 2se',
                   '207/235_after rejection','207/235_after rejection 2se%','207/235_after rejection 2se',
                    
                   '206/238 corrected', '206/238_int corrected', '206/238 corrected BB error', '206/238 corrected BB error%', '206/238 corrected Sx%',
                    '208/232 corrected', '208/232 corrected BB error','208/232 corrected BB error%', '208/232 corrected Sx%',
                    '207/206 corrected', '207/206 corrected BB error', '207/206 corrected BB error%', '207/206 corrected Sx%',
                    '208/206 corrected', '208/206 corrected BB error', '208/206 corrected BB error%', '208/206 corrected Sx%',
                    '206/204 corrected', '206/204 corrected BB error', '206/204 corrected BB error%', '206/204 corrected Sx%',
                    '207/204 corrected', '207/204 corrected BB error', '207/204 corrected BB error%', '207/204 corrected Sx%',
                   '208/204 corrected',  '208/204 corrected BB error', '208/204 corrected BB error%', '208/204 corrected Sx%',
                    '207/235 corrected', '207/235 corrected BB error', '207/235 corrected BB error%', '207/235 corrected Sx%',
                     ]
        
        
    result_plotter = ['206/238 corrected',
                      '206/238_int corrected', '206/238 corrected Sx%','206/238 corrected BB error%', 
                      '206/238_after rejection 2se%', 
       '207/206 corrected', '207/206 corrected Sx%','207/206 corrected BB error%', 
                      '207/206_after rejection 2se%',  
        '208/232 corrected', '208/232 corrected Sx%','208/232 corrected BB error%',  
                      '208/232_after rejection 2se%', 
        '208/206 corrected',  '208/206 corrected Sx%', '208/206 corrected BB error%', 
                      '208/206_after rejection 2se%',
        '206/204 corrected', '206/204 corrected Sx%', '206/204 corrected BB error%', 
                      '206/204_after rejection 2se%', 
        '207/204 corrected', '207/204 corrected Sx%', '207/204 corrected BB error%', 
                      '207/204_after rejection 2se%', 
       '208/204 corrected', '208/204 corrected Sx%',  '208/204 corrected BB error%', 
                      '208/204_after rejection 2se%', 
        '207/235 corrected', '207/235 corrected Sx%', '207/235 corrected BB error%', 
                      '207/235_after rejection 2se%', 
         ]
          
    grouper_plotter = ['206/238 corrected_mean','206/238 corrected_2sd','206/238 corrected BB error_mean', '206/238 corrected Sx%_mean',
       '208/232 corrected_mean','208/232 corrected_2sd','208/232 corrected BB error_mean', '208/232 corrected Sx%_mean',
       '207/206 corrected_mean', '207/206 corrected_2sd', '207/206 corrected BB error_mean', '207/206 corrected Sx%_mean',
        '208/206 corrected_mean', '208/206 corrected_2sd','208/206 corrected BB error_mean', '208/206 corrected Sx%_mean',
        '206/204 corrected_mean','206/204 corrected_2sd', '206/204 corrected BB error_mean', '206/204 corrected Sx%_mean',
      '207/204 corrected_mean','207/204 corrected_2sd', '207/204 corrected BB error_mean', '207/204 corrected Sx%_mean',
      '208/204 corrected_mean', '208/204 corrected_2sd', '208/204 corrected BB error_mean', '208/204 corrected Sx%_mean',
     '207/235 corrected_mean','207/235 corrected_2sd', '207/235 corrected BB error_mean', '207/235 corrected Sx%_mean',                 
        ]
    
   
    result = result[result_order]
    grouper_comb_final = grouper_comb_final[grouper_order]
                                                        
    result_plot = result[result_plotter]
    grouper_plot = grouper_comb_final[grouper_plotter]                                                    
                                                                                                           
    
    return result, grouper_comb_final, result_plot, grouper_plot

excv_dict = {}
def stat_rank_and_correct(calc_dict, rank_perc = 0, std = 'glass', excv_dict = excv_dict):
    ''' Function for automating minimization and correction of isotope ratio data.'''
    stat_tester = statistics_ranktest2(calc_dict, rank_perc)
    
    
    print('')
    #print('Error minimization rejection percentage= ', rank_perc)
    result, grouper_comb_final, result_plot, grouper_plot = stat_correcter2(stat_tester, std, excv_dict)
    
    #print('Debug')
    stat1_order = ['238_CPS_mean', '238_CPS_2σ%', '232_CPS_mean', '232_CPS_2σ%',
       '208_CPS_mean', '208_CPS_2σ%', '207_CPS_mean', '207_CPS_2σ%',
       '206_CPS_mean', '206_CPS_2σ%', '204_CPS_mean', '204_CPS_2σ%',
       '202_CPS_mean', '202_CPS_2σ%', 'OPZ_238_mean', 'OPZ_232_mean',
       'OPZ_208_mean', 'OPZ_207_mean', 'OPZ_206_mean', 'OPZ_204_mean',
       'OPZ_202_mean','OPZ_208/206_mean','OPZ_207/206_mean','OPZ_206/204_Hg-corrected_mean',
        'SNR_238_mean', 'SNR_232_mean',
       'SNR_208_mean', 'SNR_207_mean', 'SNR_206_mean', 'SNR_204_mean',
       'SNR_202_mean', '206/238_before rejection',
       '206/238_before rejection 2se%',  '208/232_before rejection',
       '208/232_before rejection 2se%',  '207/206_before rejection',
       '207/206_before rejection 2se%',  '208/206_before rejection',
       '208/206_before rejection 2se%',  '206/204_before rejection',
       '206/204_before rejection 2se%',  '208/204_before rejection',
       '208/204_before rejection 2se%',  '207/204_before rejection',
       '207/204_before rejection 2se%',  'Time (s)']
    stat1 = stat_tester[stat1_order]
    
    stat2 = pd.concat([stat1, result], axis=1)
    stat_new = stat2
    group_new = grouper_comb_final
    new_grouper_plot = grouper_plot
    
    stat2_order = ['Time (s)','238_CPS_mean', '238_CPS_2σ%', '232_CPS_mean', '232_CPS_2σ%',
    '208_CPS_mean', '208_CPS_2σ%', '207_CPS_mean', '207_CPS_2σ%',
       '206_CPS_mean', '206_CPS_2σ%', '204_CPS_mean', '204_CPS_2σ%',
       '202_CPS_mean', '202_CPS_2σ%', 'OPZ_238_mean', 'OPZ_232_mean',
       'OPZ_208_mean', 'OPZ_207_mean', 'OPZ_206_mean', 'OPZ_204_mean',
       'OPZ_202_mean', 'OPZ_208/206_mean','OPZ_207/206_mean','OPZ_206/204_Hg-corrected_mean',
        'SNR_238_mean', 'SNR_232_mean',
       'SNR_208_mean', 'SNR_207_mean', 'SNR_206_mean', 'SNR_204_mean',
       'SNR_202_mean', '206/238_before rejection',
       '206/238_before rejection 2se%',  '208/232_before rejection',
       '208/232_before rejection 2se%',  '207/206_before rejection',
       '207/206_before rejection 2se%',  '208/206_before rejection',
       '208/206_before rejection 2se%',  '206/204_before rejection',
       '206/204_before rejection 2se%',  '208/204_before rejection',
       '208/204_before rejection 2se%',  '207/204_before rejection', '207/204_before rejection 2se%',  
       '206/238_after rejection', '206/238_after rejection 2se%','206/238_after rejection 2se',
       '208/232_after rejection', '208/232_after rejection 2se%','208/232_after rejection 2se',
       '207/206_after rejection', '207/206_after rejection 2se%','207/206_after rejection 2se',
       '208/206_after rejection', '208/206_after rejection 2se%','208/206_after rejection 2se',
       '206/204_after rejection', '206/204_after rejection 2se%','206/204_after rejection 2se',
       '208/204_after rejection', '208/204_after rejection 2se%','208/204_after rejection 2se',
       '207/204_after rejection', '207/204_after rejection 2se%','207/204_after rejection 2se',
       '206/238 corrected','206/238_int corrected', '206/238 corrected BB error', '206/238 corrected BB error%', '206/238 corrected Sx%',
       '208/232 corrected', '208/232 corrected BB error','208/232 corrected BB error%', '208/232 corrected Sx%',
       '207/206 corrected', '207/206 corrected BB error', '207/206 corrected BB error%', '207/206 corrected Sx%',
       '208/206 corrected', '208/206 corrected BB error', '208/206 corrected BB error%', '208/206 corrected Sx%',
       '206/204 corrected', '206/204 corrected BB error', '206/204 corrected BB error%', '206/204 corrected Sx%',
       '207/204 corrected', '207/204 corrected BB error', '207/204 corrected BB error%', '207/204 corrected Sx%',
       '208/204 corrected',  '208/204 corrected BB error', '208/204 corrected BB error%', '208/204 corrected Sx%',         
        ]
  
    
    stat2 = stat_new[stat2_order]
    
    return stat2, group_new, result_plot, new_grouper_plot
            

# + code_folding=[19, 68, 211, 214, 254, 257, 293, 296]
### Functions for actually running correction-math, calculating ages, and exporting: calc_all_ages(), correction_to_dict(),

def calc_all_ages(plot):
    ''' Calculates isotope decay system ages for 206/238, 207/206, 207/235, and 208/232'''
     #Calculating Ages
    plot['206/238 Age [Ma]'] = plot.apply(lambda x: calc_t68(x['206/238 corrected'])/ 1E6, axis=1)
    plot['207/206 Age [Ma]'] = plot.apply(lambda x: calc_t76_2(x['207/206 corrected'])/ 1E6, axis=1)
    plot['207/235 Age [Ma]'] = plot.apply(lambda x: calc_t75(x['207/235 corrected'])/ 1E6, axis=1)
    plot['208/232 Age [Ma]'] = plot.apply(lambda x: calc_t82(x['208/232 corrected'])/ 1E6, axis=1)
    plot['238/206'] = plot.apply(lambda x: 1 / x['206/238 corrected'], axis=1)
    plot['238/206_int'] = plot.apply(lambda x: 1 / x['206/238_int corrected'], axis=1)
    
    first_column = plot.pop('238/206')
    plot.insert(0, '238/206', first_column)
    second_column = plot.pop('238/206_int')
    plot.insert(1, '238/206_int', second_column)
    
    return plot

def df_correction_207Pb(ttn_zero_full_plot, iterations = 5):
    '''Performs 5 iterations of 207-Pb correction using Stacey-Kramers on the input plot dataframe.''' #Need to be able to do this with different initial Pb values
    
    df =  ttn_zero_full_plot


    for idx in range (iterations):
        if idx == 0:
            df['206/238_Age'] = df.apply(lambda x: calc_t68(x['206/238 corrected']) / 1E6, axis=1)
            
        else:
            df['206/238_Age'] = df.apply(lambda x: calc_t68(x['206/238_corr']) / 1E6, axis=1)
             
        df['207/206_initial_SK'] = df.apply(lambda x: SK_values(x['206/238_Age'])['SK_207_206'], axis=1)
        df['207/206_rad_expected'] = df.apply(lambda x: calc_concordant_ratios(x['206/238_Age'])['207*/206*'], axis=1)
        df['f_206'] = df.apply(lambda x: (x['207/206 corrected'] - x['207/206_rad_expected']) / (x['207/206_initial_SK'] - x['207/206_rad_expected']), axis=1)
        df['206/238_corr'] = df.apply(lambda x: (1 - x['f_206']) * x['206/238 corrected'], axis=1)
        
    df['208/206_initial_SK'] = df.apply(lambda x: SK_values(x['206/238_Age'])['SK_208_206'], axis=1)   
    df['206/238_Age'] = df.apply(lambda x: calc_t68(x['206/238_corr']) / 1E6, axis=1)
    
    df['207/206_corr'] = df.apply(lambda x: (x['207/206 corrected'] -  x['207/206_initial_SK'] * x['f_206']) / (1 - x['f_206']), axis=1)
    df['208/206_corr'] = df.apply(lambda x: (x['208/206 corrected'] -  x['208/206_initial_SK'] * x['f_206']) / (1 - x['f_206']), axis=1)
    df['238/232_calc'] = df.apply(lambda x: (x['208/232 corrected']) *  (1 / x['208/206 corrected']) * (1 / x['206/238 corrected']), axis=1)
    df['208/232_corr'] = df.apply(lambda x: (x['208/232 corrected']) -  (x['206/238 corrected'] * x['208/206_initial_SK'] * x['238/232_calc'] * x['f_206']),  axis=1)                                       
    df['207/235_corr'] = df.apply(lambda x: x['207/235 corrected']  -  (x['206/238 corrected'] * x['207/206_initial_SK'] * U85r * x['f_206']),  axis=1)

    
    to_keep = ['206/238 corrected','206/238_int corrected', '206/238 corrected Sx%',
       '206/238 corrected BB error%', '206/238_after rejection 2se%',
       '207/206 corrected', '207/206 corrected Sx%',
       '207/206 corrected BB error%', '207/206_after rejection 2se%',
       '208/232 corrected', '208/232 corrected Sx%',
       '208/232 corrected BB error%', '208/232_after rejection 2se%',
       '208/206 corrected', '208/206 corrected Sx%',
       '208/206 corrected BB error%', '208/206_after rejection 2se%',
       '206/204 corrected', '206/204 corrected Sx%',
       '206/204 corrected BB error%', '206/204_after rejection 2se%',
       '207/204 corrected', '207/204 corrected Sx%',
       '207/204 corrected BB error%', '207/204_after rejection 2se%',
       '208/204 corrected', '208/204 corrected Sx%',
       '208/204 corrected BB error%', '208/204_after rejection 2se%',
       '207/235 corrected', '207/235 corrected Sx%',
       '207/235 corrected BB error%', '207/235_after rejection 2se%',
        'f_206',  '206/238_corr', '207/206_corr', '208/206_corr', '208/232_corr', '207/235_corr'     
              ]
    
    return df[to_keep]

def calc_all_ages_corr207(plot):
    ''' Calculates isotope decay system ages for 206/238, 207/206, 207/235, and 208/232 using uncorrected and corrected values.'''
    
    
    
    plot = df_correction_207Pb(plot)
    
    plot = plot.copy()
    
     #Calculating Ages
    plot['206/238 Age [Ma]_uncorr'] = plot.apply(lambda x: calc_t68(x['206/238 corrected'])/ 1E6, axis=1)
    plot['207/206 Age [Ma]_uncorr'] = plot.apply(lambda x: calc_t76_2(x['207/206 corrected'])/ 1E6, axis=1)
    plot['207/235 Age [Ma]_uncorr'] = plot.apply(lambda x: calc_t75(x['207/235 corrected'])/ 1E6, axis=1)
    plot['208/232 Age [Ma]_uncorr'] = plot.apply(lambda x: calc_t82(x['208/232 corrected'])/ 1E6, axis=1)

    plot['206/238 Age [Ma]_corr207'] = plot.apply(lambda x: calc_t68(x['206/238_corr'])/ 1E6, axis=1)
    plot['207/206 Age [Ma]_corr207'] = plot.apply(lambda x: calc_t76_2(x['207/206_corr'])/ 1E6, axis=1)
    plot['207/235 Age [Ma]_corr207'] = plot.apply(lambda x: calc_t75(x['207/235_corr'])/ 1E6, axis=1)
    plot['208/232 Age [Ma]_corr207'] = plot.apply(lambda x: calc_t82(x['208/232_corr'])/ 1E6, axis=1)
    
    
    plot['238/206_corr207'] = plot.apply(lambda x: 1 / x['206/238_corr'], axis=1)
    
    plot['238/206'] = plot.apply(lambda x: 1 / x['206/238 corrected'], axis=1)
    
    first_column = plot.pop('238/206')
    plot.insert(0, '238/206', first_column)
    
    return plot

def correction_to_dict(tester, date, standard):
    '''Performs correction with NIST 612 and MKED-1, at 0%, 20%, error minimization. Outputs a list of dicts: full_export, plot_export, and chronologic_export.'''
    
    print('Primary Standard= ', 'SRM NIST 612')
    print('Error minimization rejection percentage= ', 0)
    excV_dict1 = stat_to_excV(tester, 'SRM NIST 612', 0)
    glass612_zero_full, glass612_zero_groups, glass612_zero_full_plot, glass612_zero_groups_plot  = stat_rank_and_correct(tester, 0, '612',excV_dict1 )
    
    print('Primary Standard= ', 'SRM NIST 612')
    print('Error minimization rejection percentage= ', 20)
    excV_dict2 = stat_to_excV(tester, 'SRM NIST 612', 20)
    glass612_20_full, glass612_20_groups, glass612_20_full_plot, glass612_20_groups_plot = stat_rank_and_correct(tester, 20, '612', excV_dict2)
    #glass612_40_full, glass612_40_groups, glass612_40_full_plot, glass612_40_groups_plot = stat_rank_and_correct(tester, 40, '612')
    
    
    
    print('Primary Standard= ', standard)
    print('Error minimization rejection percentage= ', 0)
    excV_dict3 = stat_to_excV(tester, standard, 0)
    ttn_zero_full, ttn_zero_groups, ttn_zero_full_plot, ttn_zero_groups_plot  = stat_rank_and_correct(tester, 0, standard, excV_dict3)
    
    print('Primary Standard= ', standard)
    print('Error minimization rejection percentage= ', 20)
    excV_dict4 = stat_to_excV(tester, standard, 20)
    ttn_20_full, ttn_20_groups, ttn_20_full_plot, ttn_20_groups_plot = stat_rank_and_correct(tester, 20, standard, excV_dict4)
    #ttn_40_full, ttn_40_groups, ttn_40_full_plot, ttn_40_groups_plot = stat_rank_and_correct(tester, 40, 'ttn')

    
    export_dict = {'glass612_zero_full':glass612_zero_full, 
               'glass612_zero_groups':glass612_zero_groups,
               'glass612_20_full': glass612_20_full, 
               'glass612_20_groups': glass612_20_groups,
               #'glass612_40_full':glass612_40_full, 
               #'glass612_40_groups':glass612_40_groups,
               
               
               'ttn_zero_full': ttn_zero_full, 
               'ttn_zero_groups': ttn_zero_groups,
               'ttn_20_full': ttn_20_full, 
               'ttn_20_groups': ttn_20_groups,
               #'ttn_40_full': ttn_40_full, 
               #'ttn_40_groups':ttn_40_groups
              }

    plot_dict = {
               'glass612_zero_full':calc_all_ages_corr207(glass612_zero_full_plot), 
               'glass612_zero_groups':glass612_zero_groups_plot,
               'glass612_20_full': calc_all_ages_corr207(glass612_20_full_plot), 
               'glass612_20_groups': glass612_20_groups_plot,
               #'glass612_40_full':calc_all_ages(glass612_40_full_plot), 
               #'glass612_40_groups':glass612_40_groups_plot,
             
             
               'ttn_zero_full': calc_all_ages_corr207(ttn_zero_full_plot), 
               'ttn_zero_groups': ttn_zero_groups_plot,
               'ttn_20_full': calc_all_ages_corr207(ttn_20_full_plot), 
               'ttn_20_groups': ttn_20_groups_plot,
               #'ttn_40_full': calc_all_ages(ttn_40_full_plot), 
               #'ttn_40_groups':ttn_40_groups_plot
              }

    chrono_dict = {
            'glass612_zero_full': chronologic(glass612_zero_full, chrono_order), 
            'glass612_zero_groups': glass612_zero_groups,
                }

    return [export_dict, plot_dict, chrono_dict]

def export_corr_full(export_dict, date, bracket = False):
    '''Export for full data worksheets'''

    
    if bracket:
        excel_name = 'Processed_LASS_data_difRejects_difStds_BRACKETS_' + date + '.xlsx'
    else:    
        excel_name = 'Processed_LASS_data_difRejects_difStds_' + date + '.xlsx'
    with pd.ExcelWriter(excel_name, engine = 'xlsxwriter') as writer:
          # Get the xlsxwriter workbook objects.
        workbook  = writer.book
        col_format1 = workbook.add_format()
        col_format2 = workbook.add_format()
        col_format3 = workbook.add_format()
        col_format1.set_bg_color('#E6E6FA') #Lavender
        col_format2.set_bg_color('#98FB98') #Pale Green
        col_format3.set_bg_color('#98FB98') #Antique White


        for sheet in export_dict:
            export_dict[sheet].to_excel(writer, sheet_name = sheet, index = True)
            worksheet = writer.sheets[sheet]  # pull worksheet object
            test_list = list(export_dict[sheet].keys())
            worksheet.set_column(0, 0, 20)
            worksheet.freeze_panes(1,1)
            for idx in enumerate(test_list):
                if 'full' in sheet:
                    if 'BB' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20)
                    elif 'corrected' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format1)
                    elif 'SNR_20' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format2)

                    else:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20)
                elif 'groups' in sheet:
                    if 'BB' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20)
                    elif 'mean' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format2)
                    elif 'MSWD' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format3)
                    elif 'Weighted' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format1)        
                    else:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20)

def export_corr_plot(plot_dict, date, bracket = False):
    '''Export for plotting data worksheets''' 

    if bracket:
        excel_name = 'Processed_LASS_data_difRejects_difStds_BRACKETS_' + date + '_plot.xlsx'
    else: 
        excel_name = 'Processed_LASS_data_difRejects_difStds_' + date + '_plot.xlsx'

    with pd.ExcelWriter(excel_name, engine = 'xlsxwriter') as writer:
        # Get the xlsxwriter workbook objects.
        workbook  = writer.book
        col_format1 = workbook.add_format()
        col_format2 = workbook.add_format()
        col_format3 = workbook.add_format()
        col_format1.set_bg_color('#E6E6FA') #Lavender
        col_format2.set_bg_color('#98FB98') #Pale Green
        col_format3.set_bg_color('#98FB98') #Antique White

        for sheet in plot_dict:

            plot_dict[sheet].to_excel(writer, sheet_name = sheet, index = True)
            worksheet = writer.sheets[sheet]  # pull worksheet object
            test_list = list(plot_dict[sheet].keys())
            worksheet.set_column(0, 0, 20)
            worksheet.freeze_panes(1,1)
            for idx in enumerate(test_list):
                if 'full' in sheet:
                    if 'BB' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20)
                    else:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format1)
                elif 'groups' in sheet:
                    if 'BB' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20)
                    elif 'mean' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format1)
                    elif 'MSWD' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format3)
                    elif 'Weighted' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format2)         
                    else:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20) 
                        
def export_corr_chrono(chrono_dict, date):
    '''Export for chronologic data worksheets'''
    
    excel_name = 'Processed_LASS_data_difRejects_difStds_' + date +'_chrono.xlsx'
    with pd.ExcelWriter(excel_name, engine = 'xlsxwriter') as writer:
          # Get the xlsxwriter workbook objects.
        workbook  = writer.book
        col_format1 = workbook.add_format()
        col_format2 = workbook.add_format()
        col_format1.set_bg_color('#E6E6FA') #Lavender
        col_format2.set_bg_color('#98FB98') #Pale Green

        for sheet in chrono_dict:

            chrono_dict[sheet].to_excel(writer, sheet_name = sheet, index = True)
            worksheet = writer.sheets[sheet]  # pull worksheet object
            test_list = list(chrono_dict[sheet].keys())
            worksheet.set_column(0, 0, 20)
            worksheet.freeze_panes(1,1)
            for idx in enumerate(test_list):
                if 'full' in sheet:
                    if 'BB' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20)
                    elif 'corrected' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format1)
                    elif 'SNR_20' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format2)
                    elif 'OPZ' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format1)   
                    else:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20)
                elif 'groups' in sheet:
                    if 'BB' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20)
                    elif 'mean' in idx[1]:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20, col_format2)
                    else:
                        worksheet.set_column(idx[0]+1, idx[0]+1, 20)   
                                               
def corr_to_excel(tester, date, standard):
    corr_list = correction_to_dict(tester, date, standard)
    print('\nSTEP 6. correction_to_dict() is complete.\n')
    export = corr_list[0]
    plot = corr_list[1]
    #chrono = corr_list[2]
   
    export_corr_full(export, date)
    print('STEP 7. export_corr_full() is complete.\n')
    export_corr_plot(plot, date)
    print('STEP 8. export_corr_plot() is complete.\n')
    #export_corr_chrono(chrono, date)
    
    return plot['ttn_zero_full']



# + code_folding=[2, 26, 61]
### Miscellaneous functions: print_chrono(), export_isoplot(), file_to_output()
chrono_order = []
def print_chrono(chrono_order = chrono_order):
    '''This just prints lists in three's so it is easier to read as output.'''
    for i in enumerate(chrono_order):
        co = chrono_order
        x = 0 + 3 * i[0]
        y = 1 + 3 * i[0]
        z = 2 + 3 * i[0]
        limit = len(chrono_order)
        #print(i)
        
        if (z == limit):
            print(co[x],co[y], sep = ', ')
            print('\n')
            break
        elif (y == limit):
            print(co[x])
            print('\n')
            break
        elif (x == limit):
            print('\n')
            break
        else:
            print(co[x],co[y],co[z], sep = ', ')
                       
def export_isoplot(filename, isoplot, date, bracket = False, Pbcorr = False):
    '''Export for isoplot worksheet; MKED-1 std, 0% error minimization.'''
    
    
    if Pbcorr:
        iso_order = ['238/206','238/206_int','238/206_corr207','206/238_int corrected', '206/238 corrected Sx%',  '207/206 corrected', '207/206_corr', '207/206 corrected Sx%', '206/238 corrected BB error%', 
            '207/206 corrected BB error%', '206/238_after rejection 2se%','207/206_after rejection 2se%', 
             '206/238 Age [Ma]_uncorr', '207/206 Age [Ma]_uncorr', '207/235 Age [Ma]_uncorr','208/232 Age [Ma]_uncorr',
                '206/238 Age [Ma]_corr207', '207/206 Age [Ma]_corr207', '207/235 Age [Ma]_corr207','208/232 Age [Ma]_corr207']
    else:    
        iso_order = ['238/206', '238/206_int','206/238 corrected Sx%',  '207/206 corrected', '207/206 corrected Sx%', '206/238 corrected BB error%', 
            '207/206 corrected BB error%', '206/238_after rejection 2se%','207/206_after rejection 2se%', '206/238_int corrected', 
            '206/238 Age [Ma]', '207/206 Age [Ma]', '207/235 Age [Ma]','208/232 Age [Ma]',
                ]  
            
    isoplot = isoplot[iso_order]
    indices = isoplot.index.values.tolist()

    for name in indices:
        if 'SRM NIST' in name:
            isoplot = isoplot.drop(name, axis=0)
    if bracket:
        excel_name = str(filename.split('.')[0]) + '_isoplot_BRACKETS_' + date + '.xlsx'
    else:    
        excel_name = str(filename.split('.')[0]) + '_isoplot_' + date + '.xlsx'
    with pd.ExcelWriter(excel_name, engine = 'xlsxwriter') as writer:
          # Get the xlsxwriter workbook objects.
        workbook  = writer.book
        isoplot.to_excel(writer, sheet_name = 'Isoplot', index = True)
        worksheet = writer.sheets['Isoplot']  # pull worksheet object
        
        for idx in range(20):
            worksheet.set_column(idx, idx, 22)
        worksheet.freeze_panes(1,1)
                   
def file_to_output(filename, date, order, standard, report_choice = True):
    ''' Just combines all of the steps for data reduction.'''
    print(f'File chosen = {filename} \n\nDate chosen = {date} \n\nChronologic order defined as: ')

    print_chrono(order)

    test_df = read_np2_timeseries(filename)   #ACTUAL CODE

    print('STEP 1. read_np2_timeseries() is complete.\n')
    global tester
    tester = calc_CPS(test_df)     #ACTUAL CODE

    print('STEP 2. calc_CPS() is complete.\n')

    #Visual check correct baseline-subtracted file is being used.
    print('Check that these are correct time-series: \n')

    print_chrono(list(tester.keys()))

    if report_choice: 
        #Fractionation output
        fractionation_output(filename)   #ACTUAL CODE
        print('STEP 3. fractionation_output() is complete.\n')

    #Processing data
    excel_name = str(filename.split('.')[0]) + '_processed.xlsx'    #ACTUAL CODE
    files_process_toEXCEL(tester, excel_name, order)    #ACTUAL CODE
    print('STEP 4. files_process_toEXCEL() is complete.\n')

        
    if report_choice:    
        #Visual report of each line scan
        report_name = 'Splitstream_' + date + '_results'    #ACTUAL CODE
        U_Pb_report(tester, 'SS2_Oct14.pdf', False, report_name)     #ACTUAL CODE; Placeholder .pdf name might cause problems later.
        print('\nSTEP 5. U_Pb_report() is complete.\n')

    #This is a time-consuming function
     
    
    isoplot = corr_to_excel(tester, date, standard)  #ACTUAL CODE            
    
    export_isoplot(filename, isoplot, date)
    print('STEP 9. export_isoplot() is complete.\n')
    

# + code_folding=[]
### Functions for bracketing: bracket(), order_to_brackets(), file_to_blockBrackets(), file_Bracket_to_output()

def bracket(order, input_std = 'MKED-1'):
    '''Function for creating block-brackets based on chronologic order of LASS-ICPMS run and selected bracketing standard. Outputs list of indices for brackets.'''
    idx_list = []
    for idx in range(len(order)):
        if input_std in order[idx]:
            idx_list.append(idx)

    groups_mked = []
    group_list = []
    #print(idx_list)
    for entry in idx_list:
        #print(entry)
        if len(group_list) == 0:
            #print(entry)
            group_list.append(entry)
            continue
        if entry == idx_list[-1]:
            group_list.append(entry)
            groups_mked.append(group_list)
            continue
        if (entry - 1) == group_list[-1]:
            #print(entry)
            group_list.append(entry)
        
        else:
            groups_mked.append(group_list)
            #print(entry)
            group_list = [entry]
    #print(groups_mked)
    brack_list = []

    for idx in range(len(groups_mked) - 1):

        bracket = [groups_mked[idx][0],groups_mked[idx + 1][-1]]
        brack_list.append(bracket)
        
    full_bracket = [groups_mked[0][0],groups_mked[-1][-1]]
    print('full bracket = ', full_bracket)
    return brack_list, full_bracket  

def order_to_brackets(order, input_std = 'MKED-1'):
    '''Takes in the chronologic order and primary standard, and chops the full analysis list into block-brackets. Outputs dict of each bracket.'''
    brackets, full_bracket = bracket(order, input_std)
    
    brack_dict = {}

    for idx in range(len(brackets)):
        #print(idx, len(brackets))
        name = 'bracket_' + str(idx + 1)
        
        if idx == 0:
            
            brack_dict[name] = order[0:(brackets[idx][1] + 1)]
            
        elif idx == (len(brackets) - 1):
            #print('DEBUG')
            brack_dict[name] = order[brackets[idx][0]:]
            #print(brack_dict[name])
        else:
            brack_dict[name] = order[brackets[idx][0]:(brackets[idx][1] + 1)]
        
    #full_brack = order[full_bracket[0]:(full_bracket[1]+1)]
    #print(brack_dict)
    return brack_dict

def file_to_blockBrackets(filename, order, input_std = 'MKED-1'):
    '''Function to take .xlsx time-series file and export a dictionary of the different block-brackets based on input chronologic order and primary standard.'''
    test_df = read_np2_timeseries(filename)   #ACTUAL CODE
    
    stat_test_order = statistics_NP2(test_df).index.values.tolist()   

    new_order = stat_test_order
    
    
    brack_dict = order_to_brackets(order, input_std)
    
    #test_df_brack = {key: test_df[key] for key in full_bracket}
    
    
    bracket_df_dict = {}
    
    for entry in brack_dict:
        
        bracket_df_dict[entry] = {key: test_df[key] for key in brack_dict[entry]}

    return bracket_df_dict, new_order

def file_Bracket_to_output(filename, date, order, primary_std, Pb_corr = False, debug = False):
    ''' Just combines all of the steps for clock-bracket data reduction.'''
    print(f'File chosen = {filename} \n\nDate chosen = {date} \n\nChronologic order defined as: ')
    print_chrono(order)
    brack_df_dict, new_order = file_to_blockBrackets(filename, order, primary_std)
    global tester
    
    #Getting excess variance for all of the standard analyses
    
    test_df = read_np2_timeseries(filename)   #ACTUAL CODE
    tester = calc_CPS(test_df)     #ACTUAL CODE
    #Processing data
    excel_name = str(filename.split('.')[0]) + '_processed.xlsx'    #ACTUAL CODE
    files_process_toEXCEL(tester, excel_name, order)    #ACTUAL CODE
    
    print('Primary Standard= ', primary_std)
    print('Error minimization rejection percentage= ', 0)
    
    excV_dict = stat_to_excV(tester, primary_std, reject_percentage = 0)
    
    print('\nSTEP 1. stat_to_excV() is complete.\n')
    

    ### Make loop that does each bracket
    print('STEP 2. file_to_blockBrackets() is complete.\n')

    concat_export_full = []
    concat_export_group = []
    concat_plot_full = []
    concat_plot_group = []
    
    
    
    for bracket in brack_df_dict:
        test_df = brack_df_dict[bracket]
        
        tester = calc_CPS(test_df)     #ACTUAL CODE

        print('\nSTEP 3. calc_CPS() is complete for ' + bracket + '.\n')

        if primary_std == 'SRM NIST 612':
            std_entry = '612'
        elif primary_std == 'BLR-1':
            std_entry = 'BLR-1'
        else:
            std_entry = 'ttn'
        
   
        ttn_zero_full, ttn_zero_groups, ttn_zero_full_plot, ttn_zero_groups_plot  = stat_rank_and_correct(tester, 0, std_entry, excV_dict)

        concat_export_full.append(ttn_zero_full)
        concat_export_group.append(ttn_zero_groups)
        concat_plot_full.append(ttn_zero_full_plot)
        concat_plot_group.append(ttn_zero_groups_plot)   
    
    conc_export_full = pd.concat(concat_export_full)
    conc_export_group = pd.concat(concat_export_group)
    conc_plot_full = pd.concat(concat_plot_full)
    conc_plot_group = pd.concat(concat_plot_group)
    
    export_dict = {std_entry + '_zero_full': conc_export_full[~conc_export_full.index.duplicated(keep='first')], 
               std_entry + '_zero_groups': conc_export_group.sort_index(ascending=False),
              }

    export_corr_full(export_dict, date, True)
    
    if Pb_corr:
        plot_dict = { std_entry + '_zero_full': calc_all_ages_corr207(conc_plot_full), 
               std_entry + '_zero_groups': conc_plot_group.sort_index(ascending=False),
              }
    else:
        plot_dict = { std_entry + '_zero_full': calc_all_ages(conc_plot_full), 
               std_entry + '_zero_groups': conc_plot_group.sort_index(ascending=False),
              }

    print('STEP 4. Separate brackets have been concatenated.\n') 

   
    print('STEP 5. export_corr_full() is complete.\n')
    export_corr_plot(plot_dict, date, True)
    print('STEP 6. export_corr_plot() is complete.\n')

    isoplot = plot_dict[std_entry + '_zero_full']
    isoplot_new = isoplot[~isoplot.index.duplicated(keep='first')]
    
    #print(new_order)
    new_list = []
    for entry in new_order:
        if entry not in isoplot_new.index.values.tolist():
            continue
        else:
            new_list.append(entry)
    isoplot_new = isoplot_new.loc[new_list]
    
    export_isoplot(filename, isoplot_new, date, True, False)
    print('STEP 7. export_isoplot() is complete.\n')
    
    if debug:
        return isoplot

# + code_folding=[3, 201, 322, 575]
#Regression Nonsense


def regress_indiv2(tester, sample_name, choice = False, start_clip = 0, end_clip = 40, date = ''):
    '''Does a linear regression for each individual analysis and outputs a dictionary with mean, intercept, slope, coefficient of determination, and (intercept/average).'''
    sample = tester[sample_name]
    
    new_string = sample_name + '_' + date
    
    regress_data = sample[['Elapsed Time', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]
    x_orig = np.array(regress_data['Elapsed Time']).reshape((-1, 1))
    mked_full = regress_data[regress_data['Elapsed Time'] > start_clip]
    mked_full = mked_full[mked_full['Elapsed Time'] < end_clip]
    x = np.array(mked_full['Elapsed Time']).reshape((-1, 1))
    
    
    # PLot 1 #
    y = mked_full['206/238']

    model = LinearRegression().fit(x, y)

    model_dict = {}
    model_dict['average 206/238'] = y.mean()
    model_dict['intercept 206/238'] = model.intercept_
    model_dict['slope 206/238'] = float(model.coef_)
    model_dict['coefficient of determination 206/238'] = model.score(x, y)
    model_dict['intercept / average 206/238'] =  model_dict['intercept 206/238'] / model_dict['average 206/238']
    
    y_pred = model.predict(x_orig)

    fig, ax = plt.subplots(1, 2, sharex = True, figsize = (24, 8))
    fig.suptitle(new_string, fontsize=24)
    
    
    ax[0].scatter(x,y, color="black")
    ax[0].plot(x_orig, y_pred, color='blue', linewidth = 3)
    ax[0].axhline(y= y.mean(), color='green', linestyle='-')
    ax[0].set_title('206/238', fontsize = 20)
    
    # these are matplotlib.patch.Patch properties
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

    # place a text box in upper left in axes coords
    ax[0].text(0.97, 0.97, sample_name, transform=ax[0].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation = 'y = ' + str(round_sig(model_dict['slope 206/238'], 3)) + 'x + ' + str(round_sig(model_dict['intercept 206/238']))
    ax[0].text(0.97, 0.88, equation, transform=ax[0].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    # Plot 1 END #
    
    # PLot 2 #
    y2 = mked_full['208/232']

    model2 = LinearRegression().fit(x, y2)

    model_dict2 = {}
    model_dict2['average 208/232'] = y2.mean()
    model_dict2['intercept 208/232'] = model2.intercept_
    model_dict2['slope 208/232'] = float(model2.coef_)
    model_dict2['coefficient of determination 208/232'] = model2.score(x, y2)
    model_dict2['intercept / average 208/232'] =  model_dict2['intercept 208/232'] / model_dict2['average 208/232']
    
    y_pred2 = model2.predict(x_orig)

    
    ax[1].scatter(x,y2, color="black")
    ax[1].plot(x_orig, y_pred2, color='blue', linewidth = 3)
    ax[1].axhline(y = y2.mean(), color='green', linestyle='-')
    ax[1].set_title('208/232', fontsize = 20)
    
    # these are matplotlib.patch.Patch properties
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

    # place a text box in upper left in axes coords
    ax[1].text(0.97, 0.97, sample_name, transform=ax[1].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation = 'y = ' + str(round_sig(model_dict2['slope 208/232'], 3)) + 'x + ' + str(round_sig(model_dict2['intercept 208/232']))
    ax[1].text(0.97, 0.88, equation, transform=ax[1].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    # Plot 2 END #
    
    
    plt.xlim([0, end_clip])
    
    
    
    
   
    
    MYDIR = ("Regression Figures" + date)
    CHECK_FOLDER = os.path.isdir(MYDIR)

    # If folder doesn't exist, then create it.
    if not CHECK_FOLDER:
        os.makedirs(MYDIR)
        #print("created folder : ", MYDIR)
    
    #new_string = sample.replace('time series data', '').rstrip()
    
    filename = os.path.join(MYDIR, new_string + '.pdf')
    plt.savefig(filename)
    print('Plot for ', new_string, " is complete.")
    
    if choice == False:
        plt.close()
    
    
    
    
    #plt.close()
    
    #if choice:
        #print('\nRegression for', sample_name)
        #print('average:', model_dict['average'])
        #print('intercept:', new)
        #print('slope:', float(model.coef_))
       #print('coefficient of determination:', model.score(x, y))
       #print('intercept / average', new/orig)
    
    return model_dict, model_dict2


### Code for outputting full regression report

# tester_list = [tester,tester2, tester3,tester4]
# dates = ['Nov16', 'Nov17','Nov22','Nov23']



# full_dict = {}
# for i in range(len(tester_list)):
# #     if i == 1:
# #         break
        
#     int_dict_206 = {}
#     int_dict_208 = {}
#     z_dict = {}
#     testy = tester_list[i]
    
#     output_name = 'Regress_' + dates[i]
#     #MYDIR = ("Regression Figures"+dates[i])
#     #mergedObject = PdfFileMerger()
    
    
    
    
#     for sample in testy:

        
#         first, second  = regress_indiv2(testy, sample, choice = False, start_clip = 2, date = dates[i])
#         int_dict_206[sample] = first
#         int_dict_208[sample] = second
#         #fit = model.fit(x,y)
    
#         z_dict[sample] = {**first, **second}
    

#         new_string = sample + '_' + dates[i]
#         #filename = os.path.join(MYDIR, new_string + '.pdf')

#         #mergedObject.append(PdfFileReader(filename, 'rb'))

#     if '.pdf' in output_name:
#         pass
#     else:
#         output_name = output_name + '.pdf'
    
#     #output_name = "U-Pb_output.pdf"  
#     #mergedObject.write(output_name)

#     #print(f'PDF file named: {output_name} is complete.')    

#     full_dict[dates[i]] = z_dict
        
        
# nov16 = pd.DataFrame(full_dict['Nov16']).transpose()
# nov17 = pd.DataFrame(full_dict['Nov17']).transpose()
# nov22 = pd.DataFrame(full_dict['Nov22']).transpose()
# nov23 = pd.DataFrame(full_dict['Nov23']).transpose()
# nov_dict = {'nov16':nov16,'nov17': nov17,'nov22': nov22,'nov23': nov23}


# excel_name = 'regress1.xlsx'
# with pd.ExcelWriter(excel_name, engine = 'xlsxwriter') as writer:
#       # Get the xlsxwriter workbook objects.
#     workbook  = writer.book
    
#     for sheet in nov_dict:
#         nov_dict[sheet].to_excel(writer, sheet_name = sheet, index = True)

#         worksheet = writer.sheets[sheet]  # pull worksheet object

#         for idx in range(20):
#             worksheet.set_column(idx, idx, 22)
#         worksheet.freeze_panes(1,1)
        
        
def regress_compare(tester, std1, std2, clip1 = 0, clip2  = 0, endclip = 35, choice = False):    
    '''This function takes in the output from calc_CPS(), and two selected sample groups, plus clip positions for start and end. It will plot things.'''
    ###MKED###
    #start_clip = 2
    mked = group_samples2(tester)[std1]

    tester_list = []
    for entry in mked:
        tester_list.append(tester[entry])

    mked_full = pd.concat(tester_list)
    regress_data = mked_full[['Elapsed Time', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]

    x_orig = np.array(mked_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    mked_full = regress_data[regress_data['Elapsed Time'] > clip1]

    mked_full = mked_full[mked_full['Elapsed Time'] < endclip]
    
    x = np.array(mked_full['Elapsed Time']).reshape((-1, 1))
    y = mked_full['206/238']
    model = LinearRegression().fit(x, y)
    model_dict = {}
    model_dict['average'] = y.mean()
    model_dict['intercept'] = model.intercept_
    model_dict['slope'] = float(model.coef_)
    model_dict['coefficient of determination'] = model.score(x, y)

    fit = model.fit(x,y)
    y_pred = model.predict(x_orig)

    ### MKED ###

    fig, ax = plt.subplots(1, 3, sharex = True, figsize = (24, 6))
    ax[0].scatter(x,y, color="black")
    ax[0].plot(x_orig, y_pred, color='blue', linewidth = 3)
    
    # these are matplotlib.patch.Patch properties
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

    # place a text box in upper left in axes coords
    ax[0].text(0.97, 0.97, std1, transform=ax[0].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation = 'y = ' + str(round_sig(model_dict['slope'], 3)) + 'x + ' + str(round_sig(model_dict['intercept']))
    ax[0].text(0.97, 0.88, equation, transform=ax[0].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    

    ### BLR-1 ###
    #start_clip = 8
    blr = group_samples2(tester)[std2]

    tester_list = []
    for entry in blr:
        tester_list.append(tester[entry])

    blr_full = pd.concat(tester_list)
    regress_data = blr_full[['Elapsed Time', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]


    x_orig2 = np.array(blr_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    blr_full = regress_data[regress_data['Elapsed Time'] > clip2]
    blr_full = blr_full[blr_full['Elapsed Time'] < endclip]
    x2 = np.array(blr_full['Elapsed Time']).reshape((-1, 1))
    y2 = blr_full['206/238']
    model = LinearRegression().fit(x2, y2)

    model_dict2 = {}
    model_dict2['average'] = y2.mean()
    model_dict2['intercept'] = model.intercept_
    model_dict2['slope'] = float(model.coef_)
    model_dict2['coefficient of determination'] = model.score(x, y)


    fit2 = model.fit(x2,y2)
    y_pred2 = model.predict(x_orig2)

    ### BLR-1 ###
    ax[1].scatter(x2,y2, color="black")
    ax[1].plot(x_orig2, y_pred2, color='blue', linewidth = 3)

    # place a text box in upper left in axes coords
    ax[1].text(0.97, 0.97, std2, transform=ax[1].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation2 = 'y = ' + str(round_sig(model_dict2['slope'], 3)) + 'x + ' + str(round_sig(model_dict2['intercept']))
    ax[1].text(0.97, 0.88, equation2, transform=ax[1].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    
    y3 = y_pred2 / (y_pred[:len(y_pred2)])
    y4 = model_dict2['intercept'] / (y_pred[:len(y_pred2)])

    ax[2].plot(y3[0:50], color='green', linewidth = 3)
    ax[2].plot(y4[0:50], color='blue', linewidth = 3)
    
    
#     #Still in Development
#     ax[0].spines['left'].set_position(('data', 0))
#     ax[1].spines['left'].set_position(('data', 0))
#     ax[2].spines['left'].set_position(('data', 0))
    
    ###
    #print(equation)
    #print(model_dict)
    #print(model_dict2)
    print('int2/int1: ',model_dict2['intercept']/model_dict['intercept'])
    print('int2/avg1 ', model_dict2['intercept']/model_dict['average'])
    print(std1, ' chopped at: ',clip1, '   ', model_dict2['intercept']/y_pred[clip1])
    print()
    if choice:
        return model_dict, model_dict2
    
    
################################
def regress_compare2(tester, std1, std2,std3, std4, clip1 = 0, clip2  = 0, clip3 = 0, clip4 = 0, endclip = 35, choice = False):    
    '''This function takes in the output from calc_CPS(), and plots 4 standards with std1 = primary, plus clip positions for start and end. It will plot things.'''
    
   
    #d1 = tester['MKED-1 1.10']['238_CPS']

    #dt = 200
    #dt_206_fac = (1 / (1 - t1*(dt/1E9)))
    
    
    ###MKED###
    #start_clip = 2
    mked = group_samples2(tester)[std1]

    tester_list = []
    for entry in mked:
        tester_list.append(tester[entry])

    mked_full = pd.concat(tester_list)
    regress_data = mked_full[['Elapsed Time','206_CPS', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]

    x_orig = np.array(mked_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    mked_full = regress_data[regress_data['Elapsed Time'] > clip1]

    mked_full = mked_full[mked_full['Elapsed Time'] < endclip]
    
    
    
    
    x = np.array(mked_full['Elapsed Time']).reshape((-1, 1))
    y = mked_full['206/238']
    t1 = mked_full['206_CPS']
    #y = y1 * (1 / (1 - t1*(dt/1E9)))
    model = LinearRegression().fit(x, y)
    model_dict = {}
    model_dict['average'] = y.mean()
    model_dict['intercept'] = model.intercept_
    model_dict['slope'] = float(model.coef_)
    model_dict['coefficient of determination'] = model.score(x, y)

    fit = model.fit(x,y)
    y_pred = model.predict(x_orig)

    ### MKED ###

    fig, ax = plt.subplots(1, 4, sharex = True, figsize = (28, 6))
    ax[0].scatter(x,y, color="black")
    ax[0].plot(x_orig, y_pred, color='blue', linewidth = 3)
    ax[0].axhline(y= y.mean(), color='green', linestyle='-')
    
    
    # these are matplotlib.patch.Patch properties
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

    # place a text box in upper left in axes coords
    ax[0].text(0.97, 0.97, std1, transform=ax[0].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation = 'y = ' + str(round_sig(model_dict['slope'], 3)) + 'x + ' + str(round_sig(model_dict['intercept']))
    ax[0].text(0.97, 0.88, equation, transform=ax[0].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    

    ### BLR-1 ###
    #start_clip = 8
    blr = group_samples2(tester)[std2]

    tester_list = []
    for entry in blr:
        tester_list.append(tester[entry])

    blr_full = pd.concat(tester_list)
    regress_data = blr_full[['Elapsed Time','206_CPS', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]


    x_orig2 = np.array(blr_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    blr_full = regress_data[regress_data['Elapsed Time'] > clip2]
    blr_full = blr_full[blr_full['Elapsed Time'] < endclip]
    x2 = np.array(blr_full['Elapsed Time']).reshape((-1, 1))
    y2 = blr_full['206/238']
    model = LinearRegression().fit(x2, y2)

    model_dict2 = {}
    model_dict2['average'] = y2.mean()
    model_dict2['intercept'] = model.intercept_
    model_dict2['slope'] = float(model.coef_)
    model_dict2['coefficient of determination'] = model.score(x, y)


    fit2 = model.fit(x2,y2)
    y_pred2 = model.predict(x_orig2)

    ### BLR-1 ###
    ax[1].scatter(x2,y2, color="black")
    ax[1].plot(x_orig2, y_pred2, color='blue', linewidth = 3)
    ax[1].axhline(y= y2.mean(), color='green', linestyle='-')
    
    
    
    
    # place a text box in upper left in axes coords
    ax[1].text(0.97, 0.97, std2, transform=ax[1].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation2 = 'y = ' + str(round_sig(model_dict2['slope'], 3)) + 'x + ' + str(round_sig(model_dict2['intercept']))
    ax[1].text(0.97, 0.88, equation2, transform=ax[1].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    ### Bancroft ###
    #start_clip = 8
    banc = group_samples2(tester)[std3]

    tester_list = []
    for entry in banc:
        tester_list.append(tester[entry])

    banc_full = pd.concat(tester_list)
    regress_data = banc_full[['Elapsed Time','206_CPS', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]


    x_orig3 = np.array(banc_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    banc_full = regress_data[regress_data['Elapsed Time'] > clip3]
    banc_full = banc_full[banc_full['Elapsed Time'] < endclip]
    x3 = np.array(banc_full['Elapsed Time']).reshape((-1, 1))
    y3 = banc_full['206/238']
    model = LinearRegression().fit(x3, y3)

    model_dict3 = {}
    model_dict3['average'] = y3.mean()
    model_dict3['intercept'] = model.intercept_
    model_dict3['slope'] = float(model.coef_)
    model_dict3['coefficient of determination'] = model.score(x3, y3)


    fit3 = model.fit(x3,y3)
    y_pred3 = model.predict(x_orig3)

    ### Bancroft ###
    ax[2].scatter(x3,y3, color="black")
    ax[2].plot(x_orig3, y_pred3, color='blue', linewidth = 3)
    ax[2].axhline(y= y3.mean(), color='green', linestyle='-')
    
    # place a text box in upper left in axes coords
    ax[2].text(0.97, 0.97, std3, transform=ax[2].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation2 = 'y = ' + str(round_sig(model_dict3['slope'], 3)) + 'x + ' + str(round_sig(model_dict3['intercept']))
    ax[2].text(0.97, 0.88, equation2, transform=ax[2].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    ### Bear Lake ###
    #start_clip = 8
    bear = group_samples2(tester)[std4]

    tester_list = []
    for entry in bear:
        tester_list.append(tester[entry])

    bear_full = pd.concat(tester_list)
    regress_data = bear_full[['Elapsed Time', '206_CPS','206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]

    x_orig4 = np.array(bear_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    bear_full = regress_data[regress_data['Elapsed Time'] > clip4]
    bear_full = bear_full[bear_full['Elapsed Time'] < endclip]
    x4 = np.array(bear_full['Elapsed Time']).reshape((-1, 1))
    y4 = bear_full['206/238']
    model = LinearRegression().fit(x4, y4)

    model_dict4 = {}
    model_dict4['average'] = y4.mean()
    model_dict4['intercept'] = model.intercept_
    model_dict4['slope'] = float(model.coef_)
    model_dict4['coefficient of determination'] = model.score(x4, y4)


    fit4 = model.fit(x4,y4)
    y_pred4 = model.predict(x_orig4)

    ### Bear Lake ###
    ax[3].scatter(x4,y4, color="black")
    ax[3].plot(x_orig4, y_pred4, color='blue', linewidth = 3)
    
    
    ax[3].axhline(y= y4.mean(), color='green', linestyle='-')
    
    
    
    #ax[3].plot(y4.mean(), color='green', linewidth = 3)
    
    # place a text box in upper left in axes coords
    ax[3].text(0.97, 0.97, std4, transform=ax[3].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation2 = 'y = ' + str(round_sig(model_dict4['slope'], 3)) + 'x + ' + str(round_sig(model_dict4['intercept']))
    ax[3].text(0.97, 0.88, equation2, transform=ax[3].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    
    
#     y3 = y_pred2 / (y_pred[:len(y_pred2)])
#     y4 = model_dict2['intercept'] / (y_pred[:len(y_pred2)])

#     ax[2].plot(y3[0:50], color='green', linewidth = 3)
#     ax[2].plot(y4[0:50], color='blue', linewidth = 3)
    
    
#     #Still in Development
#     ax[0].spines['left'].set_position(('data', 0))
#     ax[1].spines['left'].set_position(('data', 0))
#     ax[2].spines['left'].set_position(('data', 0))
    
    ###
    #print(equation)
    #print(model_dict)
    #print(model_dict2)
    space = '                          '
   
    blr_val = model_dict2['intercept']/model_dict['intercept'] * published('MKED-1', '206/238')
    banc_val = model_dict3['intercept']/model_dict['intercept'] * published('MKED-1', '206/238')
    bear_val = model_dict4['intercept']/model_dict['intercept'] * published('MKED-1', '206/238')
    
    blr_age = round_sig(calc_t68(blr_val) /1E6, 6)
    banc_age = round_sig(calc_t68(banc_val) /1E6, 6)
    bear_age = round_sig(calc_t68(bear_val) /1E6, 6)
    
    plt.xlim([0, endclip])
    
    
    print('int/MKED-1 int: ',space,round_sig(blr_val,6),space,round_sig(banc_val,6),space,round_sig(bear_val,6) )
    print('Ages [Ma]:           ',space,round_sig(blr_age,6),space,round_sig(banc_age,6),space,round_sig(bear_age,6) )
    
    #print('int2/MKED-1 int: ',model_dict3['intercept']/model_dict['intercept'], )
    #print('int3/MKED-1 int: ',model_dict4['intercept']/model_dict['intercept'], )
    print() 
    #print('int2/avg1 ', model_dict2['intercept']/model_dict['average'])
    #print(std1, ' chopped at: ',clip1, '   ', model_dict2['intercept']/y_pred[clip1])
    #print()
    if choice:
        return model_dict, model_dict2

    
def regress_compare3(tester, std1, std2,std3, std4, clip1 = 0, clip2  = 0, clip3 = 0, clip4 = 0, endclip = 35, choice = False):    
    '''This function takes in the output from calc_CPS(), and plots 4 standards with std1 = primary, plus clip positions for start and end. It incorporates dead time adjustments on 206 and/or 238.'''
     #Dead Time tests

    dt = 0
    #dt_206_fac = (1 / (1 - t1*(dt/1E9)))
    ut = -70
    
    ###MKED###
    #start_clip = 2
    mked = group_samples2(tester)[std1]

    tester_list = []
    for entry in mked:
        tester_list.append(tester[entry])

    mked_full = pd.concat(tester_list)
    regress_data = mked_full[['Elapsed Time','206_CPS','238_CPS', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]

    x_orig = np.array(mked_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    mked_full = regress_data[regress_data['Elapsed Time'] > clip1]

    mked_full = mked_full[mked_full['Elapsed Time'] < endclip]
    
    
    
    
    x = np.array(mked_full['Elapsed Time']).reshape((-1, 1))
    y1 = mked_full['206/238']
    t1 = mked_full['206_CPS']#DT
    u1 = mked_full['238_CPS']#DT
    y = y1 * (1 / (1 - t1*(dt/1E9))) * (1 /(1 / (1 - u1*(ut/1E9))) )   #DT
    model = LinearRegression().fit(x, y)
    model_dict = {}
    model_dict['average'] = y.mean()
    model_dict['intercept'] = model.intercept_
    model_dict['slope'] = float(model.coef_)
    model_dict['coefficient of determination'] = model.score(x, y)

    fit = model.fit(x,y)
    y_pred = model.predict(x_orig)

    ### MKED ###

    fig, ax = plt.subplots(1, 4, sharex = True, figsize = (28, 6))
    ax[0].scatter(x,y, color="black")
    ax[0].plot(x_orig, y_pred, color='blue', linewidth = 3)
    ax[0].axhline(y= y.mean(), color='green', linestyle='-')
    
    
    # these are matplotlib.patch.Patch properties
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

    # place a text box in upper left in axes coords
    ax[0].text(0.97, 0.97, std1, transform=ax[0].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation = 'y = ' + str(round_sig(model_dict['slope'], 3)) + 'x + ' + str(round_sig(model_dict['intercept']))
    ax[0].text(0.97, 0.88, equation, transform=ax[0].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    

    ### BLR-1 ###
    #start_clip = 8
    blr = group_samples2(tester)[std2]

    tester_list = []
    for entry in blr:
        tester_list.append(tester[entry])

    blr_full = pd.concat(tester_list)
    regress_data = blr_full[['Elapsed Time','206_CPS','238_CPS', '206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]


    x_orig2 = np.array(blr_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    blr_full = regress_data[regress_data['Elapsed Time'] > clip2]
    blr_full = blr_full[blr_full['Elapsed Time'] < endclip]
    x2 = np.array(blr_full['Elapsed Time']).reshape((-1, 1))
    
  
    
    y2 = blr_full['206/238']
    t2 = blr_full['206_CPS']      #DT
    #y2 = y2 * (1 / (1 - t2*(dt/1E9)))    #DT
    u2 = blr_full['238_CPS']#DT
    y2 = y2 * (1 / (1 - t2*(dt/1E9))) * (1 /(1 / (1 - u2*(ut/1E9))) )   #DT
    
    
    model = LinearRegression().fit(x2, y2)

    model_dict2 = {}
    model_dict2['average'] = y2.mean()
    model_dict2['intercept'] = model.intercept_
    model_dict2['slope'] = float(model.coef_)
    model_dict2['coefficient of determination'] = model.score(x, y)


    fit2 = model.fit(x2,y2)
    y_pred2 = model.predict(x_orig2)

    ### BLR-1 ###
    ax[1].scatter(x2,y2, color="black")
    ax[1].plot(x_orig2, y_pred2, color='blue', linewidth = 3)
    ax[1].axhline(y= y2.mean(), color='green', linestyle='-')
    
    
    
    
    # place a text box in upper left in axes coords
    ax[1].text(0.97, 0.97, std2, transform=ax[1].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation2 = 'y = ' + str(round_sig(model_dict2['slope'], 3)) + 'x + ' + str(round_sig(model_dict2['intercept']))
    ax[1].text(0.97, 0.88, equation2, transform=ax[1].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    ### Bancroft ###
    #start_clip = 8
    banc = group_samples2(tester)[std3]

    tester_list = []
    for entry in banc:
        tester_list.append(tester[entry])

    banc_full = pd.concat(tester_list)
    regress_data = banc_full[['Elapsed Time','206_CPS', '238_CPS','206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]


    x_orig3 = np.array(banc_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    banc_full = regress_data[regress_data['Elapsed Time'] > clip3]
    banc_full = banc_full[banc_full['Elapsed Time'] < endclip]
    x3 = np.array(banc_full['Elapsed Time']).reshape((-1, 1))
    
    
    y3 = banc_full['206/238']
    
    t3 = banc_full['206_CPS']      #DT
    
    u3 = banc_full['238_CPS']#DT
    y3 = y3 * (1 / (1 - t3*(dt/1E9))) * (1 /(1 / (1 - u3*(ut/1E9))) )   #DT
    #y3 = y3 * (1 / (1 - t3*(dt/1E9)))    #DT
    
    model = LinearRegression().fit(x3, y3)

    model_dict3 = {}
    model_dict3['average'] = y3.mean()
    model_dict3['intercept'] = model.intercept_
    model_dict3['slope'] = float(model.coef_)
    model_dict3['coefficient of determination'] = model.score(x3, y3)


    fit3 = model.fit(x3,y3)
    y_pred3 = model.predict(x_orig3)

    ### Bancroft ###
    ax[2].scatter(x3,y3, color="black")
    ax[2].plot(x_orig3, y_pred3, color='blue', linewidth = 3)
    ax[2].axhline(y= y3.mean(), color='green', linestyle='-')
    
    # place a text box in upper left in axes coords
    ax[2].text(0.97, 0.97, std3, transform=ax[2].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation2 = 'y = ' + str(round_sig(model_dict3['slope'], 3)) + 'x + ' + str(round_sig(model_dict3['intercept']))
    ax[2].text(0.97, 0.88, equation2, transform=ax[2].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    ### Bear Lake ###
    #start_clip = 8
    bear = group_samples2(tester)[std4]

    tester_list = []
    for entry in bear:
        tester_list.append(tester[entry])

    bear_full = pd.concat(tester_list)
    regress_data = bear_full[['Elapsed Time', '206_CPS','238_CPS','206/238', '208/232', '207/206','208/206', '206/204', '208/204', '207/204', '207/235']]

    x_orig4 = np.array(bear_full['Elapsed Time']).reshape((-1, 1))

    ### Regression Math ###
    #start_clip = 1.0
    bear_full = regress_data[regress_data['Elapsed Time'] > clip4]
    bear_full = bear_full[bear_full['Elapsed Time'] < endclip]
    x4 = np.array(bear_full['Elapsed Time']).reshape((-1, 1))
    y4 = bear_full['206/238']
    
    t4 = bear_full['206_CPS']      #DT
    u4 = bear_full['238_CPS']#DT
    y4 = y4 * (1 / (1 - t4*(dt/1E9))) * (1 /(1 / (1 - u4*(ut/1E9))) )   #DT
    #y4 = y4 * (1 / (1 - t4*(dt/1E9)))    #DT
    
    model = LinearRegression().fit(x4, y4)

    model_dict4 = {}
    model_dict4['average'] = y4.mean()
    model_dict4['intercept'] = model.intercept_
    model_dict4['slope'] = float(model.coef_)
    model_dict4['coefficient of determination'] = model.score(x4, y4)


    fit4 = model.fit(x4,y4)
    y_pred4 = model.predict(x_orig4)

    ### Bear Lake ###
    ax[3].scatter(x4,y4, color="black")
    ax[3].plot(x_orig4, y_pred4, color='blue', linewidth = 3)
    
    
    ax[3].axhline(y= y4.mean(), color='green', linestyle='-')
    
    
    
    #ax[3].plot(y4.mean(), color='green', linewidth = 3)
    
    # place a text box in upper left in axes coords
    ax[3].text(0.97, 0.97, std4, transform=ax[3].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    equation2 = 'y = ' + str(round_sig(model_dict4['slope'], 3)) + 'x + ' + str(round_sig(model_dict4['intercept']))
    ax[3].text(0.97, 0.88, equation2, transform=ax[3].transAxes, fontsize=14,
    verticalalignment='top',horizontalalignment = 'right', bbox=props)
    
    

    
    space = '                          '
    blr_val = model_dict2['intercept']/model_dict['intercept'] * published('MKED-1', '206/238')
    banc_val = model_dict3['intercept']/model_dict['intercept'] * published('MKED-1', '206/238')
    bear_val = model_dict4['intercept']/model_dict['intercept'] * published('MKED-1', '206/238')
    
    blr_age = round_sig(calc_t68(blr_val) /1E6, 6)
    banc_age = round_sig(calc_t68(banc_val) /1E6, 6)
    bear_age = round_sig(calc_t68(bear_val) /1E6, 6)
    
    
    plt.xlim([0, 50])
    
    print('int/MKED-1 int: ',space,round_sig(blr_val,6),space,round_sig(banc_val,6),space,round_sig(bear_val,6) )
    
    #print('int2/MKED-1 int: ',model_dict3['intercept']/model_dict['intercept'], )
    #print('int3/MKED-1 int: ',model_dict4['intercept']/model_dict['intercept'], )
    print() 
    #print('int2/avg1 ', model_dict2['intercept']/model_dict['average'])
    #print(std1, ' chopped at: ',clip1, '   ', model_dict2['intercept']/y_pred[clip1])
    #print()
    if choice:
        return model_dict, model_dict2

# ### This is the normal regression comparison    
# space = '                          '
# print('         ','MKED-1 ',space,'BLR-1',space,'Bancroft',space,'Bear Lake')

# print('       ','Expected:',space,'0.17921',space,'0.178523',space,'0.1870\n')

# #Nov16, Nov17, Nov22, Nov23



# regress_compare2(tester, 'MKED-1', 'BLR-1','Bancroft', 'Bear Lake', clip1 = 2, clip2  =9,clip3 = 7, clip4 = 14, endclip = 34, choice = False)
# regress_compare2(tester2, 'MKED-1', 'BLR-1','Bancroft', 'Bear Lake', clip1 = 2, clip2  = 9,clip3 = 7, clip4 = 14, endclip = 34, choice = False)
# regress_compare2(tester3, 'MKED-1', 'BLR-1','Bancroft', 'Bear Lake', clip1 = 2, clip2  = 9,clip3 = 7, clip4 = 2, endclip = 34, choice = False)
# regress_compare2(tester4, 'MKED-1', 'BLR-1','Bancroft', 'Bear Lake', clip1 = 2, clip2  =9,clip3 = 7, clip4 = 5, endclip = 34, choice = False)


# ### This is the (unjustified) deadtime adjustment regression comparisons 

# space = '                          '
# print('         ','MKED-1 ',space,'BLR-1',space,'Bancroft',space,'Bear Lake')

# print('       ','Expected:',space,'0.17921',space,'0.178523',space,'0.1870\n')

# #Nov16, Nov17, Nov22, Nov23



# regress_compare3(tester, 'MKED-1', 'BLR-1','Bancroft', 'Bear Lake', clip1 = 5, clip2  =8,clip3 = 8, clip4 = 2, endclip = 30, choice = False)
# regress_compare3(tester2, 'MKED-1', 'BLR-1','Bancroft', 'Bear Lake', clip1 = 10, clip2  = 8,clip3 = 5, clip4 = 2, endclip = 37, choice = False)
# regress_compare3(tester3, 'MKED-1', 'BLR-1','Bancroft', 'Bear Lake', clip1 = 2, clip2  = 10,clip3 = 5, clip4 = 2, endclip = 30, choice = False)
# regress_compare3(tester4, 'MKED-1', 'BLR-1','Bancroft', 'Bear Lake', clip1 = 5, clip2  = 5,clip3 = 5, clip4 = 2, endclip = 30, choice = False)
            
# -

# ****Start Here****
#
# Input Filename of Iolite-baseline corrected Time-Series .xlsx, the date, and the order in which things were ablated.
#
# Potential work-flow to consider is doing a full-run average, identifying standards that are outliers, removing those, correcting the order, and re-running with block-brackets.

# + code_folding=[2]
#Check order in Iolite

chrono_order2 = [    #This order is currently manually input based on analysis order - KD 11/3/2021
 'SRM NIST 612 1.1',
 'SRM NIST 612 1.2',
 'SRM NIST 612 1.3',
 'SRM NIST 612 1.4',
 'SRM NIST 612 1.5',
 'MKED-1 1.1',
 'MKED-1 1.2',
 'MKED-1 1.3',
 'MKED-1 1.4',
 'MKED-1 1.5',  
 'BLR-1 1.1',
 'BLR-1 1.2', 
 'BLR-1 1.3',
 'BLR-1 1.4', 
 'BLR-1 1.5',
 '95-65_17 1.1',  
 '95-65_17 1.2',
 '95-65_17 1.3',    
 '95-65_18 1.1',
 '95-65_18 1.2',  
 '95-65_18 1.3', 
 '95-65_21 1.1',  
 '95-65_21 1.2',  
 'Bancroft 1.1',
 'Bancroft 1.2',
 'Bancroft 1.3',
 'Yates 1.1',
 'Yates 1.2', 
 'Yates 1.3',
 'Bear Lake 1.1',
 'Bear Lake 1.2', 
 'Bear Lake 1.3',  
 'SRM NIST 612 1.6',
 'SRM NIST 612 1.7',
 'SRM NIST 612 1.8',
 'MKED-1 1.6',
 'MKED-1 1.7',
 #'MKED-1 1.8',
 'BLR-1 1.6',
 'BLR-1 1.7', 
 'BLR-1 1.8',
 '95-65_29 1.1',  
 '95-65_29 1.2',   
 '95-65_29 1.3',  
 '95-65_29 1.4', 
 '95-65_27 1.1',   
 '95-65_27 1.2',
 '95-65_27 1.3',     
 '95-65_28 1.1',
 '95-65_28 1.2',  
 '95-65_28 1.3',  
 '95-65_28 1.4', 
 'Bancroft 1.4',
 'Bancroft 1.5',
 'Bancroft 1.6',
 'Yates 1.4',
 'Yates 1.5', 
 'Yates 1.6',
 'Bear Lake 1.4',
 'Bear Lake 1.5', 
 'Bear Lake 1.6',    
 'SRM NIST 612 1.9',
 'SRM NIST 612 1.10',
 'SRM NIST 612 1.11',
 'MKED-1 1.9',
 #'MKED-1 1.10',
 'MKED-1 1.11',
 'BLR-1 1.9',
 'BLR-1 1.10', 
 'BLR-1 1.11',
 '95-65_25 1.1',  
 '95-65_25 1.2',   
 '95-65_25 1.3',  
 '95-65_26 1.1',
 '95-65_26 1.2',  
 '95-65_26 1.3', 
 'SRM NIST 612 1.12',
 'SRM NIST 612 1.13',
 'SRM NIST 612 1.14',
 'MKED-1 1.12',
 'MKED-1 1.13',
 'MKED-1 1.14',
 'BLR-1 1.12',
 'BLR-1 1.13', 
 'BLR-1 1.14',  
 'Bancroft 1.7',
 'Bancroft 1.8',
 'Bancroft 1.9',
 'Yates 1.7',
 'Yates 1.8', 
 'Yates 1.9',
 'Bear Lake 1.7',
 'Bear Lake 1.8', 
 'Bear Lake 1.9',    
 '95-69_42 1.1',   
 '95-69_42 1.2',
 '95-69_41 1.1',   
 '95-69_41 1.2',   
 '95-69_41 1.3',
 '95-69_43 1.1',  
 '95-69_43 1.2',  
 'BLR-1 1.15',
 'BLR-1 1.16', 
 'BLR-1 1.17',   
 'BLR-1 1.18',
 'BLR-1 1.19',  
 'MKED-1 1.15',
 'MKED-1 1.16',
 'MKED-1 1.17',
 'MKED-1 1.18',
 'MKED-1 1.19', 
 'SRM NIST 612 1.15',
 'SRM NIST 612 1.16',
 'SRM NIST 612 1.17',
 'SRM NIST 612 1.18',
 'SRM NIST 612 1.19',
 ]


#Load in Baseline-subtracted file and set Date
filename = 'NP2_20211117_baseline-subtracted.xlsx'
date = '17Nov2021'
#order = chrono_order2
primary_std = 'MKED-1'
openIsoplotR = True

# + code_folding=[3, 107, 225, 386, 554, 744]
### DO NOT DELETE ###
#Order of analyses#

oct22_order = [
 'SRM NIST 612 1.1',
 'SRM NIST 612 1.2',
 'SRM NIST 612 1.3',
 'SRM NIST 612 1.4',
 'SRM NIST 612 1.5',
 'MKED-1 1.1',
 'MKED-1 1.2',
 'MKED-1 1.3',
 'MKED-1 1.4',
 'MKED-1 1.5',  
 'BLR-1 1.1',
 'BLR-1 1.2', 
 'BLR-1 1.3',
 'BLR-1 1.4', 
 'BLR-1 1.5', 
 '95-69_t30 1.1',
 '95-69_t30 1.2',
 '95-69_t30 1.3',
 '95-69_t1 1.1',
 '95-69_t1 1.2',
 '95-69_t1 1.3', 
 '95-69_t2 1.1',
 '95-69_t2 1.2',
 '95-69_t2 1.3',   
 '95-69_t6 1.1',
 '95-69_t6 1.2',
 'SRM NIST 612 1.6',
 'SRM NIST 612 1.7',
 'SRM NIST 612 1.8',
 'MKED-1 1.6',
 'MKED-1 1.7',
 'MKED-1 1.8',
 'BLR-1 1.6',
 'BLR-1 1.7', 
 'BLR-1 1.8',
 '95-69_t7 1.1',
 '95-69_t7 1.2',
 '95-69_t7 1.3',  
 '95-69_t8 1.1',
 '95-69_t8 1.2',
 '95-69_t8 1.3',    
 '95-69_t15 1.1',
 '95-69_t15 1.2',
 '95-69_t15 1.3', 
 '95-69_t3 1.1',
 '95-69_t3 1.2',
 '95-69_t3 1.3',    
 'Bancroft 1.1',
 'Bancroft 1.2',
 'Bancroft 1.3',   
 'SRM NIST 612 1.9',
 'SRM NIST 612 1.10',
 'SRM NIST 612 1.11',
 'MKED-1 1.9',
 'MKED-1 1.10',
 #'MKED-1 1.11',
 'BLR-1 1.9',
 'BLR-1 1.10', 
 'BLR-1 1.11',
 '95-69_t4 1.1',
 '95-69_t4 1.2',
 '95-69_t4 1.3', 
 '95-69_t5 1.1',
 '95-69_t5 1.2',
 '95-69_t5 1.3', 
 '95-65_t1 1.1',
 '95-65_t1 1.2',
 '95-65_t1 1.3',  
 '95-65_t2 1.1',
 '95-65_t2 1.2',
 '95-65_t2 1.3',   
 'Bancroft 1.4',
 'Bancroft 1.5',
 'Bancroft 1.6',  
 'SRM NIST 612 1.12',
 'SRM NIST 612 1.13',
 'SRM NIST 612 1.14',
 'MKED-1 1.12',
 'MKED-1 1.13',
 'MKED-1 1.14',
 'BLR-1 1.12',
 'BLR-1 1.13', 
 'BLR-1 1.14',
 '95-65_t3 1.1',
 '95-65_t3 1.2',
 '95-65_t3 1.3',   
 '95-65_t4 1.1',
 '95-65_t4 1.2',
 '95-65_t4 1.3',   
 '95-65_t15 1.1',
 '95-65_t15 1.2',
 '95-65_t15 1.3',   
 '95-65_t16 1.1',
 '95-65_t16 1.2',
 '95-65_t17 1.1',  
 'SRM NIST 612 1.15',  # Might be bad
 'SRM NIST 612 1.16',  # Might be bad
 #'SRM NIST 612 1.17',  # Might be bad
 #'MKED-1 1.15', # Might be bad
 #'MKED-1 1.16', # Might be bad
 #'MKED-1 1.17', # Might be bad
    
 ]
nov4_order = [
 'SRM NIST 612 1.1',
 'SRM NIST 612 1.2',
 'SRM NIST 612 1.3',
 'SRM NIST 612 1.4',
 'SRM NIST 612 1.5',
 'MKED-1 1.1',
 'MKED-1 1.2',
 'MKED-1 1.3',
 'MKED-1 1.4',
 'MKED-1 1.5',  
 'BLR-1 1.1',
 'BLR-1 1.2', 
 'BLR-1 1.3',
 'BLR-1 1.4', 
 'BLR-1 1.5',
 '95-65_17 1.1',  
 '95-65_17 1.2',
 '95-65_17 1.3',    
 '95-65_18 1.1',
 '95-65_18 1.2',  
 '95-65_18 1.3', 
 '95-65_21 1.1',  
 '95-65_21 1.2',  
 'Bancroft 1.1',
 'Bancroft 1.2',
 'Bancroft 1.3',
 'Yates 1.1',
 'Yates 1.2', 
 'Yates 1.3',
 'Bear Lake 1.1',
 'Bear Lake 1.2', 
 'Bear Lake 1.3',  
 'SRM NIST 612 1.6',
 'SRM NIST 612 1.7',
 'SRM NIST 612 1.8',
 'MKED-1 1.6',
 'MKED-1 1.7',
 #'MKED-1 1.8',
 'BLR-1 1.6',
 'BLR-1 1.7', 
 'BLR-1 1.8',
 '95-65_29 1.1',  
 '95-65_29 1.2',   
 '95-65_29 1.3',  
 '95-65_29 1.4', 
 '95-65_27 1.1',   
 '95-65_27 1.2',
 '95-65_27 1.3',     
 '95-65_28 1.1',
 '95-65_28 1.2',  
 '95-65_28 1.3',  
 '95-65_28 1.4', 
 'Bancroft 1.4',
 'Bancroft 1.5',
 'Bancroft 1.6',
 'Yates 1.4',
 'Yates 1.5', 
 'Yates 1.6',
 'Bear Lake 1.4',
 'Bear Lake 1.5', 
 'Bear Lake 1.6',    
 'SRM NIST 612 1.9',
 'SRM NIST 612 1.10',
 'SRM NIST 612 1.11',
 'MKED-1 1.9',
 #'MKED-1 1.10',
 'MKED-1 1.11',
 'BLR-1 1.9',
 'BLR-1 1.10', 
 'BLR-1 1.11',
 '95-65_25 1.1',  
 '95-65_25 1.2',   
 '95-65_25 1.3',  
 '95-65_26 1.1',
 '95-65_26 1.2',  
 '95-65_26 1.3', 
 'SRM NIST 612 1.12',
 'SRM NIST 612 1.13',
 'SRM NIST 612 1.14',
 'MKED-1 1.12',
 'MKED-1 1.13',
 'MKED-1 1.14',
 'BLR-1 1.12',
 'BLR-1 1.13', 
 'BLR-1 1.14',  
 'Bancroft 1.7',
 'Bancroft 1.8',
 'Bancroft 1.9',
 'Yates 1.7',
 'Yates 1.8', 
 'Yates 1.9',
 'Bear Lake 1.7',
 'Bear Lake 1.8', 
 'Bear Lake 1.9',    
 '95-69_42 1.1',   
 '95-69_42 1.2',
 '95-69_41 1.1',   
 '95-69_41 1.2',   
 '95-69_41 1.3',
 '95-69_43 1.1',  
 '95-69_43 1.2',  
 'BLR-1 1.15',
 'BLR-1 1.16', 
 'BLR-1 1.17',   
 'BLR-1 1.18',
 'BLR-1 1.19',  
 'MKED-1 1.15',
 'MKED-1 1.16',
 'MKED-1 1.17',
 'MKED-1 1.18',
 'MKED-1 1.19', 
 'SRM NIST 612 1.15',
 'SRM NIST 612 1.16',
 'SRM NIST 612 1.17',
 'SRM NIST 612 1.18',
 'SRM NIST 612 1.19',
 ]
nov16_order = [
     'SRM NIST 612 1.1',
 'SRM NIST 612 1.2',
 'SRM NIST 612 1.3',
 'SRM NIST 612 1.4',
 'SRM NIST 612 1.5',
 'MKED-1 1.1',
 'MKED-1 1.2',
 'MKED-1 1.3',
 'MKED-1 1.4',
 'MKED-1 1.5',  
 'BLR-1 1.1',
 'BLR-1 1.2', 
 'BLR-1 1.3',
 'BLR-1 1.4', 
 'BLR-1 1.5',
 '94-77_1 1.1',  
 '94-77_1 1.2',
 '94-77_1 1.3',
 '94-77_1 1.4', 
 '94-77_1 1.5',  
 '94-77_t2 1.1',  
 '94-77_t2 1.2',
 '94-77_t2 1.3',  
 '94-77_t2 1.4',  
 '94-77_t3 1.1',  
 '94-77_t3 1.2',
 '94-77_t3 1.3',  
 '94-77_t3 1.4',  
    
 'SRM NIST 612 1.6',
 'SRM NIST 612 1.7',
 'SRM NIST 612 1.8',
 'MKED-1 1.6',
 'MKED-1 1.7',
 'MKED-1 1.8',
 'BLR-1 1.6',
 'BLR-1 1.7', 
 'BLR-1 1.8', 
 '94-77_t4 1.1',  
 '94-77_t4 1.2',
 '94-77_t4 1.3',  
 '94-77_t5 1.1',  
 '94-77_t5 1.2',
 '94-77_t5 1.3',  
 '94-77_t5 1.4',  
 '94-77_t5 1.5',  
 '94-77_t6 1.1',  
 '94-77_t6 1.2',
 '94-77_t6 1.3',   
 '94-77_t7 1.1',  
 '94-77_t7 1.2',
 '94-77_t7 1.3', 
    
  
 'Yates 1.1',
 'Yates 1.2', 
 'Yates 1.3',
 'Bear Lake 1.1',
 'Bear Lake 1.2', 
 'Bear Lake 1.3',    
 'Bancroft 1.1',
 'Bancroft 1.2',
 'Bancroft 1.3',  
 'SRM NIST 612 1.9',
 'SRM NIST 612 1.10',
 'SRM NIST 612 1.11',
 'SRM NIST 612 1.12',  
 'MKED-1 1.9',
 'MKED-1 1.10',
 'MKED-1 1.11',
 'MKED-1 1.12',   
 'BLR-1 1.9',
 'BLR-1 1.10', 
 'BLR-1 1.11',  
 'BLR-1 1.12',  
    
 '94-77_t8 1.1',  
 '94-77_t8 1.2',
 '94-77_t8 1.3',   
 '94-77_t9 1.1',  
 '94-77_t9 1.2',
 '94-77_t9 1.3',  
 '94-77_t10 1.1',  
 '94-77_t10 1.2',
 '94-77_t10 1.3',   
 '94-77_t10 1.4',  
 '94-77_t11 1.1',  
 '94-77_t11 1.2',
 '94-77_t11 1.3',   
 '94-77_t11 1.4', 
 '94-77_t12 1.1',  
 '94-77_t12 1.2',
 '94-77_t12 1.3',  
 '94-77_t13 1.1',  
 '94-77_t13 1.2', 
    
 'Yates 1.4',
 'Yates 1.5', 
 'Yates 1.6',
 'Bear Lake 1.4',
 'Bear Lake 1.5', 
 'Bear Lake 1.6',    
 'Bancroft 1.4',
 'Bancroft 1.5',
 'Bancroft 1.6',  
 'SRM NIST 612 1.13',
 'SRM NIST 612 1.14',
 'SRM NIST 612 1.15',
 'SRM NIST 612 1.16',  
 'MKED-1 1.13',
 'MKED-1 1.14',
 'MKED-1 1.15',
 'MKED-1 1.16',   
 'BLR-1 1.13',
 'BLR-1 1.14', 
 'BLR-1 1.15',  
 'BLR-1 1.16',     
    
 '91-39_t7 1.1',  
 '91-39_t7 1.2',
 '91-39_t7 1.3',  
 '91-39_t8 1.1',  
 '91-39_t8 1.2',
 '91-39_t8 1.3',  
 '91-39_t4 1.1',  
 '91-39_t4 1.2',  
 '91-39_t3 1.1',  
 '91-39_t3 1.2',
 '91-39_t3 1.3',  
 '91-39_t3 1.4',  
 '91-39_t2 1.1',  
 '91-39_t2 1.2',
 '91-39_t2 1.3',    
 '91-39_t1 1.1',   
 '91-39_t6 1.1',  
  
 'Yates 1.7',
 'Yates 1.8', 
 'Yates 1.9',
 'Bear Lake 1.7',
 'Bear Lake 1.8', 
 'Bear Lake 1.9',    
 'Bancroft 1.7',
 'Bancroft 1.8',
 'Bancroft 1.9',
 'BLR-1 1.17',
 'BLR-1 1.18', 
 'BLR-1 1.19',  
 'BLR-1 1.20',     
 'BLR-1 1.21',  
 'MKED-1 1.17',
 'MKED-1 1.18', 
 'MKED-1 1.19',  
 'MKED-1 1.20',     
 'MKED-1 1.21',   
 'SRM NIST 612 1.17',
 'SRM NIST 612 1.18',
 'SRM NIST 612 1.19',
 'SRM NIST 612 1.20',   
 'SRM NIST 612 1.21',   ]
nov17_order = [
 'SRM NIST 612 1.1',
 'SRM NIST 612 1.2',
 'SRM NIST 612 1.3',
 'SRM NIST 612 1.4',
 'SRM NIST 612 1.5',
 'MKED-1 1.1',
 'MKED-1 1.2',
 'MKED-1 1.3',
 'MKED-1 1.4',
 'MKED-1 1.5',  
 'BLR-1 1.1',
 'BLR-1 1.2', 
 'BLR-1 1.3',
 'BLR-1 1.4', 
 'BLR-1 1.5',
    
 'W15-1_t5 1.1',
 'W15-1_t5 1.2',
 'W15-1_t5 1.3',
 'W15-1_t5 1.4', 
 'W15-1_t4 1.1',
 'W15-1_t4 1.2',
 'W15-1_t4 1.3',
 'W15-1_t4 1.4',      
 'W15-1_t2 1.1',
 'W15-1_t2 1.2',
 'W15-1_t2 1.3',
 'W15-1_t3a 1.1',
 'W15-1_t3a 1.2',
 'W15-1_t3a 1.3',
 'W15-1_t3b 1.1',
 'W15-1_t3b 1.2',
 'W15-1_t3b 1.3',    
 
 'SRM NIST 612 1.6',
 'SRM NIST 612 1.7',
 'SRM NIST 612 1.8',
 'MKED-1 1.6',
 'MKED-1 1.7',
 'MKED-1 1.8',
 'BLR-1 1.6',
 'BLR-1 1.7', 
 'BLR-1 1.8',
  
 'W15-1_t1 1.1',
 'W15-1_t1 1.2',
 'W15-1_t1 1.3',
 'W15-1_t1 1.4',   
 'W15-1_t11 1.1',
 'W15-1_t11 1.2',
 'W15-1_t11 1.3',
 'W15-1_t12a 1.1',
 'W15-1_t12a 1.2',
 'W15-1_t12a 1.3',  
 'W15-1_t12b 1.1',
 'W15-1_t8 1.1',
 'W15-1_t8 1.2',
 'W15-1_t8 1.3',
 'W15-1_t7 1.1',
 'W15-1_t7 1.2',
 'W15-1_t6 1.1',
 'W15-1_t6 1.2', 
    
 'Yates 1.1',
 'Yates 1.2', 
 'Yates 1.3',
 'Bear Lake 1.1',
 'Bear Lake 1.2', 
 'Bear Lake 1.3',    
 'Bancroft 1.1',
 'Bancroft 1.2',
 'Bancroft 1.3',  
 'SRM NIST 612 1.9',
 'SRM NIST 612 1.10',
 'SRM NIST 612 1.11',
 'SRM NIST 612 1.12',  
 'MKED-1 1.9',
 'MKED-1 1.10',
 'MKED-1 1.11',
 'MKED-1 1.12',   
 'BLR-1 1.9',
 'BLR-1 1.10', 
 'BLR-1 1.11',  
 'BLR-1 1.12',    
    
 #'W15-1_t9 1.1',
 'W15-1_t9 1.2',
 #'W15-1_t9 1.3',   
 'W15-1_t10 1.1',
 'W15-1_t10 1.2',   
 'W15-3_t11 1.1',
 'W15-3_t11 1.2',
 'W15-3_t11 1.3',
 'W15-3_t10 1.1',
 'W15-3_t10 1.2',
 'W15-3_t10 1.3',    
 'W15-3_t8 1.1',
 'W15-3_t8 1.2',  
 'W15-3_t9 1.1',
 'W15-3_t9 1.2',
 'W15-3_t9 1.3',
 'W15-3_t5 1.1',
 'W15-3_t5 1.2',   
    
 'Yates 1.4',
 'Yates 1.5', 
 'Yates 1.6',
 'Bear Lake 1.4',
 'Bear Lake 1.5', 
 'Bear Lake 1.6',    
 'Bancroft 1.4',
 'Bancroft 1.5',
 'Bancroft 1.6',  
 'SRM NIST 612 1.13',
 'SRM NIST 612 1.14',
 'SRM NIST 612 1.15',
 'SRM NIST 612 1.16',  
 'MKED-1 1.13',
 'MKED-1 1.14',
 'MKED-1 1.15',
 'MKED-1 1.16',   
 'BLR-1 1.13',
 'BLR-1 1.14', 
 'BLR-1 1.15',  
 'BLR-1 1.16',   
    
 'W15-3_t4 1.1',
 'W15-3_t4 1.2',
 'W15-3_t4 1.3', 
 'W15-3_t4 1.4',   
 'W15-3_t1 1.1',
 'W15-3_t1 1.2',
 'W15-3_t1 1.3', 
 'W15-3_t1 1.4',
 'W15-3_t2 1.1',
 'W15-3_t2 1.2',
 'W15-3_t2 1.3', 
 'W15-3_t2 1.4',
 'W15-3_t3 1.1',
 'W15-3_t3 1.2',
 'W15-3_t3 1.3', 
 'W15-3_t3 1.4', 
    
 'Yates 1.7',
 'Yates 1.8', 
 'Yates 1.9',
 'Bear Lake 1.7',
 'Bear Lake 1.8', 
 'Bear Lake 1.9',    
 'Bancroft 1.7',
 'Bancroft 1.8',
 'Bancroft 1.9',
 'BLR-1 1.17',
 'BLR-1 1.18', 
 'BLR-1 1.19',  
 #'BLR-1 1.20',     
 #'BLR-1 1.21',  
 'MKED-1 1.17',
 'MKED-1 1.18', 
 'MKED-1 1.19',  
 'MKED-1 1.20',     
 'MKED-1 1.21',   
 'SRM NIST 612 1.17',
 'SRM NIST 612 1.18',
 'SRM NIST 612 1.19',
 'SRM NIST 612 1.20',   
 'SRM NIST 612 1.21',   ]
nov22_order = [
    'SRM NIST 612 1.1',
 'SRM NIST 612 1.2',
 'SRM NIST 612 1.3',
 'SRM NIST 612 1.4',
 'SRM NIST 612 1.5',
 'MKED-1 1.1',
 'MKED-1 1.2',
 'MKED-1 1.3',
 'MKED-1 1.4',
 'MKED-1 1.5',  
 'BLR-1 1.1',
 'BLR-1 1.2', 
 'BLR-1 1.3',
 'BLR-1 1.4', 
 'BLR-1 1.5',
 'SAB-94-134_t8 1.1',   
 'SAB-94-134_t8 1.2',
 'SAB-94-134_t8 1.3',
 'SAB-94-134_t10 1.1',   
 'SAB-94-134_t10 1.2',
 'SAB-94-134_t10 1.3',   
 'SAB-94-134_t10 1.4',  
 'SAB-94-134_t10 1.5',  
 'SAB-94-134_t7 1.1',  
 'SAB-94-134_t7 1.2',  
 'SAB-94-134_t7 1.3',  
 'SAB-94-134_t6 1.1',  
 'SAB-94-134_t6 1.2',  
 'SAB-94-134_t6 1.3',   
 'SAB-94-134_t5 1.1',  
 'SAB-94-134_t5 1.2',  
 'SAB-94-134_t5 1.3',   
 
 'SRM NIST 612 1.6',
 'SRM NIST 612 1.7',
 'SRM NIST 612 1.8',
 'MKED-1 1.6',
 'MKED-1 1.7',
 'MKED-1 1.8',
 'BLR-1 1.6',
 'BLR-1 1.7', 
 'BLR-1 1.8',  
  
 'SAB-94-134_t3 1.1',  
 'SAB-94-134_t3 1.2',  
 'SAB-94-134_t3 1.3',   
 'SAB-94-134_t3 1.4',   
 'SAB-94-134_t4 1.1',  #ln46
 'SAB-94-134_t4 1.2',  #ln47
 'SAB-94-134_t4 1.3',  
 'SAB-94-134_t4 1.4', 
 'SAB-94-134_t2 1.1',  #ln50
 'SAB-94-134_t2 1.2',  #ln51
 'SAB-94-134_t2 1.3',   
 'SAB-94-134_t2 1.4',  
 'SAB-94-134_t1 1.1',  
 'SAB-94-134_t1 1.2',  
 'SAB-94-134_t1 1.3',   
 'SAB-94-134_t1 1.4',  
 'SAB-94-134_t1 1.5',  
  
 'Bear Lake 1.1',
 'Bear Lake 1.2', 
 'Bear Lake 1.3',    
 'Bancroft 1.1',
 'Bancroft 1.2',
 'Bancroft 1.3',
 'Yates 1.1',
 'Yates 1.2', 
 'Yates 1.3',   
 'SRM NIST 612 1.9',
 'SRM NIST 612 1.10',
 'SRM NIST 612 1.11',
 'SRM NIST 612 1.12',  
 'MKED-1 1.9',
 'MKED-1 1.10',
 'MKED-1 1.11',
 'MKED-1 1.12',   
 'BLR-1 1.9',
 'BLR-1 1.10', 
 'BLR-1 1.11',  
 'BLR-1 1.12', 
  
 'SAB-95-65c_t5 1.1',   
 'SAB-95-65c_t5 1.2',    
 'SAB-95-65c_t6 1.1',    
 'SAB-95-65c_t7 1.1',   
 'SAB-95-65c_t7 1.2',   
 'SAB-95-65c_t8 1.1',   
 'SAB-95-65c_t9 1.1',   
 'SAB-95-65c_t9 1.2',   
 'SAB-95-65c_t9 1.3',  
 'SAB-95-65c_t10 1.1',   
 'SAB-95-65c_t10 1.2',   
 'SAB-95-65c_t10 1.3',   
  
 'Bear Lake 1.4',
 'Bear Lake 1.5', 
 'Bear Lake 1.6',    
 'Bancroft 1.4',
 'Bancroft 1.5',
 'Bancroft 1.6',   
 'Yates 1.4',
 'Yates 1.5', 
 'Yates 1.6',  
 'SRM NIST 612 1.13',
 'SRM NIST 612 1.14',
 'SRM NIST 612 1.15',
 'SRM NIST 612 1.16',  
 'MKED-1 1.13',
 'MKED-1 1.14',
 'MKED-1 1.15',
 'MKED-1 1.16',   
 'BLR-1 1.13',
 'BLR-1 1.14', 
 'BLR-1 1.15',  
 'BLR-1 1.16',   
 
 'SAB-95-65c_t11 1.1',   
 'SAB-95-65c_t11 1.2',   
 'SAB-95-65c_t14 1.1',   
 'SAB-95-65c_t14 1.2', 
 'SAB-95-65c_t12 1.1',   
 'SAB-95-65c_t12 1.2',  
 'SAB-95-65c_t13 1.1',   
 'SAB-95-65c_t13 1.2',   
 'SAB-95-65c_t13 1.3',   
 'SAB-95-65b_t1 1.1',   
 'SAB-95-65b_t2 1.1',   
 'SAB-95-65b_t2 1.2', 
 #'SAB-95-65b_t3 1.1',   
 'SAB-95-65b_t4 1.1',   
 'SAB-95-65b_t4 1.2',   
 'SAB-95-65b_t5 1.1',   
 'SAB-95-65b_t5 1.2',  
 'SAB-95-65b_t5 1.3',  
    
 'Bear Lake 1.7',
 'Bear Lake 1.8', 
 'Bear Lake 1.9',    
 'Bancroft 1.7',
 'Bancroft 1.8',
 'Bancroft 1.9',
 'Yates 1.7',
 'Yates 1.8', 
 'Yates 1.9',   
 'SRM NIST 612 1.17',
 'SRM NIST 612 1.18',
 'SRM NIST 612 1.19',
 'SRM NIST 612 1.20',  
 'MKED-1 1.17',
 'MKED-1 1.18',
 'MKED-1 1.19',
 'MKED-1 1.20',   
 'BLR-1 1.17',
 'BLR-1 1.18', 
 'BLR-1 1.19',  
 'BLR-1 1.20',   
  
 'SAB-95-65b_t11 1.1',   
 'SAB-95-65b_t11 1.2',  
 'SAB-95-65b_t9 1.1',
 'SAB-95-65b_t9 1.2', 
 'SAB-95-65b_t9 1.3',  
 'SAB-95-65b_t10 1.1',
 'SAB-95-65b_t10 1.2',  
 'SAB-95-65b_t10 1.3',  
 'SAB-95-65b_t12 1.1',
 'SAB-95-65b_t12 1.2',  
 'SAB-95-65b_t12 1.3',  
 'SAB-95-65b_t13 1.1',
 'SAB-95-65b_t13 1.2',  
 'SAB-95-65b_t13 1.3',
  
 'BLR-1 1.21',
 'BLR-1 1.22', 
 'BLR-1 1.23',  
 'BLR-1 1.24',
 'BLR-1 1.25',  
 'MKED-1 1.21',
 'MKED-1 1.22',
 'MKED-1 1.23',
 'MKED-1 1.24',   
 'MKED-1 1.25',   
 'SRM NIST 612 1.21',
 'SRM NIST 612 1.22',
 'SRM NIST 612 1.23',
 'SRM NIST 612 1.24',  
 'SRM NIST 612 1.25',    ]
nov23_order = [
    'SRM NIST 612 1.1',
 'SRM NIST 612 1.2',
 'SRM NIST 612 1.3',
 'SRM NIST 612 1.4',
 'SRM NIST 612 1.5',
 'MKED-1 1.1',
 'MKED-1 1.2',
 'MKED-1 1.3',
 'MKED-1 1.4',
 'MKED-1 1.5',  
 'BLR-1 1.1',
 'BLR-1 1.2', 
 'BLR-1 1.3',
 'BLR-1 1.4', 
 'BLR-1 1.5',
    
 'SB-91-53_t6 1.1',
 'SB-91-53_t6 1.2',
 'SB-91-53_t6 1.3',   
 'SB-91-53_t4 1.1',
 'SB-91-53_t4 1.2',
 'SB-91-53_t4 1.3',   
 'SB-91-53_t4 1.4',   
 'SB-91-53_t5 1.1',
 'SB-91-53_t5 1.2',  
 'SB-91-53_t3 1.1',
 'SB-91-53_t3 1.2',
    
 'Bear Lake 1.1',
 'Bear Lake 1.2', 
 'Bear Lake 1.3',    
 'Bancroft 1.1',
 'Bancroft 1.2',
 'Bancroft 1.3',
 'Yates 1.1',
 'Yates 1.2', 
 'Yates 1.3',   
 'SRM NIST 612 1.6',
 'SRM NIST 612 1.7',
 'SRM NIST 612 1.8',
 'SRM NIST 612 1.8',  
 'MKED-1 1.6',
 'MKED-1 1.7',
 'MKED-1 1.8',
 'MKED-1 1.9',   
 'BLR-1 1.6',
 'BLR-1 1.7', 
 'BLR-1 1.8',  
 'BLR-1 1.9',
    
 'SB-91-53_t2 1.1',
 'SB-91-53_t2 1.2',   
 'SB-91-53_t1 1.1',
 'SB-91-53_t1 1.2',
 'SB-91-53_t1 1.3',  
 'SB-91-53_t12 1.1',
 'SB-91-53_t12 1.2',
 'SB-91-53_t13 1.1',   
 'SB-91-53_t7 1.1',
 'SB-91-53_t7 1.2',
 'SB-91-53_t7 1.3',   
 'SB-91-53_t14 1.1',   
 'SB-91-53_t10 1.1',
 'SB-91-53_t10 1.2',
 'SB-91-53_t10 1.3',   
    
 'Bear Lake 1.4',
 'Bear Lake 1.5', 
 'Bear Lake 1.6',    
 'Bancroft 1.4',
 'Bancroft 1.5',
 'Bancroft 1.6',   
 'Yates 1.4',
 'Yates 1.5', 
 'Yates 1.6', 
 'SRM NIST 612 1.10',
 'SRM NIST 612 1.11',
 'SRM NIST 612 1.12',
 'SRM NIST 612 1.13',  
 'MKED-1 1.10',
 'MKED-1 1.11',
 'MKED-1 1.12',
 'MKED-1 1.13',   
 'BLR-1 1.10',
 'BLR-1 1.11', 
 'BLR-1 1.12',  
 'BLR-1 1.13',
      
 'SB-91-53_t8 1.1',
 'SB-91-53_t8 1.2',   
 'SB-91-53_t8 1.3',
 'SB-91-53_t8 1.4',  
 'SB-91-53_t8 1.5',
 'SB-91-53_t11 1.1',
 'SB-91-53_t11 1.2',   
 'SB-91-53_t11 1.3',  
 'SB-91-53_t15 1.1',
 'SB-91-53_t15 1.2',   
 'SB-91-53_t15 1.3', 
 'SB-91-53_t16 1.1',
 'SB-91-53_t16 1.2',   
 'SB-91-53_t16 1.3',  
    
 'Bear Lake 1.7',
 'Bear Lake 1.8', 
 'Bear Lake 1.9',    
 'Bancroft 1.7',
 'Bancroft 1.8',
 'Bancroft 1.9',
 'Yates 1.7',
 'Yates 1.8', 
 'Yates 1.9',  
 'BLR-1 1.14',
 'BLR-1 1.15', 
 'BLR-1 1.16',  
 'BLR-1 1.17',
 'BLR-1 1.18',  
 'MKED-1 1.14',
 #'MKED-1 1.15',
 'MKED-1 1.16',
 #'MKED-1 1.17',   
# 'MKED-1 1.18',   
 'SRM NIST 612 1.14',
 'SRM NIST 612 1.15',
 'SRM NIST 612 1.16',
 'SRM NIST 612 1.17',  
 'SRM NIST 612 1.18',     ]

dec09_order = [
 #'SRM NIST 612 1.1',
 'SRM NIST 612 1.2',
 'SRM NIST 612 1.3',
 'SRM NIST 612 1.4',
 'SRM NIST 612 1.5',
 'MKED-1 1.1',
 'MKED-1 1.2',
 'MKED-1 1.3',
 'MKED-1 1.4',
 'MKED-1 1.5',  
 'BLR-1 1.1',
 'BLR-1 1.2', 
 'BLR-1 1.3',
 'BLR-1 1.4', 
 'BLR-1 1.5', 
    
 'Bear Lake 1.1',   
 'Bear Lake 1.2',
 'Bear Lake 1.3',   
 '95-69_t1 1.1',  
 '95-69_t1 1.2',  
 '95-69_t1 1.3',    
 '95-69_t8 1.1',  
 '95-69_t8 1.2',  
 '95-69_t8 1.3',     
 #'Yates 1.1', 
 #'Yates 1.2', 
 #'Yates 1.3',    
 'Bancroft 1.1',
 'Bancroft 1.2',
 'Bancroft 1.3',  
    
 'SRM NIST 612 1.6',
 'SRM NIST 612 1.7',
 'SRM NIST 612 1.8',
 'SRM NIST 612 1.9',   
 'MKED-1 1.6',
 'MKED-1 1.7',
 'MKED-1 1.8',
 'MKED-1 1.9',
 'BLR-1 1.6',
 'BLR-1 1.7', 
 'BLR-1 1.8',
 'BLR-1 1.9',
    
 'Bear Lake 1.4',   
 'Bear Lake 1.5',
 'Bear Lake 1.6',     
 '95-65_t3 1.1',  
 '95-65_t3 1.2',  
 '95-65_t3 1.3',   
 '95-65_t4 1.1',  
 '95-65_t4 1.2',  
 '95-65_t4 1.3',     
 #'Yates 1.4', 
 #'Yates 1.5', 
 #'Yates 1.6',    
 'Bancroft 1.4',
 'Bancroft 1.5',
 #'Bancroft 1.6',    
   
 'BLR-1 1.10',
 'BLR-1 1.11', 
 'BLR-1 1.12',
 'BLR-1 1.13', 
 'BLR-1 1.14',   
 'MKED-1 1.10',
 'MKED-1 1.11',
 'MKED-1 1.12',
 'MKED-1 1.13',
 'MKED-1 1.14',     
 'SRM NIST 612 1.10',
 'SRM NIST 612 1.11',
 'SRM NIST 612 1.12',
 'SRM NIST 612 1.13',
 'SRM NIST 612 1.14',

 ]

# +
# Need to fix how the rows are ordered...
#date = '29Nov2021'
filename = 'NP2_20211022_baseline-subtracted_mod1.xlsx'
filename2 = 'NP2_20211104_baseline-subtracted_mod1.xlsx'
filename3 = 'NP2_20211116_baseline-subtracted.xlsx'
filename4 = 'NP2_20211122_baseline-subtracted.xlsx'
filename5 = 'NP2_20211123_baseline-subtracted_mod1.xlsx'

filename6 = 'NP2_20211208_baseline-subtracted.xlsx'


# -

file_Bracket_to_output(filename6, '09Dec', dec09_order, 'MKED-1')

file_Bracket_to_output(filename, '22Oct', oct22_order, 'BLR-1')

file_Bracket_to_output(filename2, '4Nov', nov4_order, 'MKED-1')

file_Bracket_to_output(filename2, '4Nov', nov4_order, 'BLR-1')

file_Bracket_to_output(filename3, '16Nov', nov16_order, 'MKED-1')

file_Bracket_to_output(filename4, '22Nov', nov22_order, 'MKED-1')

file_Bracket_to_output(filename5, '23Nov', nov23_order, 'MKED-1')

#Currently will open IsoplotR in Safari Browser.
if openIsoplotR:
    webbrowser.get("safari").open_new("http://isoplotr.geol.ucsb.edu/isoplotr/")

# +
### Code to run everything in a FULL GROUP right now ###

file_to_output(filename, date, order, primary_std, True)

#Currently will open IsoplotR in Safari Browser.
if openIsoplotR:
    webbrowser.get("safari").open_new("http://isoplotr.geol.ucsb.edu/isoplotr/")

# + code_folding=[0]
# #### DEPRECATED FUNCTIONS as of 3 November 2021 #####

# #These might be useful later, but are currently not being used.

#ratio_report(plot_dict, 'fake.pdf', False, 'ratio_output.pdf')

# def files_ranked_toEXCEL(calc_dict, excel_name):
#     stats = statistics_ranktest(calc_dict)
#     with pd.ExcelWriter(excel_name) as writer:
#         for sheet in calc_dict:
#             calc_dict[sheet].to_excel(writer, sheet_name = sheet, index = False)
        
#         stats.to_excel(writer, sheet_name = 'Statistics', index = True)
#     new_filename = str(excel_name.split('.')[0]) + '_statistics.xlsx'
#     with pd.ExcelWriter(new_filename) as writer:
#         stats.to_excel(writer, sheet_name = 'Statistics', index = True) 



# def ranked_minimization(sheet, ratio, reject_percentage = 20):

#     mytest = tester[sheet].copy(deep=True)

#     df_mean_before = mytest[ratio].mean()
#     df_1std_before = mytest[ratio].std()
#     df_count_before = mytest[ratio].count()
#     df_2se_perc_before = (2 * mytest[ratio].sem()) / df_mean_before * 100

#     dif_mean = ratio + '_dif_from_mean'
#     dif_1SD = ratio + '_dif_from_1SD'
#     mytest[dif_mean] = mytest.apply(lambda x: abs(x[ratio] - df_mean_before), axis=1)
#     mytest[dif_1SD] = mytest.apply(lambda x: x[dif_mean] - df_1std_before, axis=1)


#     mytest2 = mytest.sort_values(by = dif_1SD, ascending = False)
#     #mytest2.head()

#     ratios_to_reject = int(mytest[ratio].count() * reject_percentage / 100)
#     #print(ratios_to_reject)

#     after_rejection = mytest2[ratios_to_reject:]

#     df_mean_after = after_rejection[ratio].mean()
#     df_1std_after = after_rejection[ratio].std()
#     df_count_after = after_rejection[ratio].count()
    
#     #This is not 2SE%, should probably fix labels... KD 14 June 2021
#     df_2se_perc_after = (2 * after_rejection[ratio].std()) / df_mean_after * 100

#     # print(df_mean_after)
#     # print(df_1std_after)
#     # print(df_2se_perc_after)

#     results_dict = {}
    
#     results_dict['avg_before'] = df_mean_before
#     results_dict['1sd_before'] = df_1std_before
#     results_dict['2se%_before'] = df_2se_perc_before
#     results_dict['avg_after'] = df_mean_after
#     results_dict['1sd_after'] = df_1std_after
#     results_dict['2σ%_after'] = df_2se_perc_after
    
#     return results_dict

# def statistics_ranktest(calc_dict, reject_percentage = 20):
#     calc_list = ['238_CPS', '232_CPS',
#            '208_CPS', '207_CPS', '206_CPS', '204_CPS', '202_CPS', '206/238',
#            '208/232', '207/206', '208/206', '206/204','208/204','207/204', '207/235'  ]
#     mega_dict = {}

#     for sheet in calc_dict:
#         tester = calc_dict[sheet]
#         stats_dict = {}
#         for col in tester:

#             if col in calc_list:
#                 #print(col)
#                 if '/' in col:
#                     key_bf = col + '_before rejection'
#                     key_af = col + '_after rejection'
#                     key_bf_se = col + '_before rejection 2se%'
#                     key_af_se = col + '_after rejection 2σ%'
                    
#                     ranked_dict = ranked_minimization(sheet, col, reject_percentage)
#                     stats_dict[key_bf] = ranked_dict['avg_before']
#                     stats_dict[key_bf_se] = ranked_dict['2se%_before']
#                     stats_dict[key_af] = ranked_dict['avg_after']
#                     #This is not 2SE%, should probably fix labels... KD 14 June 2021
#                     stats_dict[key_af_se] = ranked_dict['2σ%_after']
#                 else:
#                     key = col + '_mean'
#                     df_mean = tester[col].mean()
#                     stats_dict[key] = df_mean
#                     #This is not 2SE%, should probably fix labels... KD 14 June 2021
#                     df_precision = (2 * tester[col].std()) / df_mean * 100
#                     stats_dict[col + '_2σ%'] = df_precision
#             if 'OPZ' in col:
#                  stats_dict[col + '_mean'] = tester[col].mean()
#             if 'SNR' in col:
#                  stats_dict[col + '_mean'] = tester[col].mean()
                    
#         stats_dict['Time (s)'] = tester['Elapsed Time'].max()
        
#         #new_string = sheet.replace('time series data', '')
#         new_string = sheet.split('time')[0].rstrip()
#         mega_dict[new_string] = stats_dict

#     df_1 = pd.DataFrame(mega_dict)
#     df_flip = pd.DataFrame.transpose(df_1)
#     return df_flip


# def stat_grouper(stats_df):
#     ''' Input dataframe that was output from statistics_ranktest2 function.'''
#     headers_to_keep = [
#         '206/238_after rejection',
#         '206/238_after rejection 2se%',
#         '208/232_after rejection',
#         '208/232_after rejection 2se%',
#         '207/206_after rejection',
#         '207/206_after rejection 2se%',
#         '208/206_after rejection',
#         '208/206_after rejection 2se%',
#         '206/204_after rejection',
#         '206/204_after rejection 2se%',
#         '208/204_after rejection',
#         '208/204_after rejection 2se%',
#         '207/204_after rejection',
#         '207/204_after rejection 2se%',
#         '207/235_after rejection',
#         '207/235_after rejection 2se%' 
#     ]

#     headers_to_math = [
#         '206/238_after rejection',   
#         '208/232_after rejection',  
#         '207/206_after rejection',    
#         '208/206_after rejection',    
#         '206/204_after rejection',    
#         '208/204_after rejection',   
#         '207/204_after rejection',
#         '207/235_after rejection',
#     ]


#     short_tester = stats_df[headers_to_keep]

#     samples = short_tester.index.values.tolist()
#     group = None
#     group_names = []
#     start = [0]
#     end = []
#     for idx in range(len(samples)):
#         if group == None:
#             group = samples[idx].split(' 1')[0]
#             group_names.append(group)
#         #print(samples[idx])
#         if group in samples[idx]:
#             pass
#         else:
#             end.append(idx)
#             start.append(idx + 1)
#             name = samples[idx].split(' 1')[0]
#             group = name
#             group_names.append(name)
#     end.append(len(samples))

#     # print('start:', start)
#     # print('end:', end)
#     # print('group names:', group_names)

#     stats_dict = {}
#     for idx in range(len(group_names)):
#         group_dict = {}
#         #print(idx)
#         group_df = short_tester.iloc[start[idx]:end[idx]]

#         for col in headers_to_math:
#             #print(group_df[col])
#             name = col.split('_')[0]
#             df_mean = group_df[col].mean()
#             #print(df_mean)
#             group_dict[name + '_mean'] = df_mean
#             df_2se_perc = (2 * group_df[col].std())/df_mean * 100
#             group_dict[name + '_2SE%'] = df_2se_perc 

#         stats_dict[group_names[idx]] = group_dict

#     return pd.DataFrame(stats_dict)

# #Functions for ratio plots
# ###Still some hardcoded stuff in this...

# def ratio_plot(group_namer, primary_std, choice, plot_dict, published_df):

#     # Making groups of the standards for plotting
#     worksheets = list(plot_dict.keys())
#     group = None
#     std_names = []
#     start = [0]
#     end = []
#     for idx in range(len(worksheets)):
#         if group == None:
#             group = worksheets[idx].split('_')[0]
#             std_names.append(group)
#         #print(samples[idx])
#         if worksheets[idx].split('_')[0] in group:
#             pass
#         else:
#             end.append(idx)
#             start.append(idx) 
#             name = worksheets[idx].split('_')[0]
#             group = name
#             std_names.append(name)
#     end.append(len(worksheets))

#     #print('start:', start)
#     #print('end:', end)
#     #print('std names:', std_names)

#     std_dict = {}

#     for idx in range(len(std_names)):
#         group_dict = {}
#         std_group = worksheets[start[idx]:end[idx]]
#         #print('std_group', std_group)
#         group_dict['0%'] = std_group[0]
#         group_dict['20%'] = std_group[2]
#         group_dict['40%'] = std_group[4]

#         group_dict['0% group'] = std_group[1]
#         group_dict['20% group'] = std_group[3]
#         group_dict['40% group'] = std_group[5]

#         std_dict[std_names[idx]] = group_dict


#     #Creating loop for all entries in plot_dict

#     plot_dict_grouped_dict = {}

#     for key in plot_dict:

#         if 'full' in key:

#         #Create groups for plotting

#             samples = plot_dict[key].index.values.tolist()
#             group = None
#             group_names = []
#             start = [0]
#             end = []
#             for idx in range(len(samples)):
#                 if group == None:
#                     group = samples[idx].split(' 1')[0]
#                     group_names.append(group)
#                 #print(samples[idx])
#                 if group in samples[idx]:
#                     pass
#                 else:
#                     end.append(idx)
#                     start.append(idx) # + 1 was original
#                     name = samples[idx].split(' 1')[0]
#                     group = name
#                     group_names.append(name)
#             end.append(len(samples))

#             #print('start:', start)
#             #print('end:', end)
#             #print('group names:', group_names)

#             group_df_list = []
#             group_df_dict = {}

#             for idx in range(len(group_names)):
#                 group_dict = {}
#                 #print(idx)
#                 group_df = plot_dict[key].iloc[start[idx]:end[idx]]
#                 #print(group_df.index.values.tolist())

#                 group_df_list.append(group_df)
#                 group_df_dict[group_names[idx]] = group_df

#             plot_dict_grouped_dict[key] = group_df_dict

#             if 'group' in key:
#                  plot_dict_grouped_dict[key] = plot_dict[key]


#     #Data to be plotted


# #     group_namer = 'Mudtank'
# #     primary_std = 'ttn'

#     primary_std_name = {
#         'glass': "SRM NIST 610",
#         'glass612': "SRM NIST 612",
#         'glassBHVO': "BHVO-2G",
#         'ttn': "MKED-1",
#     }

#     test_zero = plot_dict_grouped_dict[ std_dict[primary_std]['0%']][group_namer]
#     test_20 = plot_dict_grouped_dict[ std_dict[primary_std]['20%']][group_namer]
#     test_40 = plot_dict_grouped_dict[ std_dict[primary_std]['40%']][group_namer]

#     test_group_zero = plot_dict[std_dict[primary_std]['0% group']].loc[group_namer]
#     test_group_20 = plot_dict[std_dict[primary_std]['20% group']].loc[group_namer]
#     test_group_40 = plot_dict[std_dict[primary_std]['40% group']].loc[group_namer]


#     if (group_namer == 'Bancroft') or (group_namer =='Yates'):
#         #print('debug')
#         pub_ratio = published_df.loc[['OLT-1', 'OLT-2', 'TCB']]
#     else:
#         pub_ratio = published_df.loc[group_namer]


#     color_dict = {
#         '0%': 'red',
#         '20%': 'blue',
#         '40%': 'green',
#         'pub': 'purple'
#     }


#     #Initialize plots
#     fig, axs = plt.subplots(2, 2, sharex = False, figsize = (24, 12))
#     fig.suptitle(group_namer + ', Primary Std: ' + primary_std_name[primary_std], fontsize=28)

#     #207/204 vs 206/204 plot
#         #Results for each line
#     x1_entry = '206/204 corrected'
#     y1_entry = '207/204 corrected'
#     ax = axs[0,0]
#     plt.setp(ax.get_xticklabels(), fontsize=16)
#     plt.setp(ax.get_yticklabels(), fontsize=16)

#     ax.plot(test_zero[x1_entry], test_zero[y1_entry], 'o', color= color_dict['0%']) #0% reject
#     ax.plot(test_20[x1_entry], test_20[y1_entry], 'o', color= color_dict['20%']) #20% reject
#     ax.plot(test_40[x1_entry], test_40[y1_entry], 'o', color= color_dict['40%']) #40% reject
#         #Group averages
#     x1_entry = '206/204 corrected_mean'
#     y1_entry = '207/204 corrected_mean'
#     ax.plot(test_group_zero[x1_entry], test_group_zero[y1_entry], '^', color= color_dict['0%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject
#     ax.plot(test_group_20[x1_entry], test_group_20[y1_entry], '^', color= color_dict['20%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject
#     ax.plot(test_group_40[x1_entry], test_group_40[y1_entry], '^', color= color_dict['40%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject
#         #Group 2sigma
#     x_error = '206/204 corrected_2σ'
#     y_error = '206/204 corrected_2σ'
#     ax.errorbar(test_group_zero[x1_entry], test_group_zero[y1_entry],test_group_zero[x_error], test_group_zero[y_error],color = 'black')
#     ax.errorbar(test_group_20[x1_entry], test_group_20[y1_entry],test_group_20[x_error], test_group_20[y_error],color = 'black')
#     ax.errorbar(test_group_40[x1_entry], test_group_40[y1_entry],test_group_40[x_error], test_group_40[y_error],color = 'black')
#         #Published value
#     x1_entry = '206/204'
#     y1_entry = '207/204'
#     ax.plot(pub_ratio[x1_entry], pub_ratio[y1_entry], 's', color= color_dict['pub'], markersize=18, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject

#         #Formatting
#     ax.set(xlabel = '206/204', ylabel = '207/204',)
#     ax.yaxis.label.set_size(20)
#     ax.xaxis.label.set_size(20)


#     #208/204 vs 206/204 plot

#     x2_entry = '206/204 corrected'
#     y2_entry = '208/204 corrected'
#     ax = axs[0,1]
#     plt.setp(ax.get_xticklabels(), fontsize=16)
#     plt.setp(ax.get_yticklabels(), fontsize=16)

#     ax.plot(test_zero[x2_entry], test_zero[y2_entry], 'o', color= color_dict['0%']) #0% reject
#     ax.plot(test_20[x2_entry], test_20[y2_entry], 'o', color= color_dict['20%']) #20% reject
#     ax.plot(test_40[x2_entry], test_40[y2_entry], 'o', color= color_dict['40%']) #40% reject

#     x1_entry = '206/204 corrected_mean'
#     y1_entry = '208/204 corrected_mean'
#     ax.plot(test_group_zero[x1_entry], test_group_zero[y1_entry], '^', color= color_dict['0%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject
#     ax.plot(test_group_20[x1_entry], test_group_20[y1_entry], '^', color= color_dict['20%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject
#     ax.plot(test_group_40[x1_entry], test_group_40[y1_entry], '^', color= color_dict['40%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject

#      #Group 2sigma
#     x_error = '206/204 corrected_2σ'
#     y_error = '208/204 corrected_2σ'
#     ax.errorbar(test_group_zero[x1_entry], test_group_zero[y1_entry],test_group_zero[x_error], test_group_zero[y_error],color = 'black')
#     ax.errorbar(test_group_20[x1_entry], test_group_20[y1_entry],test_group_20[x_error], test_group_20[y_error],color = 'black')
#     ax.errorbar(test_group_40[x1_entry], test_group_40[y1_entry],test_group_40[x_error], test_group_40[y_error],color = 'black')

#      #Published value
#     x1_entry = '206/204'
#     y1_entry = '208/204'
#     ax.plot(pub_ratio[x1_entry], pub_ratio[y1_entry], 's', color= color_dict['pub'], markersize=18, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject

#     #Formatting
#     ax.set(xlabel = '206/204', ylabel = '208/204',)
#     ax.yaxis.label.set_size(20)
#     ax.xaxis.label.set_size(20)



#     #207/206 vs 208/206 plot

#     x3_entry = '208/206 corrected'
#     y3_entry = '207/206 corrected'
#     ax = axs[1,0]
#     plt.setp(ax.get_xticklabels(), fontsize=16)
#     plt.setp(ax.get_yticklabels(), fontsize=16)


#     custom_lines = ax.plot(test_zero[x3_entry], test_zero[y3_entry], 'o', color= color_dict['0%']) #0% reject
#     ax.plot(test_20[x3_entry], test_20[y3_entry], 'o', color= color_dict['20%']) #20% reject
#     ax.plot(test_40[x3_entry], test_40[y3_entry], 'o', color= color_dict['40%']) #40% reject

#     x1_entry = '208/206 corrected_mean'
#     y1_entry = '207/206 corrected_mean'
#     ax.plot(test_group_zero[x1_entry], test_group_zero[y1_entry], '^', color= color_dict['0%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject
#     ax.plot(test_group_20[x1_entry], test_group_20[y1_entry], '^', color= color_dict['20%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject
#     ax.plot(test_group_40[x1_entry], test_group_40[y1_entry], '^', color= color_dict['40%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject

#      #Group 2sigma
#     x_error = '208/206 corrected_2σ'
#     y_error = '207/206 corrected_2σ'
#     ax.errorbar(test_group_zero[x1_entry], test_group_zero[y1_entry],test_group_zero[x_error], test_group_zero[y_error],color = 'black')
#     ax.errorbar(test_group_20[x1_entry], test_group_20[y1_entry],test_group_20[x_error], test_group_20[y_error],color = 'black')
#     ax.errorbar(test_group_40[x1_entry], test_group_40[y1_entry],test_group_40[x_error], test_group_40[y_error],color = 'black')

#      #Published value
#     x1_entry = '208/206'
#     y1_entry = '207/206'
#     ax.plot(pub_ratio[x1_entry], pub_ratio[y1_entry], 's', color= color_dict['pub'], markersize=18, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject

#     #Formatting
#     ax.set(xlabel = '208/206', ylabel = '207/206',)
#     ax.yaxis.label.set_size(20)
#     ax.xaxis.label.set_size(20)


#     #206/238 vs 208/232 plot

#     x4_entry = '208/232 corrected'
#     y4_entry = '206/238 corrected'
#     ax = axs[1,1]
#     plt.setp(ax.get_xticklabels(), fontsize=16)
#     plt.setp(ax.get_yticklabels(), fontsize=16)

#     ax.plot(test_zero[x4_entry], test_zero[y4_entry], 'o', color= color_dict['0%']) #0% reject
#     ax.plot(test_20[x4_entry], test_20[y4_entry], 'o', color= color_dict['20%']) #20% reject
#     ax.plot(test_40[x4_entry], test_40[y4_entry], 'o', color= color_dict['40%']) #40% reject

#     x1_entry = '208/232 corrected_mean'
#     y1_entry = '206/238 corrected_mean'
#     ax.plot(test_group_zero[x1_entry], test_group_zero[y1_entry], '^', color= color_dict['0%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject
#     ax.plot(test_group_20[x1_entry], test_group_20[y1_entry], '^', color= color_dict['20%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject
#     ax.plot(test_group_40[x1_entry], test_group_40[y1_entry], '^', color= color_dict['40%'], markersize=14, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject

#      #Group 2sigma
#     x_error = '208/232 corrected_2σ'
#     y_error = '206/238 corrected_2σ'
#     ax.errorbar(test_group_zero[x1_entry], test_group_zero[y1_entry],test_group_zero[x_error], test_group_zero[y_error],color = 'black')
#     ax.errorbar(test_group_20[x1_entry], test_group_20[y1_entry],test_group_20[x_error], test_group_20[y_error],color = 'black')
#     ax.errorbar(test_group_40[x1_entry], test_group_40[y1_entry],test_group_40[x_error], test_group_40[y_error],color = 'black')

#      #Published value
#     x1_entry = '208/232'
#     y1_entry = '206/238'
#     ax.plot(pub_ratio[x1_entry], pub_ratio[y1_entry], 's', color= color_dict['pub'], markersize=18, markeredgewidth=1.5, markeredgecolor= 'black') #0% reject

#     #Formatting
#     ax.set(xlabel = '208/232', ylabel = '206/238',)
#     ax.yaxis.label.set_size(20)
#     ax.xaxis.label.set_size(20)

#     # Make a legend

#     patch_0 = mpatches.Patch(color= color_dict['0%'], label= '0% Rejection')
#     patch_20 = mpatches.Patch(color=color_dict['20%'], label= '20% Rejection')
#     patch_40 = mpatches.Patch(color=color_dict['40%'], label= '40% Rejection')
#     patch_pub = mpatches.Patch(color=color_dict['pub'], label= 'Published Data')
#     plt.legend(handles=[patch_0, patch_20, patch_40, patch_pub])


#     MYDIR = ("Ratio_Figures")
#     CHECK_FOLDER = os.path.isdir(MYDIR)

#     # If folder doesn't exist, then create it.
#     if not CHECK_FOLDER:
#         os.makedirs(MYDIR)
#         #print("created folder : ", MYDIR)
    
#     new_string = group_namer + '_PrimaryStd_' + primary_std_name[primary_std]
    
#     filename = os.path.join(MYDIR, new_string + '.pdf')
#     plt.savefig(filename)
#     print('Plot for ', new_string, " is complete.")
    
#     if choice == False:
#         plt.close()

        
        
# #ratio_plot(group_namer, primary_std, choice = True, plot_dict = plot_dict, published_df = published_df)
# #ratio_report(plot_dict = plot_dict, intro_filename = 'fake.pdf', intro = False, output_name = 'ratio_output.pdf'):
    
# def ratio_report(plot_dict, intro_filename, intro, output_name):
    
#     MYDIR = ("Ratio_Figures")
#     mergedObject = PdfFileMerger()
    
#     if intro:
#         mergedObject.append(PdfFileReader(intro_filename, 'rb'))
#         print(f'Succesfully incorporated {intro_filename} into PDF.')
    
    
#     #Actual loop of plotting
    
#     # Making groups of the standards for plotting
#     worksheets = list(plot_dict.keys())
#     group = None
#     std_names = []
#     start = [0]
#     end = []
#     for idx in range(len(worksheets)):
#         if group == None:
#             group = worksheets[idx].split('_')[0]
#             std_names.append(group)
#         #print(samples[idx])
#         if worksheets[idx].split('_')[0] in group:
#             pass
#         else:
#             end.append(idx)
#             start.append(idx) 
#             name = worksheets[idx].split('_')[0]
#             group = name
#             std_names.append(name)
#     end.append(len(worksheets))

#     #print('start:', start)
#     #print('end:', end)
#     #print('std names:', std_names)

#     group_list = None

#     for key in plot_dict:

#         if 'full' in key:

#         #Create groups for plotting

#             samples = plot_dict[key].index.values.tolist()
#             group = None
#             group_names = []
#             start = [0]
#             end = []
#             for idx in range(len(samples)):
#                 if group == None:
#                     group = samples[idx].split(' 1')[0]
#                     group_names.append(group)
#                 #print(samples[idx])
#                 if group in samples[idx]:
#                     pass
#                 else:
#                     end.append(idx)
#                     start.append(idx) # + 1 was original
#                     name = samples[idx].split(' 1')[0]
#                     group = name
#                     group_names.append(name)
#             end.append(len(samples))

#             #print('start:', start)
#             #print('end:', end)
#             #print('group names:', group_names)
#             group_list = group_names

#     #print('group names:', group_list)

#     primary_std_name = {
#         'glass': "SRM NIST 610",
#         'glass612': "SRM NIST 612",
#         'glassBHVO': "BHVO-2G",
#         'ttn': "MKED-1",
#     }
    
# #ratio_plot(group_namer, primary_std, choice = True, plot_dict = plot_dict, published_df = published_df)    
#     for std in std_names:
#         #print(std)

#         for group in group_list:
#             if 'apa' in group:
#                 print(group +' was skipped.')
#                 continue
#             else:
#                 #print(group)
#                 ratio_plot(group, std, False, plot_dict, published_df)
#                 new_string = group + '_PrimaryStd_' + primary_std_name[std]

#                 filename = os.path.join(MYDIR, new_string + '.pdf')
#                 mergedObject.append(PdfFileReader(filename, 'rb'))
    
#     if '.pdf' in output_name:
#         pass
#     else:
#         output_name = output_name + '.pdf'
    
     
#     mergedObject.write(output_name)

#     print(f'PDF file named: {output_name} is complete.')   


# def AB_error2(result_df):
#     ''' Unfinished (and maybe broken) function toward calculating the Constant Ext Error from Isoplot?'''
#     pd.set_option('mode.chained_assignment',None)
    
#     AB_err_tester = result_df
#     result_headers = [
#         '206/238 corrected',
#         '208/232 corrected',
#         '207/206 corrected',
#         '208/206 corrected',
#         '206/204 corrected',
#         '207/204 corrected',
#         '208/204 corrected'    
#     ]

#     for ratio in result_headers:
#         #print(ratio)
#         error_string = str(ratio.split()[0]) + '_after rejection 2σ%'
#         new_string1 = ratio + ' 2σ'
#         #print(error_string)
#         new_string2 = ratio + ' x/(σ^2)'
#         new_string3 = ratio + ' 1/(σ^2)'

#         AB_err_tester[new_string1] = AB_err_tester.apply(lambda x: (x[error_string] / 100 ) * (x[ratio]), axis=1)

#         AB_err_tester[new_string2] = AB_err_tester.apply(lambda x: x[ratio] / ((x[new_string1] / 2)**2), axis=1)
#         AB_err_tester[new_string3] = AB_err_tester.apply(lambda x: (x[new_string2] / x[ratio]), axis=1)


#     samples = AB_err_tester.index.values.tolist()
#     group = None
#     group_names = []
#     start = [0]
#     end = []
#     for idx in range(len(samples)):
#         if group == None:
#             group = samples[idx].split(' 1')[0]
#             group_names.append(group)
#         #print(samples[idx])
#         if group in samples[idx]:
#             pass
#         else:
#             end.append(idx)
#             start.append(idx) # + 1 was original
#             name = samples[idx].split(' 1')[0]
#             group = name
#             group_names.append(name)
#     end.append(len(samples))

#     # print('start:', start)
#     # print('end:', end)
#     # print('group names:', group_names)

#     stats_dict = {}
#     corrected_stats_dict = {}
#     group_df_list = []
    
#     for idx in range(len(group_names)):
#         group_dict = {}
#         #print(idx)
#         group_df = AB_err_tester.iloc[start[idx]:end[idx]]
#         #print(group_df.index.values.tolist())
#         for ratio in result_headers:
#             new_string2 = ratio + ' x/(σ^2)'
#             new_string3 = ratio + ' 1/(σ^2)'
#             new_string4 = ratio + ' Weighted Mean'
#             new_string5 = ratio + ' σ^2'
#             new_string6 = ratio + ' σ'
#             wtd_mean = group_df[new_string2].sum() / group_df[new_string3].sum()
#             sigma_sqrd = 1 / group_df[new_string3].sum()
#             sigma = sigma_sqrd ** 0.5

#             group_dict[new_string4] = wtd_mean
#             group_dict[new_string5] = sigma_sqrd
#             group_dict[new_string6] = sigma
            
#             z_string = str(ratio.split()[0]) + '_Z'
#             group_df[z_string] = group_df.apply(lambda x: (((x[ratio] - group_dict[new_string4])**2) * x[new_string3] ), axis=1)
            
#             new_string5b = ratio +' S-factor'
#             #new_string6 = ratio + ' count'
#             new_string7b = ratio + ' MSWD'
#             #new_string8 = ratio + ' MSWD_difference_from_1'

#             S_factor = group_df[z_string].sum()
#             MSWD = S_factor / (group_df[ratio].count() - 1)
#             #MSWD_dif = 1 - MSWD

#             group_dict[new_string5b] = S_factor
#             group_dict[new_string7b] = MSWD
        
#         group_df_list.append(group_df)

#         stats_dict[group_names[idx]] = group_dict

#     grouper = pd.DataFrame(stats_dict) ### Will need to flip this at some point.
#     new_AB_err_tester = pd.concat(group_df_list)
    
#     grouper_flip = pd.DataFrame.transpose(grouper)

#     return new_AB_err_tester, grouper_flip 


StarcoderdataJupyterScriptsDedupFiltered