# MATLAB scripts used to analyze luminescence data
*This code is provided for the purposes of review and is not licensed for any other use or distribution.*

## Primary script (*.m) files
- **monteCarlo.m:** A script which loads dose response and fading measurements, then calls **main_MC.m** to iterate *numMCruns* (line 17) times. The outputs are an array of predicted ages (*age_MCarray*) and the greatest resolvable age (*age2D0*) which are saved in a workspace named **ageMCarrays.mat**. This is done for all positions and may take some time.

- **main_LnTn.m:** A script which fits the *L<sub>N</sub> / T<sub>N</sub>* profile of known-age sample KAY-110 using **profile.m** for a range of *&mu;* (mm<sup>-1</sup>) and *&sigma;&phi;<sub>0</sub>* (ka<sup>-1</sup>) values. Generated plots show the best-fit profile compared to the *L<sub>N</sub> / T<sub>N</sub>* measurement, the misfit over the range of simulated parameter values, and the predicted bleaching profiles over various hypothetical geologic sunlight exposures for sample 19-KP-H1, assuming the kinetic values determined from sample KAY-110.

- **main_MC_LT.m:** A function based on **main_MC.m** (below) which yields an array of predicted *L<sub>N</sub> / T<sub>N</sub>* ratios in field saturation (rather than fading-corrected ages). This was used to calculate the field saturation lines in main text Figure 2.

- **doseRateVsDepth.m:** A script which calculates and plots the environmental dose rate at different depths beneath the rock-ice interface. The results are plotted with grey boxes outlining the sampled intervals and the different line colors indicating different dose rate components as labelled in Figure S4.

## Dependency script (*.m) files
- **main_MC.m:** A function which inputs dose response and fading measurements, the beta source strength (Gy/s), the geologic dose rate and error (Gy/ka), a plot title string, and the number of Monte Carlo iterations and outputs the array of fading-corrected ages, given the uncertainties on kinetic parameters (see the *Demonstrating field saturation* section of the IRSL supplementary methods and lines 116 - 129 for more details). 

- **rhoPrimeFitter.m:** A function which fits signal fading measurements to equation 5 of Kars et al. (2008) to yield an estimate of *&rho;'*.

- **DRCfitter_GOK.m:** A function which fits dose response measurements to a general-order kinetics expression, equation 11 of Guralnik et al. (2015).

- **naturalDRC_GOK.m:** A function which simulates fractional saturation accumulation for a specified duration (s) given values of geologic dose rate (Gy/ka), *&rho;'*, *D<sub>0</sub>* (Gy), and kinetic order.

- **profile.m:** A function which simulates *L<sub>N</sub> / T<sub>N</sub>* at a given exposure time and depth beneath the surface using equation 4 of Sohbati et al. (2018).

## Workspace (*.mat) files
- **drcData.mat:** Contains *L<sub>x</sub> / T<sub>x</sub>* values (column 1) and 1&sigma; errors (column 2) for all aliquot positions (1, 4, 7, 9, 11, 14, 16, 19, 21, 24, 27, 30) and the IR<sub>50</sub> and pIRIR<sub>225</sub> signals. For example, the 9-by-2 matrix named **p9_50** contains the IR<sub>50</sub> dose response measurements from position 9. The regenerative doses are listed in **monteCarlo.m.** All aliquots are from sample 19-KP-H1. Positions *1, 4, 7* are slice count **1**; positions *9, 11, 14* are slice count **4**; positions *16, 19, 21* are slice count **5**; position *24* is slice count **18**; position *27* is slice count **19**; and position *30* is slice count **21**. See Sheet 1 of Data S2 of for the depth interval of each slice count.

- **fadingData.mat:** Contains *L<sub>x</sub> / T<sub>x</sub>* values (column 1) and 1&sigma; errors (column 2) for all aliquot positions and the IR<sub>50</sub> and pIRIR<sub>225</sub> signals (see **drcData.mat** for position details). Also contains the effective fading time (s) for each position and signal combination.
