﻿#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Sat Jul 14 12:53:59 2018 @author: Max Lloyd, UC Berkeley, mklloyd@berkeley.edu Script to calculate d13C, dD, ∆13CH3D, and ∆12CH2D2 values from IRMS measurements relative to a working reference gas of known composition. Error propagation is done using a Monte Carlo scheme and the 'true' random number generator of the local processor. """ import numpy as np import pandas as pd import random import re def calculate_randoms(row, mean, std_dev): return(row['seed'].gauss(mean, std_dev)) def CH4_bulk_comp(dD_wg, d13C_wg, d13CD_wg=np.nan, dD2_wg=np.nan, frag_rate=0.8109): ''' Solves the linear equations to accurately calculate the bulk composition of d13C and dD. ''' r13C_vpdb = 0.0112372 rD_vsmow = 0.00015576 # wg composition based on UC Davis analyses (Eldridge et al.) d13C_ref = -38.3672 dD_ref = -159.2795 # Calculate true R13 and RD wg ratios r13C_ref = (d13C_ref/1000 + 1)*r13C_vpdb rD_ref = (dD_ref/1000 + 1)*rD_vsmow # Use these to calculate measured sample RD and R13C ratios m12CH3D_sa = (dD_wg/1000 + 1)*rD_ref*4/(1 + frag_rate*3*rD_ref) rD_sa = m12CH3D_sa/(4 - 3*m12CH3D_sa*frag_rate) m13CH4_sa = (d13C_wg/1000 + 1)*r13C_ref/(1+frag_rate*(r13C_ref + 3*rD_ref)) r13C_sa = m13CH4_sa*(1 + 3*frag_rate*rD_sa)/(1 - m13CH4_sa*frag_rate) # Use wg clumped isotope compositions from # Eldridge et al. ACS Earth and Space Chemistry D13CD_ref = 2.5884 r13CD_ref_stoch = 4*r13C_ref*rD_ref r13CD_ref = (D13CD_ref/1000 + 1)*r13CD_ref_stoch DD2_ref = 5.8600 rD2_ref_stoch = 6*rD_ref*rD_ref rD2_ref = (DD2_ref/1000 + 1)*rD2_ref_stoch # Calculate sample clumped isotope compositions in stochastic ref frame r13CD_sa = ((d13CD_wg/1000 + 1)*r13CD_ref/(1 + frag_rate*( r13C_ref + 3*rD_ref)))*(1 + frag_rate*(r13C_sa + 3*rD_sa)) rD2_sa = ((dD2_wg/1000 + 1)*rD2_ref/(1 + frag_rate*( r13C_ref + 3*rD_ref)))*(1 + frag_rate*(r13C_sa + 3*rD_sa)) D13CD_stoch = (r13CD_sa/(4*r13C_sa*rD_sa)-1)*1000 DD2_stoch = (rD2_sa/(6*rD_sa*rD_sa)-1)*1000 dD_vsmow = (rD_sa/rD_vsmow-1)*1000 d13C_vpdb = (r13C_sa/r13C_vpdb-1)*1000 return(pd.DataFrame({'dD_vsmow': dD_vsmow, 'd13C_vpdb': d13C_vpdb, 'D13CD_wg': D13CD_stoch, 'DD2_wg': DD2_stoch})) if __name__ == '__main__': file_path = ('C:/Users/Administrator/Desktop/Data/summary_sheets/' 'methane_summary_sheet.xlsx') # read summary sheet d = pd.read_excel(file_path, header=[0, 1], index_col=0) # get indices of samples that are unprocessed i_to_process = d.loc[(d[('Calculated data'), 'D13CD_wg'].isnull()) | ( (d[('dD2')].notnull().all(axis=1)) & ( d[('Calculated data'), 'DD2_wg'].isnull()))].index if len(i_to_process) == 0: print('No columns to process. Goodbye! ') else: print('Processing rows: {0}'.format([j for j in i_to_process])) N = int(1e5) # create a new database for MC simulations with random seeds dmc = pd.DataFrame(data={ 'i': np.arange(N), 'seed': [random.SystemRandom() for j in range(N)]}) cols_to_make = ['dD_vsmow', 'd13C_vpdb', 'D13CD_wg', 'DD2_wg'] for i in i_to_process: # perturb each measurement according a normal distribution print('MC error propagation for row {0}... '.format(i)) dmc['dD_wg'] = dmc.apply(calculate_randoms, axis=1, args=( d.loc[i, 'dD']['mean vs. wg'], d.loc[i, 'dD']['std error'])) dmc['d13C_wg'] = dmc.apply(calculate_randoms, axis=1, args=( d.loc[i, 'd13C']['mean vs. wg'], d.loc[i, 'd13C']['std error'])) dmc['d13CD_wg'] = dmc.apply(calculate_randoms, axis=1, args=( d.loc[i, 'd13CD']['mean vs. wg'], d.loc[i, 'd13CD']['std error'])) dmc['dD2_wg'] = dmc.apply(calculate_randoms, axis=1, args=( d.loc[i, 'dD2']['mean vs. wg'], d.loc[i, 'dD2']['std error'])) # process each row as a separate measurement results = CH4_bulk_comp(dmc['dD_wg'], dmc['d13C_wg'], d13CD_wg=dmc['d13CD_wg'], dD2_wg=dmc['dD2_wg'], frag_rate=d.loc[i, 'd13C'] ['fragmentation rate']) # calculate the mean of all rows as the true measurement mean d.loc[i, ('Calculated data', cols_to_make)] = results[ cols_to_make].mean().round(3).values # calculate the std. dev. of all rows to estimate uncertainties cols_se = [re.sub(r'_\w+', '_std_error', i) for i in cols_to_make] d.loc[i, ('Calculated data', cols_se)] = results[ cols_to_make].std().round(3).values print('MC calcs complete. Writing new spreadsheet... ') # Export data d.loc[:, (['Info', 'Calculated data'])].to_excel( ('C:/Users/Administrator/Desktop/Data/summary_sheets/' 'methane_summary_sheet_calc.xlsx'))