Arrhenius_template.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Arrhenius plot and diffusion kinetic parameters for step-heating experiments\n",
    "\n",
    "This notebook uses the reduced JSON data file of each step-heating experiments from the Helium Analysis Lab (HAL) at the University of Illinois Urbana-Champaign to calculate diffusivity, create Arrhenius plots and calculate the diffusion kinetic parameters for the zircon (U-Th)/He system. This notebook was written by Ryan Sigat and is constantly under construction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import packages\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#defining some constants\n",
    "gas_constant = 8.314472   # universal gas constant in J/K*mol\n",
    "avogadro = 6.022e+23   # in 1/mol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Input sample specific values\n",
    "\n",
    "In the cell below, input all necessary information for the specific sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sample_ID = ''\n",
    "\n",
    "# slab dimensions in microns\n",
    "width = \n",
    "length_1 = \n",
    "length_2 = \n",
    "\n",
    "# insert zircon density\n",
    "zirc_rho = 4.36  # in g/cm^3\n",
    "\n",
    "# dose\n",
    "dose_UPb =  # time-integrated dose in alphas/gram\n",
    "dose_UTHe = \n",
    "dose_eq1 =   # equivalent dose (model 1) in alphas/gram\n",
    "dose_eq2 =  # equivalent dose (model 1) in alphas/gram\n",
    "\n",
    "half_width = (width/2) / 1e4 # in cm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculating volume \n",
    "\n",
    "vol = (width * length_1 * length_2) / 1e12  # in cm^3\n",
    "\n",
    "#calculating volumetric mol of He\n",
    "#make sure to choose the alpha dose value you want to use\n",
    "\n",
    "dose = \n",
    "total_fraction = 1  # fraction of gas released after total degassing (varies depending on which dose is used)\n",
    "He_mol_total = (total_fraction * (vol * zirc_rho) * dose) / avogadro   # in mol   \n",
    "\n",
    "#vol_He = (1e15 * (vol * zirc_rho) * alpha) / avogadro  # in fmol\n",
    "\n",
    "#vol_He_sci = \"{:e}\".format(vol_He)\n",
    "#print(vol_He_sci)\n",
    "\n",
    "print(He_mol_total, ' mol')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Import JSON file\n",
    "Use the script below to read a JSON from a local file. Simply input the name of the JSON file as `filename.json`. IMPORTANT: the file needs to be in the same directory as this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "He_data = pd.read_json('_diff_cell_He_data.json', orient='index')\n",
    "\n",
    "# adding column name to the respective columns\n",
    "He_data.reset_index(inplace=True)  # assigning simple index for ease of calculation\n",
    "He_data.columns =['sample_ID' , 'ncc_4He' , 'ncc_4He_1sig' , 'mol_He' , 'mol_He_1sig' , 'notes' , 'temp_degC', 'time_s', 't_err']\n",
    "\n",
    "He_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Appending temperature, time, and cumulative time\n",
    "\n",
    "The easiest way would be to compile the heating temperature and time in an excel file (with column labels) and appending it to the data frame.\n",
    "\n",
    "If temp and time were input in `He_line_calc_diff_cell.ipynb` script, then now just need to calculate `t_cum`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#adding new columns from excel\n",
    "#temp_time = pd.read_excel('G168_heating_tT.xlsx')\n",
    "\n",
    "#He_data['temp_degC'] = temp_time ['temp_degC'].values\n",
    "#He_data['time_s'] = temp_time ['time_s'].values\n",
    "\n",
    "#adding new column for cumulative time in s\n",
    "He_data['t_cum'] = He_data['time_s'].cumsum()\n",
    "\n",
    "He_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Making an Arrhenius plot\n",
    "#### Calculating the fraction released `frac`, cumulative fraction released `f_cum`, and ln(D/${a^2}$) `ln_D_a2` for all the heating steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#empty lists that we'll append to\n",
    "frac = []\n",
    "ln_D_a2 = []\n",
    "\n",
    "#calculating just the fraction released at individual steps\n",
    "\n",
    "for i in range(len(He_data)):\n",
    "    frac = He_data.mol_He / He_mol_total\n",
    "\n",
    "He_data['frac'] = frac\n",
    "\n",
    "#calculating cumulative fraction released and putting it into a new column\n",
    "He_data['f_cum'] = He_data['frac'].cumsum()\n",
    "\n",
    "# first define the first step with f_cum that is more or equal to 0.45\n",
    "fcum_45 = \n",
    "\n",
    "#now calculating ln(D/a^2)\n",
    "for i in range(len(He_data)):\n",
    "\n",
    "    if i == 0:\n",
    "        ln_Da2 = np.log((He_data.f_cum[i] ** 2) * (np.pi / (4 * He_data.time_s[i])))\n",
    "    \n",
    "    elif 0 < i < fcum_45:\n",
    "        ln_Da2 = np.log((np.pi * (He_data.f_cum[i]**2 - He_data.f_cum[i-1]**2)) / (4 * He_data.time_s[i]))\n",
    "        \n",
    "    else:\n",
    "        ln_Da2 = np.log((-4 / (np.pi**2 * He_data.time_s[i])) * np.log((1 - He_data.f_cum[i])/(1 - He_data.f_cum[i-1])))\n",
    "        \n",
    "    ln_D_a2.append(ln_Da2) \n",
    "\n",
    "#adding list into dataframe\n",
    "He_data['ln_D_a2'] = ln_D_a2\n",
    "\n",
    "He_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculating uncertainties for `D_a2`.\n",
    "\n",
    "We will use equations 20 and 24 of Ginster (2018, dissertation) to calculate the uncertainty for `D_a2`.\n",
    "\n",
    "First, calculating error with respect to uncertainty in evolved gas quantities `mol_He_1sig`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D_a2_sig_X = []\n",
    "\n",
    "#now calculating uncertainty for ln(D/a^2)\n",
    "#but first define if you want to calc using 1 or 2 sig uncertainty of the evolved gas mol\n",
    "sig = 1\n",
    "\n",
    "for i in range(len(He_data)):\n",
    "\n",
    "    if i == 0:\n",
    "        long_eq = ((1 - He_data.f_cum[i]) * (sig*He_data.mol_He_1sig[i]))**2 + (He_data.f_cum[i]**2) \\\n",
    "                    * np.sum((sig*He_data.mol_He_1sig[1:])**2)\n",
    "        Da2_sig_X = (np.pi/(2 * He_data.time_s[i])) * (He_data.f_cum[i]/He_mol_total) * np.sqrt(long_eq)\n",
    "    \n",
    "    elif 0 < i < fcum_45:\n",
    "        long_eq1 = ((He_data.f_cum[i] * (1 - He_data.f_cum[i]) - He_data.f_cum[i-1] * (1 - He_data.f_cum[i-1]))**2) \\\n",
    "                    * np.sum((sig*He_data.mol_He_1sig[:i])**2)\n",
    "        long_eq2 = ((He_data.f_cum[i] * (1 - He_data.f_cum[i]) + He_data.f_cum[i-1]**2)**2) \\\n",
    "                    * ((sig*He_data.mol_He_1sig[i])**2)\n",
    "        long_eq3 = ((He_data.f_cum[i-1]**2 - He_data.f_cum[i]**2)**2) \\\n",
    "                    * np.sum((sig*He_data.mol_He_1sig[i+1:])**2)\n",
    "        Da2_sig_X = (np.pi/(2 * He_data.time_s[i] * He_mol_total)) * np.sqrt(long_eq1 + long_eq2 + long_eq3)\n",
    "        \n",
    "    else:\n",
    "        long_eq4 = ((1 + (He_data.f_cum[i-1] / (1 - He_data.f_cum[i-1])))**2) * (sig*He_data.mol_He_1sig[i])**2\n",
    "        long_eq5 = ((He_data.f_cum[i-1]/(1 - He_data.f_cum[i-1]) - He_data.f_cum[i]/(1-He_data.f_cum[i]))**2) \\\n",
    "                    * np.sum((sig*He_data.mol_He_1sig[i+1:]**2))\n",
    "        Da2_sig_X = (4 / (np.pi * He_data.time_s[i] * He_mol_total)) * np.sqrt(long_eq4 + long_eq5)\n",
    "        \n",
    "    D_a2_sig_X.append(Da2_sig_X) \n",
    "\n",
    "#adding list into dataframe\n",
    "He_data['D_a2_sig_X'] = D_a2_sig_X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll calculate error in `D_a2` wrt to `t_err`, `D_a2_sig_t`. For now we will use the cooling time as the `t_err`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "D_a2_sig_t = []\n",
    "\n",
    "#now calculating uncertainty for ln(D/a^2)\n",
    "\n",
    "for i in range(len(He_data)):\n",
    "\n",
    "    if i == 0:\n",
    "        Da2_sig_t = np.absolute((np.pi * He_data.f_cum[i]**2 / 4) * (-1 * He_data.t_err[i] / He_data.time_s[i]**2))\n",
    "    \n",
    "    elif 0 < i < fcum_45:\n",
    "        Da2_sig_t = np.absolute((np.pi * (He_data.f_cum[i]**2 - He_data.f_cum[i-1]**2) / 4) \n",
    "                                * (-1 * He_data.t_err[i] / He_data.time_s[i]**2))\n",
    "        \n",
    "    else:\n",
    "        Da2_sig_t = np.absolute((-4 / np.pi**2 * np.log((1 - He_data.f_cum[i])/(1 - He_data.f_cum[i-1]))) \n",
    "                                * (-1 * He_data.t_err[i] / He_data.time_s[i]**2)) \n",
    "        \n",
    "    D_a2_sig_t.append(Da2_sig_t) \n",
    "\n",
    "#adding list into dataframe\n",
    "He_data['D_a2_sig_t'] = D_a2_sig_t\n",
    "\n",
    "He_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now calculating error in D/a2 propogated from both `mol_He_1sig` and `t_err`, `D_a2_sig`, and calculating `ln_D_a2_sig`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "D_a2_sig = []\n",
    "ln_D_a2_sig = []\n",
    "\n",
    "for i in range(len(He_data)):\n",
    "    \n",
    "    Da2_sig = np.sqrt(He_data.D_a2_sig_X[i]**2 + He_data.D_a2_sig_t[i]**2)\n",
    "    ln_Da2_sig = np.log(np.exp(He_data.ln_D_a2[i]) + Da2_sig) - He_data.ln_D_a2[i]\n",
    "    \n",
    "    D_a2_sig.append(Da2_sig)\n",
    "    ln_D_a2_sig.append(ln_Da2_sig)\n",
    "    \n",
    "He_data['D_a2_sig'] = D_a2_sig\n",
    "He_data['ln_D_a2_sig'] = ln_D_a2_sig\n",
    "\n",
    "He_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now plotting `ln_D_a2` vs 1e4/T(K)\n",
    "Temperature is listed as degree C, so we'd need to convert that and divide it by 1e4. In order to be able to easily calculate linear regression and visualize linear vs non-linear trends, we will plot these separately. Therefore, the index number for the last NAB step should be changed for each experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (6,4), dpi=200)\n",
    "plt.xlabel(r'$10^{4}$/T ($K^{-1}$)')\n",
    "plt.ylabel(r'ln(D/$a^{2}$)')\n",
    "plt.xlim(8, 26)\n",
    "plt.ylim(-28, -8)\n",
    "plt.xticks(fontsize=10)\n",
    "plt.yticks(fontsize=10)\n",
    "plt.title(' Arrhenius plot')\n",
    "\n",
    "# FIRST OF ALL input the index for the FIRST Arrhenius datapoint that will be used for regression\n",
    "\n",
    "arr_1 = \n",
    "arr_limit = \n",
    "\n",
    "# now we will make the Arrhenius plot. NOTE last step (index 20) is not included\n",
    "\n",
    "x_axis = 1e4 / (He_data.temp_degC + 273.15)\n",
    "x_axis1 = 1e4 / (He_data.temp_degC.iloc[:arr_1] + 273.15)\n",
    "x_axis2 = 1e4 / (He_data.temp_degC.iloc[arr_1:arr_limit] + 273.15)\n",
    "x_axis3 = 1e4 / (He_data.temp_degC.iloc[arr_limit:] + 273.15)\n",
    "\n",
    "# doing this for connecting lines mainly\n",
    "plt.plot(x_axis, He_data.ln_D_a2, color = 'gray', linewidth=0.5, zorder=0)\n",
    "\n",
    "# plotting all errorbars\n",
    "plt.errorbar(x_axis, He_data.ln_D_a2, yerr=He_data.ln_D_a2_sig, fmt='none', ecolor='black', capsize=6, mew=0.5, zorder=0)\n",
    "\n",
    "# now platting scatter points\n",
    "plt.scatter(x_axis1, He_data.ln_D_a2.iloc[:arr_1], marker='o', color='none', edgecolor='black', size=60, alpha=0.75)\n",
    "plt.scatter(x_axis2, He_data.ln_D_a2.iloc[arr_1:arr_limit], marker='o', color='black', edgecolor='black', size=60, alpha=0.75)\n",
    "plt.scatter(x_axis3, He_data.ln_D_a2.iloc[arr_limit:], marker='o', color='none', edgecolor='black', size=60, alpha=0.75)\n",
    "\n",
    "\n",
    "#plt.savefig('.pdf') \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 4. Doing weighted and unweighted regression for Arrhennius plot\n",
    "First we'll set up the linear regression equations\n",
    "\n",
    "DON'T FORGET TO CHANGE THE n-value to the number of datapoints to regress!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#THIS IS FOR UNWEIGHTED LINEAR REGRESSION\n",
    "#insert n-value here. will be used for the rest of the calc\n",
    "n = \n",
    "\n",
    "def ord_lin_regress(X_vals, Y_vals):\n",
    "    \n",
    "    #hint: it helps to create arrays for XY values and X^2 (use x**2 to raise something to a power) from your input\n",
    "    #first create some empty arrays\n",
    "    \n",
    "    XY_vals = []\n",
    "    X2_vals = []\n",
    "    \n",
    "    #then try this for loop (remember that Python uses zero indexing):\n",
    "    \n",
    "    for i in range(len(X_vals)):\n",
    "        XY = X_vals[i] * Y_vals[i]\n",
    "        X2 = X_vals[i]**2\n",
    "        \n",
    "        XY_vals.append(XY)\n",
    "        X2_vals.append(X2)\n",
    "        \n",
    "    #now calculate slope and intercept below using the sum function, I'll let you try the math here\n",
    "    \n",
    "    slope = (sum(XY_vals) - (sum(X_vals)*sum(Y_vals)/n)) / (sum(X2_vals) - ((sum(X_vals)**2)/n))\n",
    "    intercept = (sum(X_vals) * sum(XY_vals) - sum(Y_vals) * sum(X2_vals)) / ((sum(X_vals)**2) - n*sum(X2_vals))\n",
    "    \n",
    "    return slope, intercept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll calculate the slope and intercept for the Arrhenius trend. To make things easier we'll make a new dataframe only with the data we want to calculate the linear regression with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#THIS IS FOR UNWEIGHTED LINEAR REGRESSION\n",
    "\n",
    "#making new lists. don't forget to change the range of datapoints\n",
    "#first define the first index after the last data for regression, or basically the last index of data for regression plus 1\n",
    "arr_limit = \n",
    "\n",
    "x_axis2 = 1e4 / (He_data.temp_degC.iloc[arr_1:arr_limit] + 273.15)\n",
    "ln_D_a2_arr = He_data.ln_D_a2.iloc[arr_1:arr_limit]\n",
    "ln_D_a2_sig_arr = He_data.ln_D_a2_sig.iloc[arr_1:arr_limit]\n",
    "\n",
    "arrhenius_data = {'temp': x_axis2, 'ln_D_a2_arr':ln_D_a2_arr, 'ln_D_a2_sig_arr':ln_D_a2_sig_arr}\n",
    "#making new dataframe and resetting index\n",
    "arr_data = pd.DataFrame(arrhenius_data)\n",
    "arr_data.reset_index(inplace=True)\n",
    "\n",
    "#now we'll do the calculation\n",
    "arrhenius = ord_lin_regress(arr_data.temp, arr_data.ln_D_a2_arr)\n",
    "\n",
    "slope_un = arrhenius[0]\n",
    "intercept_un = arrhenius[1]\n",
    "\n",
    "print('The unweighted slope is ', arrhenius[0])\n",
    "print('The unweighted intercept is ', arrhenius[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.5 Weighted linear regression with errors in y (`D_a2`)\n",
    "Alternatively, we can calculate the slope and intercept using weighted linear regression, taking into account the errors in `ln_D_a2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#first calculate the weighing factor, w that will be used for \n",
    "w_i = []\n",
    "\n",
    "for i in range(len(arr_data)):\n",
    "    w = (n * arr_data.ln_D_a2_sig_arr[i]**(-2)) / np.sum(arr_data.ln_D_a2_sig_arr**(-2))\n",
    "    w_i.append(w)\n",
    "\n",
    "arr_data['w_i'] = w_i\n",
    "\n",
    "w_sum = np.sum(arr_data.w_i)\n",
    "print(w_sum)  # check that this is equal to the number of datapoints to regress, n\n",
    "arr_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wy_i = []\n",
    "wx_i = []\n",
    "wxy_i = []\n",
    "w_x2_i = []  # w * x^2\n",
    "    \n",
    "#then try this for loop (remember that Python uses zero indexing):\n",
    "    \n",
    "for i in range(len(arr_data)):\n",
    "    wy = arr_data.w_i[i] * arr_data.ln_D_a2_arr[i]\n",
    "    wx = arr_data.w_i[i] * arr_data.temp[i]\n",
    "    wxy = arr_data.w_i[i] * arr_data.temp[i] * arr_data.ln_D_a2_arr[i]\n",
    "    w_x2 = arr_data.w_i[i] * arr_data.temp[i]**2\n",
    "        \n",
    "    wy_i.append(wy)\n",
    "    wx_i.append(wx)\n",
    "    wxy_i.append(wxy)\n",
    "    w_x2_i.append(w_x2)\n",
    "        \n",
    "#now calculate slope and intercept below using the sum function, I'll let you try the math here\n",
    "    \n",
    "slope_wght = (n * sum(wxy_i) - sum(wx_i) * sum(wy_i)) / (n * sum(w_x2_i) - sum(wx_i)**2)\n",
    "intercept_wght = (sum(wy_i) - slope_wght * sum(wx_i)) / n\n",
    "\n",
    "print('The weighted slope is ', slope_wght)\n",
    "print('The weighted intercept is ', intercept_wght)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Calculting std dev and 95% confidence interval in `slope` and `intercept` of linear regressions.\n",
    "We will use the equation in \"Chemometrics using R: Chapter 8.1\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Uncertainty for unweighted regression\n",
    "First, we will calculate the standard deviation about the regression, `s_r`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_min_y_2 =[]\n",
    "\n",
    "#calculating std dev about the regression, s_r. uncomment the eq you want to use\n",
    "for i in range(len(arr_data)):\n",
    "    #reg_eq = slope_un * arr_data.temp[i] + intercept_un  # unweighted linear eq\n",
    "    reg_eq = slope_wght * arr_data.temp[i] + intercept_wght  # weighted linear eq\n",
    "    yi_min_yp_2 = (arr_data.ln_D_a2_arr[i] - reg_eq) ** 2\n",
    "    \n",
    "    y_min_y_2.append(yi_min_yp_2)\n",
    "    \n",
    "s_r = np.sqrt(np.sum(y_min_y_2) / (n - 2))  # same n as in section 4\n",
    "print('Standard deviation about the regression =', s_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will calculate the standard deviations for the slope, `s_slope`, and y-intercept, `s_int`. IMPORTANT: Make sure to choose which equation to use for `reg_u` above!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#first to calculate the mean of x (temp)\n",
    "mean_T = np.mean(arr_data.temp)\n",
    "\n",
    "x2 = []\n",
    "x_min_x2 = []\n",
    "\n",
    "# now to calculate x_i - x_expected in a for loop\n",
    "for i in range(len(arr_data)):\n",
    "    \n",
    "    x2_i = arr_data.temp[i] ** 2\n",
    "    x_min_x2_i = (arr_data.temp[i] - mean_T) ** 2\n",
    "    \n",
    "    x2.append(x2_i)\n",
    "    x_min_x2.append(x_min_x2_i)\n",
    "    \n",
    "#calculating std dev for slope and intercept\n",
    "s_slope = np.sqrt((s_r ** 2) / np.sum(x_min_x2))\n",
    "s_int = np.sqrt((s_r**2 * sum(x2)) / (n * sum(x_min_x2)))\n",
    "\n",
    "print(s_slope, s_int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate the confidence interval, we need to select the t-test critical values `alpha_t` and degrees of freedom `df` (simply n-2). Make sure to use the two-tailed `alpha_t`. The t value can be found online (e.g., https://www.medcalc.org/manual/t-distribution-table.php)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we want to calculate 95% confidence interval (alpha_T = 0.05)\n",
    "t_test =   # for n =  (df = n - 2)\n",
    "\n",
    "CI95_slope = s_slope * t_test\n",
    "CI95_int = s_int * t_test\n",
    "\n",
    "print(CI95_slope, CI95_int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Calculate diffusion kinetic parameters  `E_a` and `D_0` with uncertainties.\n",
    "Now we can calculate the diffusion kinetic parameters: activation energy `E_a` and frequency factor `D_0`.\n",
    "#### NOTE: only either the weighted or unweighted uncertainties can be reported at a time due to the way the code is set up in part 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#calculating activation energy and frequency factor with UNWEIGHTED lin reg\n",
    "\n",
    "E_a = - gas_constant * slope_un * 10  # in kJ/mol\n",
    "D_0 = np.exp(intercept_un) * half_width**2\n",
    "\n",
    "#uncertainty calculation.\n",
    "\n",
    "E_a_stdv = E_a - (- gas_constant * (slope_un+s_slope) * 10)\n",
    "E_a_CI95 = E_a - (- gas_constant * (slope_un+CI95_slope) * 10)\n",
    "D_0_stdev = np.exp(intercept_un+s_int) * half_width**2 - D_0\n",
    "D_0_CI95 = np.exp(intercept_un+CI95_int) * half_width**2 - D_0\n",
    "\n",
    "print(u'Unweighted E_a =', E_a, '\\u00B1', E_a_stdv, 'kJ/mol (1\\u03C3).')\n",
    "print(u'Unweighted E_a =', E_a, '\\u00B1', E_a_CI95, 'kJ/mol (95% CI).')\n",
    "print(u'Unweighted D_0 =', D_0, '\\u00B1', D_0_stdev, 'cm^2/s (1\\u03C3).')\n",
    "print(u'Unweighted D_0 =', D_0, '\\u00B1', D_0_CI95, 'cm^2/s (95% CI).')\n",
    "print('These are not the actual errors. Please see note for part 5 above.') #uncomment if necessary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculating activation energy and frequency factor with WEIGHTED lin reg\n",
    "\n",
    "Ea_w = - gas_constant * slope_wght * 10  # in kJ/mol\n",
    "D0_w = np.exp(intercept_wght) * half_width**2\n",
    "\n",
    "#uncertainty calculation.\n",
    "\n",
    "Ea_w_stdv = Ea_w - (- gas_constant * (slope_wght+s_slope) * 10)\n",
    "Ea_w_CI95 = Ea_w - (- gas_constant * (slope_wght+CI95_slope) * 10)\n",
    "D0_w_stdev = np.exp(intercept_wght+s_int) * half_width**2 - D0_w\n",
    "D0_w_CI95 = np.exp(intercept_wght+CI95_int) * half_width**2 - D0_w\n",
    "\n",
    "print(u'Weighted E_a =', Ea_w, '\\u00B1', Ea_w_stdv, 'kJ/mol (1\\u03C3).')\n",
    "print(u'Weighted E_a =', Ea_w, '\\u00B1', Ea_w_CI95, 'kJ/mol (95% CI).')\n",
    "print(u'Weighted D_0 =', D0_w, '\\u00B1', D0_w_stdev, 'cm^2/s (1\\u03C3).')\n",
    "print(u'Weighted D_0 =', D0_w, '\\u00B1', D0_w_CI95, 'cm^2/s (95% CI).')\n",
    "#print('These are not the actual errors. Please see note for part 5 above.') #uncomment if necessary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Recalculating diffusion kinetic parameters `E_a` and `D_0` and their errors.\n",
    "This time we'll use the compact equations 13(a-b) in York et al. (2004) to calculate the erros in the slope and intercept (and thus `E_a` and `D_0`. This is a little more elaborate than the method in section 5. We'll also take into account the errors in the x-axis `1e4/T`. We'll use the unweighted slope as the initial value of b."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# making a new dataframe for this calculation\n",
    "aliquot_ID = He_data.sample_ID.iloc[arr_1:arr_limit]\n",
    "\n",
    "York_arr = pd.DataFrame()\n",
    "York_arr['aliquot_ID'] = aliquot_ID\n",
    "York_arr['temp'] = arr_data['temp'].values\n",
    "York_arr['ln_D_a2_arr'] = arr_data['ln_D_a2_arr'].values\n",
    "York_arr['ln_D_a2_sig_arr'] = arr_data['ln_D_a2_sig_arr'].values\n",
    "\n",
    "# restarting index at 0 and deleting extra \"Index\" column\n",
    "York_arr.reset_index(inplace=True)\n",
    "York_arr = York_arr.drop(labels=\"index\", axis=1)\n",
    "\n",
    "# adding column for error in T from excel sheet (see spreadsheet for details)\n",
    "temp_sig = pd.read_excel('_temp_err.xlsx')\n",
    "York_arr.insert(loc = 3, column = 'temp_err', value = temp_sig.T_err)\n",
    "#York_arr['temp_err'] = temp_sig['T_err'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this next cell we'll calculate the weights for the X and Y vals (step 2), then do an iterative calculation for steps 3 - 5. This is done in a while loop (gay gods please help me)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# first we'll define the slope, b and correlation `r_i`\n",
    "slope_york = slope_un\n",
    "r_i = 0  # no correlation between errors in x and y\n",
    "    \n",
    "# determining weights `w_x` and `w_y` and adding it to the dataframe\n",
    "\n",
    "w_X = []\n",
    "w_Y = []\n",
    "\n",
    "for i in range(len(York_arr)):\n",
    "    \n",
    "    weight_X = 1 / (York_arr.temp_err[i])**2\n",
    "    weight_Y = 1 / (York_arr.ln_D_a2_sig_arr[i])**2\n",
    "    \n",
    "    w_X.append(weight_X)\n",
    "    w_Y.append(weight_Y)\n",
    "    \n",
    "York_arr['w_X'] = w_X\n",
    "York_arr['w_Y'] = w_Y\n",
    "\n",
    "# let's see if this works\n",
    "\n",
    "tolerance = 1e-10  # input any number bigger than or equal to the desired tolerance\n",
    "    \n",
    "while tolerance >= 1e-15:\n",
    "    \n",
    "    W_i =[]\n",
    "    U_i = []\n",
    "    V_i = []\n",
    "    beta_i = []    \n",
    "    slope_guess = slope_york\n",
    "\n",
    "    for i in range(len(York_arr)):\n",
    "        alpha_i = np.sqrt(York_arr.w_X[i] * York_arr.w_Y[i])\n",
    "        weight_W = (York_arr.w_X[i] * York_arr.w_Y[i]) / (York_arr.w_X[i] + (slope_guess**2) \n",
    "                                        * York_arr.w_Y[i] - (2 * slope_guess * r_i * alpha_i))\n",
    "        W_i.append(weight_W)\n",
    "        \n",
    "    York_arr['W_i'] = W_i\n",
    "        \n",
    "    for i in range(len(York_arr)):\n",
    "        \n",
    "        X_mean = np.sum(York_arr.W_i * York_arr.temp) / np.sum(York_arr.W_i)\n",
    "        Y_mean = np.sum(York_arr.W_i * York_arr.ln_D_a2_arr) / np.sum(York_arr.W_i)\n",
    "    \n",
    "        U_vals = York_arr.temp[i] - X_mean\n",
    "        V_vals = York_arr.ln_D_a2_arr[i] - Y_mean\n",
    "        \n",
    "        beta_vals = York_arr.W_i[i] * (U_vals/York_arr.w_Y[i] + (slope_guess * V_vals)/York_arr.w_X[i] \n",
    "                            - (slope_guess*U_vals + V_vals) * (r_i / alpha_i))\n",
    "        \n",
    "        U_i.append(U_vals)\n",
    "        V_i.append(V_vals)\n",
    "        beta_i.append(beta_vals)\n",
    "     \n",
    "    York_arr['U_i'] = U_i\n",
    "    York_arr['V_i'] = V_i\n",
    "    York_arr['beta_i'] = beta_i\n",
    "    \n",
    "    slope_york = np.sum(York_arr.W_i * York_arr.beta_i * York_arr.V_i) / np.sum(York_arr.W_i * York_arr.beta_i * York_arr.U_i)\n",
    "    tolerance = np.abs(slope_guess - slope_york)\n",
    "\n",
    "\n",
    "York_arr['W_i'] = W_i   \n",
    "York_arr['U_i'] = U_i\n",
    "York_arr['V_i'] = V_i\n",
    "York_arr['beta_i'] = beta_i\n",
    "\n",
    "print('Slope = ', slope_york, slope_guess)\n",
    "print('Tolerance = ', tolerance)\n",
    "\n",
    "print('The weighted mean of X is', X_mean)\n",
    "print('The weighted mean of Y is', Y_mean)\n",
    "York_arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now calculating the y-intercept `int_york`\n",
    "intercept_york = Y_mean - slope_york * X_mean\n",
    "\n",
    "# now calculating new E_a and D_0\n",
    "Ea_york = - gas_constant * slope_york * 10  # in kJ/mol\n",
    "D0_york = np.exp(intercept_york) * half_width**2\n",
    "\n",
    "print('Intercept =', intercept_york)\n",
    "print('E_a =', Ea_york, 'kJ/mol')\n",
    "print('D_0 =', D0_york, 'cm^2/s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculatiing adjusted values `adj_x_i`\n",
    "x_adj = []\n",
    "\n",
    "for i in range(len(York_arr)):\n",
    "    adj_x_i = X_mean + York_arr.beta_i[i]\n",
    "    \n",
    "    x_adj.append(adj_x_i)\n",
    "\n",
    "York_arr['x_adj'] = x_adj\n",
    "\n",
    "# calculating adjusted mean of x `mean_x_adj`\n",
    "mean_x_adj = np.sum(York_arr.W_i * York_arr.x_adj) / np.sum(York_arr.W_i)\n",
    "\n",
    "# calculating adjusted u `u_adj`\n",
    "u_adj = []\n",
    "\n",
    "for i in range(len(York_arr)):\n",
    "    adj_u_i = York_arr.x_adj[i] - mean_x_adj\n",
    "    \n",
    "    u_adj.append(adj_u_i)\n",
    "    \n",
    "York_arr['u_adj'] = u_adj\n",
    "\n",
    "# calculating error in slope `sig_slope_york` and intercept `sig_int_york`\n",
    "sig_slope_york = np.sqrt(1 / np.sum(York_arr.W_i * (York_arr.u_adj**2)))\n",
    "sig_int_york = np.sqrt((1 / np.sum(York_arr.W_i)) + (mean_x_adj**2) * (sig_slope_york**2))\n",
    "\n",
    "# finally calculating error in `Ea_york` and `D0_york`\n",
    "Ea_york_sig = Ea_york - (- gas_constant * (slope_york+sig_slope_york) * 10)\n",
    "D0_york_sig_max = (np.exp(intercept_york+sig_int_york) * half_width**2) - D0_york\n",
    "D0_york_sig_min = D0_york - (np.exp(intercept_york-sig_int_york) * half_width**2)\n",
    "\n",
    "print(dose, 'alphas/g')\n",
    "print(u'Weighted (with EVA) E_a =', Ea_york, '\\u00B1', Ea_york_sig, 'kJ/mol (1\\u03C3).')\n",
    "print(u'Weighted (with EVA) D_0 =', D0_york, '\\u00B1', '(-', D0_york_sig_min, ', +',  D0_york_sig_max, ') cm^2/s (1\\u03C3).')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Arrhenius plot with trend\n",
    "Here we will replot the Arrhenius plot with the unweighted/weighted linear reg trend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (6, 5), dpi=200)\n",
    "plt.xlabel(r'$10^{4}$/T ($K^{-1}$)')\n",
    "plt.ylabel(r'ln(D/$a^{2}$)')\n",
    "plt.xlim(12, 26)\n",
    "plt.ylim(-28, -10)\n",
    "plt.xticks(fontsize=10)\n",
    "plt.yticks(fontsize=10)\n",
    "#plt.title('G3 Arrhenius plot')\n",
    "plt.rcParams['axes.linewidth'] = 0.9\n",
    "\n",
    "# now we will make the Arrhenius plot. NOTE last step (index 20) is not included\n",
    "\n",
    "x_axis1 = 1e4 / (He_data.temp_degC.iloc[:arr_1] + 273.15)\n",
    "x_axis2 = 1e4 / (He_data.temp_degC.iloc[arr_1:arr_limit] + 273.15)\n",
    "x_axis_btwn = 1e4 / (He_data.temp_degC.iloc[arr_1-1:arr_1+1] + 273.15)\n",
    "x_regtrend = np.arange(12,28) # defining range of x for regression trend\n",
    "\n",
    "plt.plot(x_axis1, He_data.ln_D_a2.iloc[:arr_1], color = 'gray', linewidth=0.75)\n",
    "plt.plot(x_axis2, He_data.ln_D_a2.iloc[arr_1:arr_limit], color = 'gray', linewidth=0.75)\n",
    "plt.plot(x_axis_btwn, He_data.ln_D_a2.iloc[arr_1-1:arr_1+1], color = 'gray', linewidth=0.75)\n",
    "\n",
    "# overlapping bigger datapoints\n",
    "#plt.scatter(x_axis1, He_data.ln_D_a2.iloc[:arr_1], color = 'b')\n",
    "#plt.scatter(x_axis2, He_data.ln_D_a2.iloc[arr_1:], color = 'm')\n",
    "\n",
    "# plotting error bars with option to magnify errors by f_err degrees\n",
    "f_err = \n",
    "\n",
    "plt.errorbar(x_axis1, He_data.ln_D_a2.iloc[:arr_1], yerr=f_err*He_data.ln_D_a2_sig.iloc[:arr_1], fmt='o', alpha=1,\n",
    "             markeredgecolor='royalblue', markerfacecolor='white',markersize=12, color='royalblue', capsize=6)\n",
    "plt.errorbar(x_axis2, He_data.ln_D_a2.iloc[arr_1:arr_limit], yerr=f_err*He_data.ln_D_a2_sig.iloc[arr_1:arr_limit], \n",
    "             markeredgecolor='brown', markerfacecolor='firebrick', markersize=12, color='brown', alpha=0.75,\n",
    "             xerr=York_arr.temp_err, fmt='o', capsize=6)\n",
    "\n",
    "# finally plotting\n",
    "reg_wght = slope_wght * x_regtrend + intercept_wght\n",
    "reg_unwght = slope_un * x_regtrend + intercept_un\n",
    "reg_york = slope_york * x_regtrend + intercept_york\n",
    "#plt.plot(x_regtrend, reg_wght, color = 'blue')\n",
    "#plt.plot(x_regtrend, reg_unwght, color = 'green')\n",
    "plt.plot(x_regtrend, reg_york, color = 'darkslategray', linestyle = '--', zorder=5, alpha=0.6)\n",
    "\n",
    "# putting text and assigning values in the desired form\n",
    "def as_si(x, ndp):\n",
    "    s = '{x:0.{ndp:d}e}'.format(x=x, ndp=ndp)\n",
    "    m, e = s.split('e')\n",
    "    return r'{m:s}\\times 10^{{{e:d}}}'.format(m=m, e=int(e))\n",
    "\n",
    "textstr = '\\n'.join((\n",
    "    'Sample G3',\n",
    "    r'$\\alpha$-dose = ${0:s}$'.format(as_si(dose_eq2,2)),\n",
    "    '$\\mathrm{E_a}=%.2f$'%(Ea_york),\n",
    "    'D$_0$ = ${0:s}$'.format(as_si(D0_york,2))))\n",
    "\n",
    "plt.text(12.6,-27,textstr,fontsize=12)\n",
    "\n",
    "plt.savefig('G3_Arrhenius_plot.pdf', transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9. Calculating chi-squared and reduced chi-squared\n",
    "This is a simple test to evaluate the goodness of fit of the regression. First, we'll calculate the `chi_sq` and then we'll calculate the `red_chisq`. This method takes into account each datapoint's error, and assumes that the errors are Gaussian. Basically we are maximizing the likelihood function.\n",
    "\n",
    "To assess reduced chi-squared:\n",
    "1. If `red_chisq` > 1, it is considered a \"bad\" fit\n",
    "2. If `red_chisq` < 1, it is considered an overfit\n",
    "3. The closer `red_chisq` is to 1, the better the fit\n",
    "\n",
    "Additionally, we can also \"improve\" the reduced chi-squared (to 1) by calibrating the errors, but this is a problem for future me in part 9 of this very lengthy notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chisq = [] # making a list of individual terms\n",
    "\n",
    "# note that `reg_eq` is either the weighted/unweighted function. see part 5\n",
    "\n",
    "for i in range(len(arr_data)):\n",
    "    unwght_reg = slope_un * arr_data.temp[i] + intercept_un  # unweighted reg eq\n",
    "    wght_reg = slope_wght * arr_data.temp[i] + intercept_wght  # weighter reg eq\n",
    "    chi = ((arr_data.ln_D_a2_arr[i] - wght_reg)**2 / (f_err*arr_data.ln_D_a2_sig_arr[i])**2)\n",
    "    \n",
    "    chisq.append(chi)\n",
    "    \n",
    "# calculating chi-squared\n",
    "chi_sq = np.sum(chisq)\n",
    "\n",
    "# now calculating reduced chi-squared, red_chisq\n",
    "red_chisq = chi_sq / (n - 2)\n",
    "\n",
    "red_chisq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating reduced chi squared for the weighted least square solution with EVA (ie, following York et al. (2004))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# revise this to account for error in x\n",
    "\n",
    "chisq_york = [] # making a list of individual terms\n",
    "\n",
    "for i in range(len(York_arr)):\n",
    "    york_reg = slope_york * York_arr.temp[i] + intercept_york  # weighted reg eq\n",
    "    chi_york = ((York_arr.ln_D_a2_arr[i] - york_reg)**2 / (f_err*York_arr.ln_D_a2_sig_arr[i])**2)\n",
    "    \n",
    "    chisq_york.append(chi_york)\n",
    "    \n",
    "# calculating chi-squared\n",
    "chi_sq_york = np.sum(chisq_york)\n",
    "\n",
    "# now calculating reduced chi-squared, red_chisq\n",
    "red_chisq_york = chi_sq_york / (n - 2)\n",
    "\n",
    "red_chisq_york"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 10. Exporting data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exporting dataframe as an excel file\n",
    "He_data.to_excel(\"_Arr_data.xlsx\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### END (FOR NOW)\n",
    "misc stuff below"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some random but useful lines of code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adding new empty columns\n",
    "#He_data.insert(loc=4, column='temp_degC', value='')\n",
    "\n",
    "#appending new values one by one OR to replace an existing one\n",
    "#He_data.at['RS23_G168_200_1', 'temp_degC'] = 200\n",
    "#He_data.at['RS23_G168_200_1', 'time_s'] = 3600\n",
    "\n",
    "#if you need to delete a row/column for some reason\n",
    "#He_data = He_data.drop(labels=\"time_s\", axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating fraction and ln(D/${a^2}$) for the first heating step (not needed anymore but just keeping it here)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculating first fraction, f_first and ln(D/a^2)\n",
    "\n",
    "#f_first = He_data.mol_He[0] / He_mol_total\n",
    "#ln_D_a2 = np.log((f_first**2) * (np.pi / (4 * He_data.time_s[0])))\n",
    "\n",
    "#print(f_first, ln_D_a2)\n",
    "\n",
    "#adding columns for frac, f_cum, and ln_D_a2\n",
    "#He_data['frac'] = ''\n",
    "#He_data['f_cum'] = ''\n",
    "#He_data['ln_D_a2'] = ''\n",
    "\n",
    "#appending to first row of dataframe\n",
    "\n",
    "#He_data.at['RS23_G168_150_1', 'frac'] = f_first\n",
    "#He_data.at['RS23_G168_150_1', 'ln_D_a2'] = ln_D_a2\n",
    "\n",
    "#He_data"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}