He_date_calc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# (U-Th)/He data reduction notebook for the HAL at UIUC\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This notebook was written by William Guenthner in fall of 2022 for the reduction of grain size data, U, Th, Sm, and He measurements towards calculation of (U-Th)/He dates. Some of the inputs and file formats are specific to data generated in the Helium Analysis Laboratory (HAL) at the University of Illinois Urbana-Champaign (UIUC), but hopefully it has broader applicability and utilty for other lab groups. \n",
    "\n",
    "The notebook is structured to interact with 3 separate CSV files that should be colocated with each instance of the notebook in the same folder. The CSV files are related to: 1) U, Th, Sm, Zr, and Ca amount measurements obtained from ICP-MS anlaysis (obtained with an iCAP Q using Qtegra software at UIUC), 2) He amount measurements (obtained with a PrismaPlus 220 and reported as peak hops on masses 1-5 at UIUC), and 3) grain size measurements for Ft correction. Cells are grouped below roughly in that order of reduction (wet chemistry first, then He, then grain size)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "project = ''\n",
    "bomb_run = 1\n",
    "\n",
    "#packages to import and constants\n",
    "import csv\n",
    "import json\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "display(HTML(\"<style>div.output_scroll {height: 44em;}</style>\"))\n",
    "pd.options.display.max_rows = 100\n",
    "pd.options.display.max_columns = 25\n",
    "\n",
    "#set file path, uncomment based on directory structure applicable to your work\n",
    "#the first option is for the current directory, others are more verbose paths specific to HAL computers\n",
    "\n",
    "file_path = os.getcwd()\n",
    "#file_path = os.path.join('C:\\\\','Users','lab-admin','Documents','GitHub','HAL_data','projects_2023',project)\n",
    "#file_path = os.path.join('/Users','wrg','Documents','GitHub','HAL_data','projects_2023',project)\n",
    "#file_path = os.path.join('/Users','wrg','Documents','GitHub','HAL_data','sandbox')\n",
    "os.chdir(file_path)\n",
    "\n",
    "Avogadro = 6.022045e23 #atom/mol\n",
    "ideal_gas_moles = 22.414 #liter/mol\n",
    "\n",
    "#all in amu\n",
    "mass_233U = 233.0396280 \n",
    "mass_234U = 234.0409456\n",
    "mass_235U = 235.0439231\n",
    "mass_236U = 236.0455619\n",
    "mass_238U = 238.0507826\n",
    "mass_U = 238.028913\n",
    "\n",
    "mass_229Th = 229.031754\n",
    "mass_230Th = 230.033127\n",
    "mass_232Th = 232.0380504\n",
    "mass_Th = 232.0380504\n",
    "\n",
    "mass_144Sm = 143.911995\n",
    "mass_147Sm = 146.914893\n",
    "mass_148Sm = 147.914818\n",
    "mass_149Sm = 148.917180\n",
    "mass_150Sm = 149.917271\n",
    "mass_152Sm = 151.919728\n",
    "mass_154Sm = 153.922205\n",
    "mass_Sm = 150.366344\n",
    "\n",
    "mass_40Ca = 39.96259\n",
    "mass_42Ca = 41.95862\n",
    "mass_43Ca = 42.95877\n",
    "mass_44Ca = 43.95548\n",
    "mass_46Ca = 45.95369\n",
    "mass_48Ca = 47.95243\n",
    "mass_Ca = 40.08601\n",
    "\n",
    "mass_90Zr = 89.90470\n",
    "mass_91Zr = 90.90564\n",
    "mass_92Zr = 91.90504\n",
    "mass_94Zr = 93.90631\n",
    "mass_96Zr = 95.90828\n",
    "mass_Zr = 91.22365\n",
    "\n",
    "#all in 1/yr\n",
    "lambda_238 = 1.55125e-10\n",
    "lambda_235 = 9.84850e-10\n",
    "lambda_232 = 4.9475e-11\n",
    "lambda_147 = 6.54e-12\n",
    "\n",
    "#Ft constants a1 and a2 from Farley (2002) Table 1\n",
    "a1_238U_ap = -5.13\n",
    "a2_238U_ap = 6.78\n",
    "a1_Th_ap = -5.90\n",
    "a2_Th_ap = 8.99\n",
    "a1_235U_ap = a1_Th_ap\n",
    "a2_235U_ap = a2_Th_ap\n",
    "\n",
    "a1_238U_zirc = -4.31\n",
    "a2_238U_zirc = 4.92\n",
    "a1_Th_zirc = -5.00\n",
    "a2_Th_zirc = 6.80\n",
    "a1_235U_zirc = a1_Th_zirc\n",
    "a2_235U_zirc = a2_Th_zirc\n",
    "\n",
    "#g/cc\n",
    "ap_density = 3.19\n",
    "zirc_density = 4.60\n",
    "\n",
    "#chemistry constants specific to HAL\n",
    "#mL\n",
    "Vnm_UTh = 0.025\n",
    "d_Vnm_UTh = Vnm_UTh * 0.01\n",
    "Vnm_Sm = 0.025\n",
    "d_Vnm_Sm = Vnm_Sm * 0.01\n",
    "\n",
    "#abudances in normal solutions\n",
    "Ab_238U_nm = 1 - (1/137.818)\n",
    "d_Ab_238U_nm = 0\n",
    "Ab_235U_nm = 1/137.818\n",
    "d_Ab_235U_nm = 0\n",
    "Ab_147Sm_nm = 0.1499\n",
    "d_Ab_147Sm_nm = 0\n",
    "Ab_149Sm_nm = 0.1382\n",
    "d_Ab_149Sm_nm = 0\n",
    "Ab_152Sm_nm = 0.2675\n",
    "d_Ab_152Sm_nm = 0\n",
    "\n",
    "#ng/mL\n",
    "concnm_U = 25.3458\n",
    "d_concnm_U = 0.0181\n",
    "concnm_Th = 49.8397\n",
    "d_concnm_Th = 0.0284\n",
    "concnm_Sm = 50.0797\n",
    "d_concnm_Sm = 0.0284"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## U, Th, Sm, Zr, and Ca reduction\n",
    "\n",
    "At UIUC, we measure ratios of $^{238}$U/$^{236}$U, $^{232}$Th/$^{230}$Th, $^{152}$Sm/$^{149}$Sm, $^{90}$Zr/$^{91}$Zr, and $^{40}$Ca/$^{42}$Ca in our unknowns. For Sm, additional ratios are measured in our spike normals. The notebook is therefore designed around those specific ratios. It is also designed for the particular output format and column headers reported by an iCAP Q ICP-MS running the Qtegra software\n",
    "\n",
    "First, we open our Qtegra CSV file and extract the relevant ratios. The U_Th_file needs to be updated with the appropriate name used for the iCAP run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#csv file for U Th data\n",
    "U_Th_file = project + '_' + str(bomb_run) + '_UTh_data.csv'\n",
    "\n",
    "#function to extract relevant columns from csv file\n",
    "def get_csv_data(file_name, type_string, title_string):\n",
    "    \n",
    "    with open(file_name, 'r', newline = '') as in_file:\n",
    "        csv_reader_UTh = csv.reader(in_file)\n",
    "    \n",
    "        #get type of columns (Raw.Average, Raw.Ratio.STD, etc) and title (238/236 (KED), etc.)\n",
    "        run_type = next(csv_reader_UTh) \n",
    "        next(csv_reader_UTh) #throw away empty row\n",
    "        run_title = next(csv_reader_UTh)\n",
    "        next(csv_reader_UTh) #throw away this row too\n",
    "        \n",
    "        col_num = 0\n",
    "        while col_num < len(run_type) and (run_type[col_num] != type_string or run_title[col_num] != title_string):\n",
    "            col_num = col_num + 1\n",
    "        if col_num == len(run_type):\n",
    "            print('Error: column type or title not found')\n",
    "        else:\n",
    "             list_col = [0 if line[col_num]=='N/A' or line[col_num]=='' else float(line[col_num]) for line in csv_reader_UTh]\n",
    "                \n",
    "        return list_col\n",
    "\n",
    "#extract relevant columns from csv file\n",
    "list_149intensity = get_csv_data(U_Th_file, 'Raw.Average', '149Sm (KED)')\n",
    "list_152intensity = get_csv_data(U_Th_file, 'Raw.Average', '152Sm (KED)')\n",
    "list_230intensity = get_csv_data(U_Th_file, 'Raw.Average', '230Th (KED)')\n",
    "list_232intensity = get_csv_data(U_Th_file, 'Raw.Average', '232Th (KED)')\n",
    "list_236intensity = get_csv_data(U_Th_file, 'Raw.Average', '236U (KED)')\n",
    "list_236intensity = get_csv_data(U_Th_file, 'Raw.Average', '236U (KED)')\n",
    "list_238intensity = get_csv_data(U_Th_file, 'Raw.Average', '238U (KED)')\n",
    "list_238_236 = get_csv_data(U_Th_file, 'Raw.Ratio.Average', '238U (KED) / 236U (KED)')\n",
    "list_238_236_1s = get_csv_data(U_Th_file, 'Raw.Ratio.STD', '238U (KED) / 236U (KED)')\n",
    "list_232_230 = get_csv_data(U_Th_file, 'Raw.Ratio.Average', '232Th (KED) / 230Th (KED)')\n",
    "list_232_230_1s = get_csv_data(U_Th_file, 'Raw.Ratio.STD', '232Th (KED) / 230Th (KED)')\n",
    "list_152_149 = get_csv_data(U_Th_file, 'Raw.Ratio.Average', '152Sm (KED) / 149Sm (KED)')\n",
    "list_152_149_1s = get_csv_data(U_Th_file, 'Raw.Ratio.STD', '152Sm (KED) / 149Sm (KED)')\n",
    "list_144_152 = get_csv_data(U_Th_file, 'Raw.Ratio.Average', '144Sm (KED) / 152Sm (KED)')\n",
    "list_144_152_1s = get_csv_data(U_Th_file, 'Raw.Ratio.STD', '144Sm (KED) / 152Sm (KED)')\n",
    "list_147_152 = get_csv_data(U_Th_file, 'Raw.Ratio.Average', '147Sm (KED) / 152Sm (KED)')\n",
    "list_147_152_1s = get_csv_data(U_Th_file, 'Raw.Ratio.STD', '147Sm (KED) / 152Sm (KED)')\n",
    "list_148_152 = get_csv_data(U_Th_file, 'Raw.Ratio.Average', '148Sm (KED) / 152Sm (KED)')\n",
    "list_148_152_1s = get_csv_data(U_Th_file, 'Raw.Ratio.STD', '148Sm (KED) / 152Sm (KED)')\n",
    "list_150_152 = get_csv_data(U_Th_file, 'Raw.Ratio.Average', '150Sm (KED) / 152Sm (KED)')\n",
    "list_150_152_1s = get_csv_data(U_Th_file, 'Raw.Ratio.STD', '150Sm (KED) / 152Sm (KED)')\n",
    "list_154_152 = get_csv_data(U_Th_file, 'Raw.Ratio.Average', '154Sm (KED) / 152Sm (KED)')\n",
    "list_154_152_1s = get_csv_data(U_Th_file, 'Raw.Ratio.STD', '154Sm (KED) / 152Sm (KED)')\n",
    "\n",
    "\n",
    "#extract sample names from csv file\n",
    "with open(U_Th_file, 'r') as in_file:\n",
    "    col_num = 1\n",
    "    csv_reader_UTh = csv.reader(in_file)\n",
    "    for i in range (4):\n",
    "        next(csv_reader_UTh) #throw away the junk rows\n",
    "    aliquot_list = [str(line[col_num]) for line in csv_reader_UTh]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assign different aliquot types (acid blank, spike blank, spike normal, sample, etc.) in the next cell. We also want to check blank intensities and spike normal ratio consistency throughout the run before moving on, which are reported out below the cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#spike normal ratio lists\n",
    "SN_U_vals = []\n",
    "SN_U_err = []\n",
    "SN_Th_vals = [] \n",
    "SN_Th_err = []\n",
    "SN_Sm_vals = []\n",
    "SN_Sm_err = []\n",
    "SN_names = []\n",
    "\n",
    "#spike blank ratio lists\n",
    "SB_U_vals = []\n",
    "SB_U_err = []\n",
    "SB_Th_vals = [] \n",
    "SB_Th_err = []\n",
    "SB_Sm152_149_vals = []\n",
    "SB_Sm152_149_err = []\n",
    "SB_Sm154_152_vals = []\n",
    "SB_Sm154_152_err = []\n",
    "SB_Sm150_152_vals = []\n",
    "SB_Sm150_152_err = []\n",
    "SB_Sm148_152_vals = []\n",
    "SB_Sm148_152_err = []\n",
    "SB_Sm147_152_vals = []\n",
    "SB_Sm147_152_err = []\n",
    "SB_Sm144_152_vals = []\n",
    "SB_Sm144_152_err = []\n",
    "\n",
    "#sample ratio lists\n",
    "sample_U_vals = []\n",
    "sample_U_err = []\n",
    "sample_Th_vals = []\n",
    "sample_Th_err = []\n",
    "sample_Sm_vals = []\n",
    "sample_Sm_err = []\n",
    "sample_names = []\n",
    "\n",
    "#populate lists depending on aliquot type\n",
    "#also report out the blanks values for 152, 232, and 238\n",
    "for i in range(len(aliquot_list)):\n",
    "    if aliquot_list[i].startswith('AB'):\n",
    "        print('Blank levels for acid blank ',aliquot_list[i], ' are:\\n',list_152intensity[i], ' for 152Sm, ',list_232intensity[i],' for 232Th, and ',list_238intensity[i],' for 238U\\n')\n",
    "    elif aliquot_list[i].startswith('Bomb') or aliquot_list[i].startswith('BB'):\n",
    "        print('Blank levels for bomb blank ',aliquot_list[i],' are:\\n',list_152intensity[i],' for 152Sm, ',list_232intensity[i],' for 232Th, and ',list_238intensity[i],' for 238U\\n')\n",
    "    elif aliquot_list[i].startswith('Empty') or aliquot_list[i].startswith('Blank'):\n",
    "        print('Blank levels for Nb blank ',aliquot_list[i],' are:\\n',list_152intensity[i],' for 152Sm, ',list_232intensity[i],' for 232Th, and ',list_238intensity[i],' for 238U\\n')\n",
    "    elif aliquot_list[i].startswith('SN'):\n",
    "        SN_names.append(aliquot_list[i])\n",
    "        SN_U_vals.append(float(list_238_236[i]))\n",
    "        SN_U_err.append(float(list_238_236_1s[i]))\n",
    "        SN_Th_vals.append(float(list_232_230[i]))\n",
    "        SN_Th_err.append(float(list_232_230_1s[i]))\n",
    "        SN_Sm_vals.append(float(list_152_149[i]))\n",
    "        SN_Sm_err.append(float(list_152_149_1s[i]))\n",
    "    elif aliquot_list[i].startswith('SBS'):\n",
    "        SB_Sm152_149_vals.append(float(list_152_149[i]))\n",
    "        SB_Sm152_149_err.append(float(list_152_149_1s[i]))\n",
    "        SB_Sm154_152_vals.append(float(list_154_152[i]))\n",
    "        SB_Sm154_152_err.append(float(list_154_152_1s[i]))\n",
    "        SB_Sm150_152_vals.append(float(list_150_152[i]))\n",
    "        SB_Sm150_152_err.append(float(list_150_152_1s[i]))\n",
    "        SB_Sm148_152_vals.append(float(list_148_152[i]))\n",
    "        SB_Sm148_152_err.append(float(list_148_152_1s[i]))\n",
    "        SB_Sm147_152_vals.append(float(list_147_152[i]))\n",
    "        SB_Sm147_152_err.append(float(list_147_152_1s[i]))\n",
    "        SB_Sm144_152_vals.append(float(list_144_152[i]))\n",
    "        SB_Sm144_152_err.append(float(list_144_152_1s[i]))\n",
    "    elif aliquot_list[i].startswith('SBU'):\n",
    "        SB_U_vals.append(float(list_238_236[i]))\n",
    "        SB_U_err.append(float(list_238_236_1s[i]))\n",
    "        SB_Th_vals.append(float(list_232_230[i]))\n",
    "        SB_Th_err.append(float(list_232_230_1s[i]))\n",
    "    else:\n",
    "        sample_U_vals.append(float(list_238_236[i]))\n",
    "        sample_U_err.append(float(list_238_236_1s[i]))\n",
    "        sample_Th_vals.append(float(list_232_230[i]))\n",
    "        sample_Th_err.append(float(list_232_230_1s[i]))\n",
    "        sample_Sm_vals.append(float(list_152_149[i]))\n",
    "        sample_Sm_err.append(float(list_152_149_1s[i]))\n",
    "        sample_names.append(aliquot_list[i])\n",
    "        \n",
    "fig, axs = plt.subplots(3)\n",
    "\n",
    "#this one plots 238/236\n",
    "axs[0].set_xlabel('SN name')\n",
    "axs[0].set_ylabel('$^{238}$U/$^{236}$U')\n",
    "\n",
    "#this one plots 232/230\n",
    "axs[1].set_xlabel('SN name')\n",
    "axs[1].set_ylabel('$^{232}$Th/$^{230}$Th')\n",
    "\n",
    "#this one plots 152/149\n",
    "axs[2].set_xlabel('SN name')\n",
    "axs[2].set_ylabel('$^{152}$Sm/$^{149}$Sm')\n",
    "\n",
    "axs[0].errorbar(SN_names, SN_U_vals, SN_U_err, fmt='o')\n",
    "axs[1].errorbar(SN_names, SN_Th_vals, SN_Th_err, fmt='o')\n",
    "axs[2].errorbar(SN_names, SN_Sm_vals, SN_Sm_err, fmt='o')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell below calculates the spike mass of U, Th, Sm, Zr, and Ca delivered to each of our samples using the spiked normal solutions as a reference. These values are assigned to the variables Sw_U, Sw_Th, Sw_Sm, Sw_Ca, and Sw_Zr that will be used in the next set of cells below this one. The fundamental equation for each element looks like this:\n",
    "\n",
    "### $N_w = S_w * \\frac{W_N}{W_S} * \\frac{Ab_S^A - R_m*Ab_S^B}{R_m*Ab_N^B - Ab_N^A}$\n",
    "\n",
    "which is rearranged to solve for $S_w$ such that:\n",
    "\n",
    "### $S_w = N_w * \\frac{W_S}{W_N} * \\frac{R_m*Ab_N^B - Ab_N^A}{Ab_S^A - R_m*Ab_S^B}$\n",
    "\n",
    "where: \n",
    "$N_w$ is the weight of the element in the normal solution,\n",
    "$S_w$ is the weight of the element in the spike solution,\n",
    "$W_N$ is the atomic weight of the element in the normal solution,\n",
    "$W_S$ is the atomic weight of the element in the spike solution,\n",
    "$Ab_S^A$ is the abundance of isotope A in the spike,\n",
    "$Ab_S^B$ is the abundance of isotope B in the spike,\n",
    "$Ab_N^A$ is the abundance of isotope A in the normal solution,\n",
    "$Ab_N^B$ is the abundance of isotope B in the spike solution, and\n",
    "$R_m$ is the measured ratio of the reference isotope over the enriched isotope. Each of these variables needs to have partial deriviatives calculated for error propagation, and this is also performed in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#start with spike calibration: calculate the mass of each spike using the spiked normals as a reference standard\n",
    "\n",
    "#function for weighted averages\n",
    "def weight_avg(ratio_list, d_ratio_list):\n",
    "    sum_W = 0\n",
    "    sum_XW = 0\n",
    "    \n",
    "    for i in range(0,len(ratio_list)):\n",
    "        W = 1/d_ratio_list[i]**2\n",
    "        XW = W * ratio_list[i]\n",
    "        sum_W = sum_W + W\n",
    "        sum_XW = sum_XW + XW\n",
    "    \n",
    "    w_avg = sum_XW/sum_W\n",
    "    w_avg_err = 1/math.sqrt(sum_W)\n",
    "    \n",
    "    return w_avg, w_avg_err\n",
    "\n",
    "#function for calculating S_w\n",
    "def spike_mass(Ab_A_nm, d_Ab_A_nm, Ab_B_nm, d_Ab_B_nm, Ab_A_spk, d_Ab_A_spk, Ab_B_spk, d_Ab_B_spk, Rm_spknm, d_Rm_spknm, Nw, d_Nw, Ws, d_Ws, Wn, d_Wn):\n",
    "    \n",
    "    #use some terms here for ease of doing error propagation: omega, gamma, kappa, zeta\n",
    "    omega = Nw*Ws*1/Wn\n",
    "    gamma = (Rm_spknm*Ab_B_nm - Ab_A_nm)/(Ab_A_spk - Rm_spknm*Ab_B_spk)\n",
    "    kappa = Ab_A_spk - Rm_spknm*Ab_B_spk\n",
    "    zeta = omega*Rm_spknm*Ab_B_nm - omega*Ab_A_nm\n",
    "    \n",
    "    #calculate all the partial derivatives for error propagation\n",
    "    d_Sw_d_Nw = Ws * gamma/Wn\n",
    "    d_Sw_d_Ws = Nw * gamma/Wn\n",
    "    d_Sw_d_Wn = -omega * gamma/Wn\n",
    "    d_Sw_d_Rm = (kappa*omega*Ab_B_nm + zeta*Ab_B_spk)/kappa**2\n",
    "    d_Sw_d_Ab_A_nm = -omega/kappa\n",
    "    d_Sw_d_Ab_B_nm = omega*Rm_spknm/kappa\n",
    "    d_Sw_d_Ab_A_spk = -zeta/kappa**2\n",
    "    d_Sw_d_Ab_B_spk = Rm_spknm*zeta/kappa\n",
    "    \n",
    "    #calculate Sw and d_Sw using terms above\n",
    "    Sw = omega * gamma\n",
    "    d_Sw = math.sqrt(d_Nw**2*d_Sw_d_Nw**2 + d_Ws**2*d_Sw_d_Ws**2 + d_Wn**2*d_Sw_d_Wn**2 + \n",
    "                     d_Rm_spknm**2*d_Sw_d_Rm**2 + d_Ab_B_nm**2*d_Sw_d_Ab_B_nm**2 + d_Ab_A_nm**2*d_Sw_d_Ab_A_nm**2 + \n",
    "                     d_Ab_B_spk**2*d_Sw_d_Ab_B_spk**2 + d_Ab_A_spk**2*d_Sw_d_Ab_A_spk**2)\n",
    "    \n",
    "    return Sw, d_Sw\n",
    "\n",
    "#calculate mass of spike in each spike normal, SW_U, SW_Th, SW_Sm\n",
    "\n",
    "#first up is U\n",
    "#natural/spike ratio in spike blank used to get abundances of isotopes in the spike\n",
    "Rs_U, d_Rs_U = weight_avg(SB_U_vals,SB_U_err)\n",
    "Ab_236U_spk = 1/(1+Rs_U/137.818+Rs_U)\n",
    "d_Ab_236U_spk = math.sqrt(d_Rs_U**2 * (-1*(1/137.818+1)/(1+Rs_U*(1/137.818+1))**2)**2)\n",
    "Ab_238U_spk = (1-Ab_236U_spk)/(1+1/137.818)\n",
    "d_Ab_238U_spk = math.sqrt(d_Ab_236U_spk**2 * (-1*1/(1+1/137.818))**2)\n",
    "Ab_235U_spk = Ab_238U_spk/137.818\n",
    "d_Ab_235U_spk = math.sqrt(d_Ab_238U_spk**2*(1/137.818)**2)\n",
    "\n",
    "#atomic weight of U in spike (g/mol)\n",
    "Ws_U = Ab_235U_spk*mass_235U + Ab_236U_spk*mass_236U + Ab_238U_spk*mass_238U\n",
    "d_Ws_U = math.sqrt(d_Ab_235U_spk**2 * mass_235U**2 + d_Ab_236U_spk**2 * mass_236U**2 + d_Ab_238U_spk**2 * mass_238U**2)\n",
    "\n",
    "#weight of U in the normal solution (ng)\n",
    "Nw_U = Vnm_UTh*concnm_U\n",
    "d_Nw_U = math.sqrt(d_Vnm_UTh**2*concnm_U**2 + d_concnm_U**2*Vnm_UTh**2)\n",
    "\n",
    "#get weighted average of spike normals\n",
    "Rm_spknm_U, d_Rm_spknm_U = weight_avg(SN_U_vals,SN_U_err)\n",
    "\n",
    "#mass of U spike in each spike normal (ng), used in subsequent calculations for mass of U in unknowns \n",
    "Sw_U, d_Sw_U = spike_mass(Ab_238U_nm, d_Ab_238U_nm, 0, 0, Ab_238U_spk, d_Ab_238U_spk, Ab_236U_spk, d_Ab_236U_spk, \n",
    "                          Rm_spknm_U, d_Rm_spknm_U, Nw_U, d_Nw_U, Ws_U, d_Ws_U, mass_U, 0)\n",
    "\n",
    "\n",
    "#next up is Th\n",
    "#natural/spike ratio in spike blank used to get abundances of isotopes in the spike\n",
    "Rs_Th, d_Rs_Th = weight_avg(SB_Th_vals,SB_Th_err)\n",
    "Ab_230Th_spk = 1/(1+Rs_Th)\n",
    "d_Ab_230Th_spk = math.sqrt(d_Rs_Th**2 * (-1/(1+Rs_Th)**2)**2)\n",
    "Ab_232Th_spk = 1 - Ab_230Th_spk\n",
    "d_Ab_232Th_spk = math.sqrt(d_Ab_230Th_spk**2 * 1**2)\n",
    "\n",
    "#atomic weight of Th in spike (g/mol)\n",
    "Ws_Th = Ab_230Th_spk*mass_230Th + Ab_232Th_spk*mass_232Th\n",
    "d_Ws_Th = math.sqrt(d_Ab_230Th_spk**2 * mass_230Th**2 + d_Ab_232Th_spk**2 * mass_232Th**2)\n",
    "\n",
    "#weight of Th in the normal solution (ng)\n",
    "Nw_Th = Vnm_UTh*concnm_Th\n",
    "d_Nw_Th = math.sqrt(d_Vnm_UTh**2*concnm_Th**2 + d_concnm_Th**2*Vnm_UTh**2)\n",
    "\n",
    "#get weighted average of spike normals\n",
    "Rm_spknm_Th, d_Rm_spknm_Th = weight_avg(SN_Th_vals,SN_Th_err)\n",
    "\n",
    "#mass of Th spike in each spike normal (ng), used in subsequent calculations for mass of U in unknowns \n",
    "Sw_Th, d_Sw_Th = spike_mass(1, 0, 0, 0, Ab_232Th_spk, d_Ab_232Th_spk, Ab_230Th_spk, d_Ab_230Th_spk, \n",
    "                          Rm_spknm_Th, d_Rm_spknm_Th, Nw_Th, d_Nw_Th, Ws_Th, d_Ws_Th, mass_Th, 0)\n",
    "\n",
    "\n",
    "#now we do Sm, but only if we're running apatites (not used for zircon)\n",
    "#initial if statment checks if we actually spiked our normal solutions with 149Sm\n",
    "#natural/spike ratio in spike blank used to get abundances of isotopes in the spike\n",
    "\n",
    "if SB_Sm152_149_vals:\n",
    "    Rs_Sm, d_Rs_Sm = weight_avg(SB_Sm152_149_vals,SB_Sm152_149_err)\n",
    "    Rs_154_152, d_Rs_154_152 = weight_avg(SB_Sm154_152_vals,SB_Sm154_152_err)\n",
    "    Rs_150_152, d_Rs_150_152 = weight_avg(SB_Sm150_152_vals,SB_Sm150_152_err)\n",
    "    Rs_148_152, d_Rs_148_152 = weight_avg(SB_Sm148_152_vals,SB_Sm148_152_err)\n",
    "    Rs_147_152, d_Rs_147_152 = weight_avg(SB_Sm147_152_vals,SB_Sm147_152_err)\n",
    "    Rs_144_152, d_Rs_144_152 = weight_avg(SB_Sm144_152_vals,SB_Sm144_152_err)\n",
    "    delta_Sm = 1 + Rs_154_152 + Rs_150_152 + Rs_148_152 + Rs_147_152 + Rs_144_152\n",
    "\n",
    "    Ab_149Sm_spk = 1/(1 + Rs_Sm*delta_Sm)\n",
    "    d_Ab_149Sm_spk = math.sqrt(d_Rs_Sm**2*(-delta_Sm/(1+Rs_Sm*delta_Sm)**2)**2 + \n",
    "                           d_Rs_144_152**2*(-Rs_Sm/(1+Rs_Sm*delta_Sm)**2)**2 + \n",
    "                           d_Rs_147_152**2*(-Rs_Sm/(1+Rs_Sm*delta_Sm)**2)**2 + \n",
    "                           d_Rs_148_152**2*(-Rs_Sm/(1+Rs_Sm*delta_Sm)**2)**2 + \n",
    "                           d_Rs_150_152**2*(-Rs_Sm/(1+Rs_Sm*delta_Sm)**2)**2 + \n",
    "                           d_Rs_154_152**2*(-Rs_Sm/(1+Rs_Sm*delta_Sm)**2)**2)\n",
    "    Ab_152Sm_spk = Rs_Sm * Ab_149Sm_spk\n",
    "    d_Ab_152Sm_spk = math.sqrt(d_Rs_Sm**2*Ab_149Sm_spk**2 + d_Ab_149Sm_spk**2*Rs_Sm**2)\n",
    "\n",
    "    Ab_144Sm_spk = Ab_152Sm_spk*Rs_144_152\n",
    "    d_Ab_144Sm_spk = math.sqrt(d_Rs_144_152**2*Ab_152Sm_spk**2 + d_Ab_152Sm_spk**2*Rs_144_152**2)\n",
    "    Ab_147Sm_spk = Ab_152Sm_spk*Rs_147_152\n",
    "    d_Ab_147Sm_spk = math.sqrt(d_Rs_147_152**2*Ab_152Sm_spk**2 + d_Ab_152Sm_spk**2*Rs_147_152**2)\n",
    "    Ab_148Sm_spk = Ab_152Sm_spk*Rs_148_152\n",
    "    d_Ab_148Sm_spk = math.sqrt(d_Rs_148_152**2*Ab_152Sm_spk**2 + d_Ab_152Sm_spk**2*Rs_148_152**2)\n",
    "    Ab_150Sm_spk = Ab_152Sm_spk*Rs_150_152\n",
    "    d_Ab_150Sm_spk = math.sqrt(d_Rs_150_152**2*Ab_152Sm_spk**2 + d_Ab_152Sm_spk**2*Rs_150_152**2)\n",
    "    Ab_154Sm_spk = Ab_152Sm_spk*Rs_154_152\n",
    "    d_Ab_154Sm_spk = math.sqrt(d_Rs_154_152**2*Ab_152Sm_spk**2 + d_Ab_152Sm_spk**2*Rs_154_152**2)\n",
    "\n",
    "    #atomic weight of Sm in spike (g/mol)\n",
    "    Ws_Sm = Ab_144Sm_spk*mass_144Sm + Ab_147Sm_spk*mass_147Sm + Ab_148Sm_spk*mass_148Sm + Ab_149Sm_spk*mass_149Sm + Ab_150Sm_spk*mass_150Sm + Ab_152Sm_spk*mass_152Sm + Ab_154Sm_spk*mass_154Sm\n",
    "    d_Ws_Sm = math.sqrt(d_Ab_144Sm_spk**2*mass_144Sm**2 + d_Ab_147Sm_spk**2*mass_147Sm**2 + d_Ab_148Sm_spk**2*mass_148Sm**2 \n",
    "                    + d_Ab_149Sm_spk**2*mass_149Sm**2 + d_Ab_150Sm_spk**2*mass_150Sm**2 + d_Ab_152Sm_spk**2*mass_152Sm**2 \n",
    "                    + d_Ab_154Sm_spk**2*mass_154Sm**2)\n",
    "\n",
    "    #weight of Sm in the normal solution (ng)\n",
    "    Nw_Sm = Vnm_Sm*concnm_Sm\n",
    "    d_Nw_Sm = math.sqrt(d_Vnm_Sm**2*concnm_Sm**2 + d_concnm_Sm**2*Vnm_Sm**2)\n",
    "\n",
    "    #get weighted average of spike normals\n",
    "    Rm_spknm_Sm, d_Rm_spknm_Sm = weight_avg(SN_Sm_vals,SN_Sm_err)\n",
    "\n",
    "    #mass of Sm spike in each spike normal (ng), used in subsequent calculations for mass of U in unknowns \n",
    "    Sw_Sm, d_Sw_Sm = spike_mass(Ab_152Sm_nm, d_Ab_152Sm_nm, Ab_149Sm_nm, d_Ab_149Sm_nm, Ab_152Sm_spk, d_Ab_152Sm_spk, Ab_149Sm_spk, d_Ab_149Sm_spk, \n",
    "                          Rm_spknm_Sm, d_Rm_spknm_Sm, Nw_Sm, d_Nw_Sm, Ws_Sm, d_Ws_Sm, mass_Sm, 0)\n",
    "else:\n",
    "    Sw_Sm = 0\n",
    "    d_Sw_Sm = 0\n",
    "\n",
    "print('Masses of each spike in an aliquot are: {0:.5f} +/- {1:.5f} ng for U, {2:.5f} +/- {3:.5f} ng for Th, and {4:.5f} +/- {5:.5f} ng for Sm.'.format(Sw_U,d_Sw_U,Sw_Th,d_Sw_Th,Sw_Sm,d_Sw_Sm))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our spike masses calibrated, we now use those spike masses to determine the mass of U, Th, and Sm in each unknown aliquot. The same fundamental equation is used:\n",
    "\n",
    "### $Nw_U = Sw_U * \\frac{Wn_U}{Ws_U} * \\frac{Ab_S^A - R_m*Ab_S^B}{R_m*Ab_U^B - Ab_U^A}$\n",
    "\n",
    "$Nw_U$ is the weight of the element (in this example U) in the unknown solution,\n",
    "$Sw_U$ is the weight of the element in the spike solution,\n",
    "$Wn_U$ is the atomic weight of the element in the unknown solution,\n",
    "$Ws_U$ is the atomic weight of the element in the spike solution,\n",
    "$Ab_S^A$ is the abundance of isotope A in the spike,\n",
    "$Ab_S^B$ is the abundance of isotope B in the spike,\n",
    "$Ab_U^A$ is the abundance of isotope A in the unknown solution,\n",
    "$Ab_U^B$ is the abundance of isotope B in the unknown solution, and\n",
    "$R_m$ is the measured ratio of the reference isotope over the enriched isotope. Each of these variables needs to have partial deriviatives calculated for error propagation, and this is also performed in the cell below. To simplify this math, once again we define:\n",
    "\n",
    "### $\\omega_U = Sw_U * \\frac{Wn_U}{Ws_U}$\n",
    "and\n",
    "### $\\gamma_U = \\frac{Ab_S^A - R_m*Ab_S^B}{R_m*Ab_U^B - Ab_U^A}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#function for calculating Nw_element\n",
    "def aliquot_mass(Ab_A_el, d_Ab_A_el, Ab_B_el, d_Ab_B_el, Ab_A_spk, d_Ab_A_spk, Ab_B_spk, d_Ab_B_spk, Rm_el, d_Rm_el, Sw, d_Sw, Ws, d_Ws, Wn, d_Wn):\n",
    "    \n",
    "    #use some terms here for ease of doing error propagation: omega, gamma, kappa, zeta\n",
    "    omega = Sw*Wn/Ws\n",
    "    gamma = (Ab_A_spk - Rm_el*Ab_B_spk)/(Rm_el*Ab_B_el - Ab_A_el)\n",
    "    kappa = Rm_el*Ab_B_el - Ab_A_el\n",
    "    zeta = omega*Ab_A_spk - omega*Rm_el*Ab_B_spk\n",
    "    \n",
    "    #calculate all the partial derivatives for error propagation\n",
    "    d_Nw_d_Sw = Wn * gamma/Ws\n",
    "    d_Nw_d_Ws = -omega*gamma/Ws\n",
    "    d_Nw_d_Wn = Sw * gamma/Ws\n",
    "    d_Nw_d_Rm = (-omega*Ab_B_spk*kappa - zeta*Ab_B_el)/kappa**2\n",
    "    d_Nw_d_Ab_A_el = -Rm_el*zeta/kappa**2\n",
    "    d_Nw_d_Ab_B_el = zeta/kappa**2\n",
    "    d_Nw_d_Ab_A_spk = omega/kappa\n",
    "    d_Nw_d_Ab_B_spk = -Rm_el*omega/kappa\n",
    "    \n",
    "    #calculate Sw and d_Sw using terms above\n",
    "    Nw = omega * gamma\n",
    "    d_Nw = math.sqrt(d_Sw**2*d_Nw_d_Sw**2 + d_Ws**2*d_Nw_d_Ws**2 + d_Wn**2*d_Nw_d_Wn**2 + \n",
    "                     d_Rm_el**2*d_Nw_d_Rm**2 + d_Ab_B_el**2*d_Nw_d_Ab_B_el**2 + d_Ab_A_el**2*d_Nw_d_Ab_A_el**2 + \n",
    "                     d_Ab_B_spk**2*d_Nw_d_Ab_B_spk**2 + d_Ab_A_spk**2*d_Nw_d_Ab_A_spk**2)\n",
    "    \n",
    "    return Nw, d_Nw\n",
    "\n",
    "#calculate U, Th, and Sm mass (ng) for each sample, add to sample mass lists\n",
    "sample_U_mass = []\n",
    "sample_U_mass_1s = []\n",
    "sample_238U_mol =[]\n",
    "sample_238U_mol_1s = []\n",
    "sample_235U_mol = []\n",
    "sample_235U_mol_1s = []\n",
    "sample_Th_mass = []\n",
    "sample_Th_mass_1s = []\n",
    "sample_Th_mol = []\n",
    "sample_Th_mol_1s = []\n",
    "sample_Sm_mass = []\n",
    "sample_Sm_mass_1s = []\n",
    "sample_Sm_mol = []\n",
    "sample_Sm_mol_1s = []\n",
    "\n",
    "\n",
    "for i in range(len(sample_names)):\n",
    "    U_mass, U_mass_1s = aliquot_mass(Ab_238U_nm, d_Ab_238U_nm, 0, 0, Ab_238U_spk, d_Ab_238U_spk, Ab_236U_spk, d_Ab_236U_spk, \n",
    "                                     sample_U_vals[i], sample_U_err[i], Sw_U, d_Sw_U, Ws_U, d_Ws_U, mass_U, 0)\n",
    "    sample_U_mass.append(U_mass)\n",
    "    sample_U_mass_1s.append(U_mass_1s)\n",
    "    sample_238U_mol.append(U_mass * 1e-9 * 1/mass_U*Ab_238U_nm)\n",
    "    sample_238U_mol_1s.append(U_mass_1s * 1e-9 * 1/mass_U*Ab_238U_nm)\n",
    "    sample_235U_mol.append(U_mass * 1e-9 * 1/mass_U*Ab_235U_nm)\n",
    "    sample_235U_mol_1s.append(U_mass_1s * 1e-9 * 1/mass_U*Ab_235U_nm)\n",
    "    \n",
    "    Th_mass, Th_mass_1s = aliquot_mass(1, 0, 0, 0, Ab_232Th_spk, d_Ab_232Th_spk, Ab_230Th_spk, d_Ab_230Th_spk, \n",
    "                                     sample_Th_vals[i], sample_Th_err[i], Sw_Th, d_Sw_Th, Ws_Th, d_Ws_Th, mass_Th, 0)\n",
    "    sample_Th_mass.append(Th_mass)\n",
    "    sample_Th_mass_1s.append(Th_mass_1s)\n",
    "    sample_Th_mol.append(Th_mass * 1e-9 * 1/(mass_Th))\n",
    "    sample_Th_mol_1s.append(Th_mass_1s * 1e-9 * 1/(mass_Th))\n",
    "    \n",
    "    if SB_Sm152_149_vals:\n",
    "        Sm_mass, Sm_mass_1s = aliquot_mass(Ab_152Sm_nm, d_Ab_152Sm_nm, Ab_149Sm_nm, d_Ab_149Sm_nm, Ab_152Sm_spk, d_Ab_152Sm_spk, Ab_149Sm_spk, d_Ab_149Sm_spk, \n",
    "                                     sample_Sm_vals[i], sample_Sm_err[i], Sw_Sm, d_Sw_Sm, Ws_Sm, d_Ws_Sm, mass_Sm, 0)\n",
    "        sample_Sm_mass.append(Sm_mass)\n",
    "        sample_Sm_mass_1s.append(Sm_mass_1s)\n",
    "        sample_Sm_mol.append(Sm_mass * 1e-9 * 1/mass_Sm*Ab_147Sm_nm)\n",
    "        sample_Sm_mol_1s.append(Th_mass_1s * 1e-9 * 1/mass_Sm*Ab_147Sm_nm)\n",
    "    else:\n",
    "        sample_Sm_mass.append(0)\n",
    "        sample_Sm_mass_1s.append(0)\n",
    "        sample_Sm_mol.append(0)\n",
    "        sample_Sm_mol_1s.append(0)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell below brings in the line summary json and calculates the uncorrected date for each sample. The age equation, which can't be solved for t using algebra, has the form:\n",
    "\n",
    "### $^{4}He_{mol} = 8 * ^{238}U_{mol}(\\exp(\\lambda_{238} * t) - 1) + 7 * ^{235}U_{mol}(\\exp(\\lambda_{235} * t) - 1) + 6 * ^{232}Th_{mol}(\\exp(\\lambda_{232} * t) - 1) + ^{147}Sm_{mol}(\\exp(\\lambda_{147} * t) - 1)$\n",
    "\n",
    "We use a Newton-Raphson method for determine date, where:\n",
    "\n",
    "### $t_1 = t_0 - \\frac{f(t_0)}{f'(t_0)}$\n",
    "\n",
    "### $f(t_0) = 8 * ^{238}U_{mol}(\\exp(\\lambda_{238} * t) - 1) + 7 * ^{235}U_{mol}(\\exp(\\lambda_{235} * t) - 1) + 6 * ^{232}Th_{mol}(\\exp(\\lambda_{232} * t) - 1) + ^{147}Sm_{mol}(\\exp(\\lambda_{147} * t) - 1) - ^{4}He_{mol}$\n",
    "\n",
    "and \n",
    "\n",
    "### $f'(t_0) = 8 * ^{238}U_{mol} * \\lambda_{238} * \\exp(\\lambda_{238} * t) + 7 * ^{235}U_{mol} * \\lambda_{235} * \\exp(\\lambda_{238} * t) + 6 * ^{232}Th_{mol} * \\lambda_{232} * \\exp(\\lambda_{232} * t) + ^{147}Sm_{mol} * \\lambda_{147} * \\exp(\\lambda_{147} * t)$\n",
    "\n",
    "The error is solved by taking the partial derivatives of terms within an approximation of the age equation using the second term of the Taylor Series expansion. This is a pretty dang good approximation as it gets the error value within less than 1% of the \"true\" value. That equation looks like:\n",
    "\n",
    "## $ t = \\frac{^{4}He_{mol}}{8 * ^{238}U_{mol} * \\lambda_{238} + 7 * ^{235}U_{mol} * \\lambda_{235} + 6 * ^{232}Th_{mol} * \\lambda_{232} + ^{147}Sm_{mol} * \\lambda_{147}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#line summary json file\n",
    "He_file = project + '_He_data.json'\n",
    "\n",
    "with open(He_file, 'r') as j_file:\n",
    "     He_data = json.loads(j_file.read())\n",
    "\n",
    "#lists of He data \n",
    "sample_He_vol = []\n",
    "sample_He_vol_1s = []\n",
    "sample_He_mol = []\n",
    "sample_He_mol_1s = []\n",
    "sample_He_notes = []\n",
    "sample_type = []\n",
    "\n",
    "#sample uncorrected date and error list\n",
    "sample_uncorrected_date = []\n",
    "sample_uncorrected_date_1s= []\n",
    "\n",
    "for index in sample_names:\n",
    "    sample_He_vol.append(He_data[index][0])\n",
    "    sample_He_vol_1s.append(He_data[index][1])\n",
    "    sample_He_mol.append(He_data[index][2])\n",
    "    sample_He_mol_1s.append(He_data[index][3])\n",
    "    sample_He_notes.append(He_data[index][4])\n",
    "    sample_type.append(He_data[index][5])\n",
    "\n",
    "for index in range(len(sample_names)):\n",
    "    \n",
    "    #this is here because there could be cases where apatites and zircons are run in the same ICP-MS routine\n",
    "    if sample_type[index] == 'zirc':\n",
    "        sample_Sm_mol[index] = 0\n",
    "        sample_Sm_mol_1s[index] = 0\n",
    "        sample_Sm_mol.append(0)\n",
    "        sample_Sm_mol_1s.append(0)\n",
    "    \n",
    "    #calculate uncorrected dates using Newton-Raphson method, get the guess to the 1 year precision\n",
    "    tolerance = 1\n",
    "    date_guess = 5e7\n",
    "    uncorrected_date = 0\n",
    "    f_date = 8*sample_238U_mol[index]*(math.exp(lambda_238*date_guess) - 1) + 7*sample_235U_mol[index]*(math.exp(lambda_235*date_guess) - 1) + 6*sample_Th_mol[index]*(math.exp(lambda_232*date_guess) - 1) + sample_Sm_mol[index]*(math.exp(lambda_147*date_guess) - 1) - sample_He_mol[index]\n",
    "    f_date_prime = 8*sample_238U_mol[index]*lambda_238*math.exp(lambda_238*date_guess)  + 7*sample_235U_mol[index]*lambda_235*math.exp(lambda_235*date_guess) + 6*sample_Th_mol[index]*lambda_232*math.exp(lambda_232*date_guess) + sample_Sm_mol[index]*lambda_147*math.exp(lambda_147*date_guess)\n",
    "    date_diff = date_guess\n",
    "    \n",
    "    while abs(date_diff) > tolerance:\n",
    "        uncorrected_date = date_guess - f_date/f_date_prime\n",
    "        date_diff = date_guess - uncorrected_date\n",
    "        date_guess = uncorrected_date\n",
    "        f_date = 8*sample_238U_mol[index]*(math.exp(lambda_238*date_guess) - 1) + 7*sample_235U_mol[index]*(math.exp(lambda_235*date_guess) - 1) + 6*sample_Th_mol[index]*(math.exp(lambda_232*date_guess) - 1) + sample_Sm_mol[index]*(math.exp(lambda_147*date_guess) - 1) - sample_He_mol[index]\n",
    "        f_date_prime = 8*sample_238U_mol[index]*lambda_238*math.exp(lambda_238*date_guess)  + 7*sample_235U_mol[index]*lambda_235*math.exp(lambda_235*date_guess) + 6*sample_Th_mol[index]*lambda_232*math.exp(lambda_232*date_guess) + sample_Sm_mol[index]*lambda_147*math.exp(lambda_147*date_guess)\n",
    "        \n",
    "    sample_uncorrected_date.append(uncorrected_date/1e6)\n",
    "    \n",
    "    #calculate uncorrected date error using the 2nd term in Taylor Series of age equation\n",
    "    alpha = 8*sample_238U_mol[index]*lambda_238 + 7*sample_235U_mol[index]*lambda_235 + 6*sample_Th_mol[index]*lambda_232 + sample_Sm_mol[index]*lambda_147\n",
    "    \n",
    "    d_t_d_4He = 1/alpha\n",
    "    d_t_d_238U = (-sample_He_mol[index]*8*lambda_238)/alpha**2\n",
    "    d_t_d_235U = (-sample_He_mol[index]*7*lambda_235)/alpha**2\n",
    "    d_t_d_232Th = (-sample_He_mol[index]*6*lambda_232)/alpha**2\n",
    "    d_t_d_147Sm = (-sample_He_mol[index]*lambda_147)/alpha**2\n",
    "    \n",
    "    uncorrected_date_1s = math.sqrt(d_t_d_4He**2 * sample_He_mol_1s[index]**2 + \n",
    "                                    d_t_d_238U**2 * sample_238U_mol_1s[index]**2 + \n",
    "                                    d_t_d_235U**2 * sample_235U_mol_1s[index]**2 + \n",
    "                                    d_t_d_232Th**2 * sample_Th_mol_1s[index]**2 + \n",
    "                                    d_t_d_147Sm**2 * sample_Sm_mol_1s[index]**2)\n",
    "    \n",
    "    sample_uncorrected_date_1s.append(uncorrected_date_1s/1e6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we bring in Ft data for each sample. Two Ft methods are used here. We follow the equations from Farley, 2002:\n",
    "\n",
    "### $F_{t} = 1 + a_{1}\\beta + a_{2}\\beta^{2}$\n",
    "\n",
    "### $^{Mean}F_{t} = a_{238} * ^{^{238}U}F_{t} + (1 - a_{238}) * ^{^{232}Th}F_{t}$\n",
    "\n",
    "### $a_{238} = (1.04 + 0.245(Th/U))^{-1}$\n",
    "\n",
    "where $a_{1}$ and $a_{2}$ are U and Th specific constants listed in Farley (2002) Table 1 and set in the first cell of this notebook. $\\beta$ is the surface area to volume ratio, which for a hexagonal prism (apatite) is given by:\n",
    "\n",
    "### $\\beta = \\frac{2.31L + 2R}{RL}$\n",
    "\n",
    "where $L$ is the length and $R$ is the half-width of the grain. For zircon grains, we follow the equations described in Reiners et al. (2005, AJS) and Hourigan et al. (2005), with $a_{1}$ and $a_{2}$ values listed in Farley (2002) Table 1. Two equations are used, depending upon the degree of abrasion: tetragonal prism with pyramidal terminations ($z$, abrasion index 1-3), or prolate spheroid ($ps$, abrasion index 4-5). The basic $F_{t}$ equation has the same form as the Farley (2002) equation used above. However, the equation for $\\beta$ is different. For the tetragonal prism $\\beta$, the equations are (where $SA$ is surface area and $V$ is volume):\n",
    "\n",
    "### $SA_{z} = 4(l - h_{1} - h_{2})(r_{1} + r_{2}) + 2r_{1}a + 2r_{2}b$\n",
    "\n",
    "### $V_{z} = 4r_{1}r_{2}[(l - h_{1} - h_{2}) + \\frac{1}{3}(h_{1} + h_{2})]$\n",
    "\n",
    "where $l$ is grain length, $r_1$ and $r_2$ are prism half-widths, $h_1$ and $h_2$ are pyramidal termination heights, and $a$ and $b$ are related to the slant heights of the terminations as determined by:\n",
    "\n",
    "### $a = \\sqrt{h_{1}^{2} + r_{2}^{2}} + \\sqrt{h_{2}^{2} + r_{2}^{2}}$\n",
    "\n",
    "### $b = \\sqrt{h_{1}^{2} + r_{1}^{2}} + \\sqrt{h_{2}^{2} + r_{1}^{2}}$\n",
    "\n",
    "For the prolate spheroid case, the $\\beta$ equations are (NOTE: Reiners et al. (2005, AJS) has a number of typos in the reported equations. I have fixed these below by adding a ^2 to the appropriate l/2 terms):\n",
    "\n",
    "### $SA_{ps} = 2\\pi{}r^{2} + \\frac{2\\pi{}r(l/2)^{2}}{\\sqrt{(l/2)^{2} - r^{2}}}\\sin^{-1}[\\frac{\\sqrt{(l/2)^{2} - r^{2}}}{(l/2)}]$\n",
    "\n",
    "### $V_{ps} = \\frac{2}{3}\\pi{}r^{2}l$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#lists for Ft data\n",
    "Ft238U_list = []\n",
    "Ft235U_list = []\n",
    "FtTh_list = []\n",
    "FtSm_list = []\n",
    "Ft_list = []\n",
    "morph_note_list = []\n",
    "Equiv_rad_list = []\n",
    "grain_mass_morph_list = []\n",
    "\n",
    "#csv file for Ft data\n",
    "Ft_file = project + '_Ft_data.csv'\n",
    "\n",
    "#load the Ft data into a dictionary to look up by sample name\n",
    "with open(Ft_file, 'r', newline = '') as Ft_in_file:\n",
    "    csv_reader_Ft = csv.reader(Ft_in_file)\n",
    "    \n",
    "    #throw away the first 8 lines\n",
    "    for i in range(8):\n",
    "        next(csv_reader_Ft)\n",
    "    \n",
    "    #fill in Ft dictionary, indexed by sample name\n",
    "    Ft_dict = {}\n",
    "    for line in csv_reader_Ft:\n",
    "        #sometimes there might be blank rows separating Ft data\n",
    "        if line:\n",
    "            Ft_dict[line[0]] = line[1:]\n",
    "\n",
    "#for loop through the sample names and calculate Ft for ap or zirc (z or ps)\n",
    "\n",
    "for i in range(len(sample_names)):    \n",
    "    Ft_sample_name =  sample_names[i]\n",
    "    morph_note_list.append(Ft_dict[Ft_sample_name][0])\n",
    "    \n",
    "    A_238 = 1/(1.04 + 0.245*(sample_Th_mol[i]/(sample_238U_mol[i]+sample_235U_mol[i])))\n",
    "    \n",
    "    if sample_type[i] == 'ap':       \n",
    "        if Ft_dict[Ft_sample_name][1] == 'DUR' or Ft_dict[Ft_sample_name][1] == 'dur':\n",
    "            Ft238U_list.append(1)\n",
    "            Ft235U_list.append(1)\n",
    "            FtTh_list.append(1)\n",
    "            FtSm_list.append(1)\n",
    "            Ft_list.append(1)\n",
    "            Equiv_rad_list.append('NaN')\n",
    "            grain_mass_morph_list.append('NaN')\n",
    "        \n",
    "        else:\n",
    "            broken = float(Ft_dict[Ft_sample_name][1])   \n",
    "            \n",
    "            grain_length_avg = (float(Ft_dict[Ft_sample_name][2]) + float(Ft_dict[Ft_sample_name][3]))/2\n",
    "            grain_radius = ((float(Ft_dict[Ft_sample_name][4]) + float(Ft_dict[Ft_sample_name][5]))/2)/2\n",
    "            beta_ap = (2.31*grain_length_avg*broken + 2*grain_radius)/(broken*grain_length_avg*grain_radius)\n",
    "            \n",
    "            Ft238U_list.append(1 + a1_238U_ap*beta_ap + a2_238U_ap*beta_ap**2)\n",
    "            Ft235U_list.append(1 + a1_235U_ap*beta_ap + a2_235U_ap*beta_ap**2)\n",
    "            FtTh_list.append(1 + a1_Th_ap*beta_ap + a2_Th_ap*beta_ap**2)\n",
    "            FtSm_list.append(-0.09158*Ft238U_list[i]**2 + 0.46974*Ft238U_list[i] + 0.61895)\n",
    "            Ft_list.append(A_238*Ft238U_list[i] + (1-A_238)*FtTh_list[i])\n",
    "            \n",
    "            Equiv_rad_list.append(3/beta_ap)\n",
    "            grain_volume = 2.598*grain_radius**2*grain_length_avg*1e-12 #in cc   \n",
    "            grain_mass_morph_list.append(grain_volume * ap_density) #in g\n",
    "            \n",
    "    elif sample_type[i] == 'zirc':\n",
    "            abrasion = float(Ft_dict[Ft_sample_name][1])\n",
    "            grain_length_avg = (float(Ft_dict[Ft_sample_name][2]) + float(Ft_dict[Ft_sample_name][3]))/2\n",
    "            grain_radius_1 = float(Ft_dict[Ft_sample_name][4])/2\n",
    "            grain_radius_2 = float(Ft_dict[Ft_sample_name][5])/2\n",
    "            grain_width_avg = (float(Ft_dict[Ft_sample_name][4]) + float(Ft_dict[Ft_sample_name][5]))/2\n",
    "            \n",
    "            # there could be cases for very rounded grains where what's a length vs. a width is not apparent\n",
    "            # this if statement ensures (so that the SA_ps math works) that the length is ALWAYS the longest dimension\n",
    "            if grain_width_avg > grain_length_avg:\n",
    "                grain_length_avg_old = grain_length_avg\n",
    "                grain_length_avg = grain_width_avg\n",
    "                grain_width_avg = grain_length_avg_old\n",
    "            \n",
    "            grain_radius_avg = grain_width_avg/2\n",
    "            \n",
    "            if abrasion <= 3:\n",
    "                a_z = math.sqrt(float(Ft_dict[Ft_sample_name][6])**2 + grain_radius_2**2) + math.sqrt(float(Ft_dict[Ft_sample_name][7])**2 + grain_radius_2**2)\n",
    "                b_z = math.sqrt(float(Ft_dict[Ft_sample_name][6])**2 + grain_radius_1**2) + math.sqrt(float(Ft_dict[Ft_sample_name][7])**2 + grain_radius_1**2)\n",
    "                SA_z = 4*(grain_length_avg - float(Ft_dict[Ft_sample_name][6]) - float(Ft_dict[Ft_sample_name][7]))*(grain_radius_1 + grain_radius_2) + 2*grain_radius_1*a_z + 2*grain_radius_2*b_z\n",
    "                grain_volume = 4*grain_radius_1*grain_radius_2*((grain_length_avg - float(Ft_dict[Ft_sample_name][6]) - float(Ft_dict[Ft_sample_name][7])) + (1/3)*(float(Ft_dict[Ft_sample_name][6]) + float(Ft_dict[Ft_sample_name][7])))  \n",
    "                beta_z = SA_z/grain_volume\n",
    "            \n",
    "                Ft238U_list.append(1 + a1_238U_zirc*beta_z + a2_238U_zirc*beta_z**2)\n",
    "                Ft235U_list.append(1 + a1_235U_zirc*beta_z + a2_235U_zirc*beta_z**2)\n",
    "                FtTh_list.append(1 + a1_Th_zirc*beta_z + a2_Th_zirc*beta_z**2)\n",
    "                FtSm_list.append(1)\n",
    "                Ft_list.append(A_238*Ft238U_list[i] + (1-A_238)*FtTh_list[i])\n",
    "            \n",
    "                Equiv_rad_list.append(3/beta_z)   \n",
    "                grain_mass_morph_list.append(grain_volume*1e-12 * zirc_density) #in g\n",
    "            else:\n",
    "                SA_ps = 2*math.pi*grain_radius_avg**2 + ((2*math.pi*grain_radius_avg*(grain_length_avg/2)**2)/(math.sqrt((grain_length_avg/2)**2 - grain_radius_avg**2)))*math.asin(math.sqrt((grain_length_avg/2)**2 - grain_radius_avg**2)/(grain_length_avg/2))\n",
    "                grain_volume = (2/3)*math.pi*grain_radius_avg**2*grain_length_avg\n",
    "                beta_ps = SA_ps/grain_volume\n",
    "            \n",
    "                Ft238U_list.append(1 + a1_238U_zirc*beta_ps + a2_238U_zirc*beta_ps**2)\n",
    "                Ft235U_list.append(1 + a1_235U_zirc*beta_ps + a2_235U_zirc*beta_ps**2)\n",
    "                FtTh_list.append(1 + a1_Th_zirc*beta_ps + a2_Th_zirc*beta_ps**2)\n",
    "                FtSm_list.append(1)\n",
    "                Ft_list.append(A_238*Ft238U_list[i] + (1-A_238)*FtTh_list[i])\n",
    "            \n",
    "                Equiv_rad_list.append(3/beta_ps)   \n",
    "                grain_mass_morph_list.append(grain_volume*1e-12 * zirc_density) #in g\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we use the Ft corrections on each of the parent isotopes and then solve for the corrected dates using the Newton-Raphson method as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sample corrected date and error list\n",
    "sample_corrected_date = []\n",
    "sample_corrected_date_1s= []\n",
    "\n",
    "for index in range(len(sample_names)):\n",
    "    \n",
    "    #calculate corrected dates using Newton-Raphson method, get the guess to the 1 year precision\n",
    "    tolerance = 1\n",
    "    date_guess = 5e7\n",
    "    uncorrected_date = 0\n",
    "    Ft_238U = sample_238U_mol[index] * Ft238U_list[index]\n",
    "    Ft_238U_1s = sample_238U_mol_1s[index] * Ft238U_list[index]\n",
    "    Ft_235U = sample_235U_mol[index] * Ft235U_list[index]\n",
    "    Ft_235U_1s = sample_235U_mol_1s[index] * Ft235U_list[index]\n",
    "    Ft_Th = sample_Th_mol[index] * FtTh_list[index]\n",
    "    Ft_Th_1s = sample_Th_mol_1s[index] * FtTh_list[index]\n",
    "    Ft_Sm = sample_Sm_mol[index] * FtSm_list[index]\n",
    "    Ft_Sm_1s = sample_Sm_mol_1s[index] * FtSm_list[index]\n",
    "    \n",
    "    f_date = 8*Ft_238U*(math.exp(lambda_238*date_guess) - 1) + 7*Ft_235U*(math.exp(lambda_235*date_guess) - 1) + 6*Ft_Th*(math.exp(lambda_232*date_guess) - 1) + Ft_Sm*(math.exp(lambda_147*date_guess) - 1) - sample_He_mol[index]\n",
    "    f_date_prime = 8*Ft_238U*lambda_238*math.exp(lambda_238*date_guess)  + 7*Ft_235U*lambda_235*math.exp(lambda_235*date_guess) + 6*Ft_Th*lambda_232*math.exp(lambda_232*date_guess) + Ft_Sm*lambda_147*math.exp(lambda_147*date_guess)\n",
    "    date_diff = date_guess\n",
    "    \n",
    "    while abs(date_diff) > tolerance:\n",
    "        corrected_date = date_guess - f_date/f_date_prime\n",
    "        date_diff = date_guess - corrected_date\n",
    "        date_guess = corrected_date\n",
    "        f_date = 8*Ft_238U*(math.exp(lambda_238*date_guess) - 1) + 7*Ft_235U*(math.exp(lambda_235*date_guess) - 1) + 6*Ft_Th*(math.exp(lambda_232*date_guess) - 1) + Ft_Sm*(math.exp(lambda_147*date_guess) - 1) - sample_He_mol[index]\n",
    "        f_date_prime = 8*Ft_238U*lambda_238*math.exp(lambda_238*date_guess)  + 7*Ft_235U*lambda_235*math.exp(lambda_235*date_guess) + 6*Ft_Th*lambda_232*math.exp(lambda_232*date_guess) + Ft_Sm*lambda_147*math.exp(lambda_147*date_guess)\n",
    "    \n",
    "    sample_corrected_date.append(corrected_date/1e6)\n",
    "    \n",
    "    #calculate corrected date error using the 2nd term in Taylor Series of age equation\n",
    "    alpha = 8*Ft_238U*lambda_238 + 7*Ft_235U*lambda_235 + 6*Ft_Th*lambda_232 + Ft_Sm*lambda_147\n",
    "    \n",
    "    d_t_d_4He = 1/alpha\n",
    "    d_t_d_238U = (-sample_He_mol[index]*8*lambda_238)/alpha**2\n",
    "    d_t_d_235U = (-sample_He_mol[index]*7*lambda_235)/alpha**2\n",
    "    d_t_d_232Th = (-sample_He_mol[index]*6*lambda_232)/alpha**2\n",
    "    d_t_d_147Sm = (-sample_He_mol[index]*lambda_147)/alpha**2\n",
    "    \n",
    "    corrected_date_1s = math.sqrt(d_t_d_4He**2 * sample_He_mol_1s[index]**2 + \n",
    "                                    d_t_d_238U**2 * Ft_238U_1s**2 + \n",
    "                                    d_t_d_235U**2 * Ft_235U_1s**2 + \n",
    "                                    d_t_d_232Th**2 * Ft_Th_1s**2 + \n",
    "                                    d_t_d_147Sm**2 * Ft_Sm_1s**2)\n",
    "    \n",
    "    sample_corrected_date_1s.append(corrected_date_1s/1e6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final step is to combine the various lists into a report that is exported to a csv file. The data lists are combined into a data frame to accomplish this task."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create the sample data frame\n",
    "sample_frame = pd.DataFrame(columns = ['analysis_notes','morph_notes','corrected_date','corrected_1s','corrected_1s_%','radius',\n",
    "                                       'Ft','uncorrected_date','uncorrected_1s','uncorrected_1s_%','U_ng','U_ng_1s','U_ng_1s_%',\n",
    "                                       'Th_ng','Th_ng_1s','Th_ng_1s_%','Sm_ng','Sm_ng_1s','Sm_ng_1s_%','He_pmol','He_pmol_1s','He_pmol_1s_%',\n",
    "                                       'mass_morph_g','U_ppm','Th_ppm','Sm_ppm','eU'], index = [])\n",
    "\n",
    "#fill in the frame\n",
    "for i in range(len(sample_names)):\n",
    "    \n",
    "    He_pmol = sample_He_mol[i] * 1e12\n",
    "    He_pmol_1s = sample_He_mol_1s[i] * 1e12\n",
    "    \n",
    "    if Ft_dict[sample_names[i]][1] == 'DUR' or Ft_dict[sample_names[i]][1] == 'dur':\n",
    "        U_ppm = 0\n",
    "        Th_ppm = 0\n",
    "        Sm_ppm = 0\n",
    "        eU = 0\n",
    "    else:\n",
    "        U_ppm = (sample_U_mass[i]*1e-3)/grain_mass_morph_list[i]\n",
    "        Th_ppm = (sample_Th_mass[i]*1e-3)/grain_mass_morph_list[i]\n",
    "        Sm_ppm = (sample_Sm_mass[i]*1e-3)/grain_mass_morph_list[i]\n",
    "        eU = U_ppm + 0.235*Th_ppm\n",
    "    \n",
    "    #for zircons, need to make sure not to attempt a divide by zero for Sm % error\n",
    "    if Sm_ppm > 0:\n",
    "        Sm_ng_1s_rel = 100*(sample_Sm_mass_1s[i]/sample_Sm_mass[i])\n",
    "    else:\n",
    "        Sm_ng_1s_rel = 0\n",
    "\n",
    "    row_data = pd.DataFrame({'analysis_notes':sample_He_notes[i],'morph_notes':morph_note_list[i],\n",
    "                            'corrected_date':sample_corrected_date[i],'corrected_1s':sample_corrected_date_1s[i],\n",
    "                            'corrected_1s_%':100*(sample_corrected_date_1s[i]/sample_corrected_date[i]),\n",
    "                            'radius':Equiv_rad_list[i],'Ft':Ft_list[i],'uncorrected_date':sample_uncorrected_date[i],\n",
    "                            'uncorrected_1s':sample_uncorrected_date_1s[i],'uncorrected_1s_%':100*(sample_uncorrected_date_1s[i]/sample_uncorrected_date[i]),\n",
    "                            'U_ng':sample_U_mass[i],'U_ng_1s':sample_U_mass_1s[i], 'U_ng_1s_%':100*(sample_U_mass_1s[i]/sample_U_mass[i]),\n",
    "                            'Th_ng':sample_Th_mass[i],'Th_ng_1s':sample_Th_mass_1s[i], 'Th_ng_1s_%':100*(sample_Th_mass_1s[i]/sample_Th_mass[i]), \n",
    "                            'Sm_ng':sample_Sm_mass[i], 'Sm_ng_1s':sample_Sm_mass_1s[i], 'Sm_ng_1s_%':Sm_ng_1s_rel,\n",
    "                            'He_pmol':He_pmol,'He_pmol_1s':He_pmol_1s, 'He_pmol_1s_%':100*(He_pmol_1s/He_pmol),\n",
    "                            'mass_morph_g':grain_mass_morph_list[i],'U_ppm':U_ppm,'Th_ppm':Th_ppm,'Sm_ppm':Sm_ppm,'eU':eU}, \n",
    "                            columns = ['analysis_notes','morph_notes','corrected_date','corrected_1s','corrected_1s_%','radius',\n",
    "                                       'Ft','uncorrected_date','uncorrected_1s','uncorrected_1s_%','U_ng','U_ng_1s','U_ng_1s_%','Th_ng','Th_ng_1s',\n",
    "                                       'Th_ng_1s_%', 'Sm_ng','Sm_ng_1s','Sm_ng_1s_%','He_pmol','He_pmol_1s','He_pmol_1s_%','mass_morph_g','U_ppm','Th_ppm',\n",
    "                                       'Sm_ppm','eU'], index = [sample_names[i]])\n",
    "    sample_frame = pd.concat([sample_frame, row_data])\n",
    "    \n",
    "sample_frame.to_csv(project + '_' + str(bomb_run) + '_report.csv')\n",
    "\n",
    "display(sample_frame)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "vscode": {
   "interpreter": {
    "hash": "1f0c7602c82e39efa19a01e5e068584db7a6d17aff8711ab06660aac81377393"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}