---
title: "Geochron figures"
author: "Miriam Primus"
date: "2025-04-11"
output: html_document
---

# **Script for making geochronology figures for my thesis project using the IsoplotR package**

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## ***Load packages and set up the document:***

```{r}
source("setup_geochron.R")
```

```{r}
# load the IsoplotR package:  
library(IsoplotR)  
```

## ***Load in geochronology data using the read.data() isoplotR function:***

```{r}
#X=07/35, err[X], Y=06/38, err[Y], rho[X,Y] #CSV format = 1
ShC1 <- IsoplotR::read.data("csv/ShC1_UPb_Z_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 1, #X=07/35, err[X], Y=06/38, err[Y], rho[X,Y]
                            ierr = 2) #2 sigma errors


b233 <- IsoplotR::read.data("csv/233b_age_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 1, #X=07/35, err[X], Y=06/38, err[Y], rho[X,Y]
                            ierr = 2) #2 sigma errors

SLM64 <- IsoplotR::read.data("csv/SLM_DM64_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 2, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z])
                            ierr = 2) #2 sigma errors

```

## ***Calculate ages and make plots of U-Pb zircon datasets:***

### 23-6-ShC1

```{r}
#ShC1 concordia
pdf("pdf/ShC1_concordia.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
ShC1_plot <- concordia(
                x = ShC1,
                type = 1,  # Wetherill (207Pb/235U vs. 206Pb/238U)
                xlim = c(0.001, 0.06),       # Adjust x-axis limits
                ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
                show.age = 0, # don;t show age
                show.numbers = FALSE,        # Suppress grain labels
                ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
                ellipse.stroke = "black",    # Ellipse outline color
                concordia.col = "darksalmon",# Concordia curve color
                exterr = FALSE,               # Include external uncertainties
                oerr = 2,                    # Error format (2 = ±2σ)
                sigdig = 2,                  # Number of significant digits
                common.Pb = 0,               # No common Pb correction
                ticks = 5,                   # Number of tick marks on concordia
                cutoff.disc = discfilter()  # Discordance filter
)
dev.off()
```

```{r}
#ShC1 WM
pdf("pdf/ShC1_WM.pdf", width = 8, height = 8)
par(pty = "s")  # 's' means square plotting region
ShC1_WM <- weightedmean(ShC1,
             random.effects = FALSE, 
             detect.outliers = TRUE,
             plot = TRUE, 
             from = NA, 
             to = NA, 
             levels = NA, 
             clabel = "", 
             rect.col = c("#00FF0080", "#FF000080"), 
             outlier.col = "#00FFFF80", 
             sigdig = 2, 
             type = 2, 
             #cutoff.76 = 1100, 
             oerr = 2, 
             cutoff.disc = discfilter(), 
             exterr = FALSE, 
             ranked = TRUE, 
             common.Pb = 0, 
             hide = NULL, 
             omit = NULL, 
             omit.col = NA)
dev.off()
```

### 24-5-MP233b

```{r}
#233b concordia all data
pdf("pdf/233b_concordia_all.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
b233_plot <- concordia(
              x = b233,
              type = 1,  # Wetherill (207Pb/235U vs. 206Pb/238U)
              #xlim = c(0.025, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              cutoff.disc = discfilter()  # Discordance filter
)
dev.off()
```

```{r}
#233b concordia reduced
pdf("pdf/233b_concordia_reduced.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
b233_plot_reduced <- concordia(
              x = b233,
              type = 1,  # Wetherill (207Pb/235U vs. 206Pb/238U)
              xlim = c(0.020, 0.05),       # Adjust x-axis limits
              ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              anchor = 0, 
              cutoff.disc = discfilter(option = 0),  # Discordance filter
              hide = c(4, 7, 13, 14, 22, 26, 27, 33, 36, 37, 41, 21, 24), 
              omit = NULL, 
              omit.fill = NA, 
              omit.stroke = NA)
dev.off()
```

```{r}
#233b WM
pdf("pdf/233b_WM.pdf", width = 8, height = 8)
par(pty = "s")  # 's' means square plotting region
b233_WM <- weightedmean( x = b233,
             random.effects = FALSE, 
             detect.outliers = TRUE,
             plot = TRUE, 
             from = NA, 
             to = NA, 
             levels = NA, 
             clabel = "", 
             rect.col = c("#00FF0080", "#FF000080"), 
             outlier.col = "#00FFFF80", 
             sigdig = 1, 
             type = 2, 
             #cutoff.76 = 1100, 
             oerr = 2, 
             cutoff.disc = discfilter(option = 0), 
             exterr = FALSE, 
             ranked = TRUE, 
             common.Pb = 0, 
             hide = c(4, 7, 13, 14, 22, 26, 27, 33, 36, 37, 41, 21, 24), 
             omit = NA, 
             omit.col = NA)
dev.off()


```
### SLM_DM64 Apatites

```{r}
#SLM64 concordia all data
pdf("pdf/SLM64_concordia_all.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
SLM64_plot <- concordia(
              x = SLM64,
              type = 2,  # T-W
              #xlim = c(0.025, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 2, # discordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              cutoff.disc = discfilter()  # Discordance filter
)
dev.off()
```


### 90913-3A - Chlorite Quartz mylonite

```{r}
a90913_3A <- IsoplotR::read.data("csv/90913-3A_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 3, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z])
                            ierr = 2) #2 sigma errors

```

```{r}
#90913-3A concordia all data
#pdf("pdf/90913_3A_concordia_all.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a90913_3A_plot <- concordia(
              x = a90913_3A,
              type = 1,  # Wetherill (207Pb/235U vs. 206Pb/238U)
              #xlim = c(0.025, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              cutoff.disc = discfilter()  # Discordance filter
)
#dev.off()
```

```{r}
#90913-3A concordia reduced
#pdf("pdf/90913_3A_concordia_reduced.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a90913_3A_plot_reduced <- concordia(
              x = a90913_3A,
              type = 1,  # Wetherill (207Pb/235U vs. 206Pb/238U)
              #xlim = c(0.020, 0.05),       # Adjust x-axis limits
             #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              anchor = 0, 
              #cutoff.disc = dscf,  # Discordance filter
              hide = c(3, 4, 5, 13, 14, 15, 19, 20, 21, 22, 23, 24), #Need to change which ones to hide
              omit = NULL, 
              omit.fill = NA, 
              omit.stroke = NA)
#dev.off()
```


```{r}
##90913-3A WM
#pdf("pdf/90913_3A_WM.pdf", width = 8, height = 8)
par(pty = "s")  # 's' means square plotting region
a90913_3A_WM <- weightedmean( x = a90913_3A,
             random.effects = FALSE, 
             detect.outliers = TRUE,
             plot = TRUE, 
             from = NA, 
             to = NA, 
             levels = NA, 
             clabel = "", 
             rect.col = c("#00FF0080", "#FF000080"), 
             outlier.col = "#00FFFF80", 
             sigdig = 1, 
             type = 2, 
             #cutoff.76 = 1100, 
             oerr = 2, 
             #cutoff.disc = dscf, 
             exterr = FALSE, 
             ranked = TRUE, 
             common.Pb = 0, 
             hide = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 21, 22, 23, 24), 
             omit = NA, 
             omit.col = NA)
#dev.off()
```

### 91013-5A - Rhyolitic porphyry, altered

```{r}
a91013_5A <- IsoplotR::read.data("csv/91013-5A_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 3, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z])
                            ierr = 2) #2 sigma errors

```

```{r}
#91013-5A concordia all data
pdf("pdf/91013_5A_concordia_all.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a91013_5A_plot <- concordia(
              x = a91013_5A,
              type = 1,  # Wetherill (207Pb/235U vs. 206Pb/238U)
              #xlim = c(0.025, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              cutoff.disc = discfilter()  # Discordance filter
)
dev.off()
```
```{r}
#91013_5A concordia reduced
pdf("pdf/91013_5A_concordia_reduced.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a91013_5A_plot_reduced <- concordia(
              x = a91013_5A,
              type = 1,  # Wetherill (207Pb/235U vs. 206Pb/238U)
              #xlim = c(0.020, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              anchor = 0, 
              cutoff.disc = dscf,  # Discordance filter
              hide = c(6, 11, 26, 27, 34, 35), #Need to change which ones to hide
              omit = NULL, 
              omit.fill = NA, 
              omit.stroke = NA)
dev.off()
```

```{r}
#91013_5A WM
pdf("pdf/a91013_5A_WM.pdf", width = 8, height = 8)
par(pty = "s")  # 's' means square plotting region
a91013_5A_WM <- weightedmean( x = a91013_5A,
             random.effects = FALSE, 
             detect.outliers = TRUE,
             plot = TRUE, 
             from = NA, 
             to = NA, 
             levels = NA, 
             clabel = "", 
             rect.col = c("#00FF0080", "#FF000080"), 
             outlier.col = "#00FFFF80", 
             sigdig = 1, 
             type = 3, 
             #cutoff.76 = 1100, 
             oerr = 2, 
             #cutoff.disc = discfilter(option = 0), 
             exterr = FALSE, 
             ranked = TRUE, 
             common.Pb = 0, 
             hide = c(6, 11, 26, 27, 34, 35), 
             omit = NA, 
             omit.col = NA)
dev.off()
```


### 91113-6B - Altered, aphanitic felsite

```{r}
a91113_6B <- IsoplotR::read.data("csv/91113-6B_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 3, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z])
                            ierr = 2) #2 sigma errors
```

```{r}
#91113-6B concordia all data

# filter data
dscf <- discfilter(option = '2', before = TRUE, cutoff = c(-5, 15))

#pdf("pdf/91113-6B_concordia_all.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a91113_6B_plot <- concordia(
              x = a91113_6B,
              type = 1,  # Wetherill (207Pb/235U vs. 206Pb/238U) because the sample is younger
              #xlim = c(0.025, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Don;t show age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                  # Number of tick marks on concordia
              cutoff.disc = dscf)  # Discordance filter

#dev.off()
```
```{r}
#a91113_6B concordia reduced - removed data points that are presumed to be inherited
#pdf("pdf/91113_6B_concordia_reduced.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a91113_6B_plot_reduced <- concordia(
              x = a91113_6B,
              type = 2,  # Wetherill (207Pb/235U vs. 206Pb/238U)
             # xlim = c(150, 250),       # Adjust x-axis limits
              #ylim = c(0.00, 0.03),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              anchor = 0, 
              #cutoff.disc = dscf,  # Discordance filter
            hide = c(1, 2, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 25),
              omit = NULL, 
              omit.fill = NA, 
              omit.stroke = NA)
#dev.off()




```

```{r}
#91113_6B WM
pdf("pdf/91113_6B_WM.pdf", width = 8, height = 8)
par(pty = "s")  # 's' means square plotting region
a91113_6B_WM <- weightedmean( x = a91113_6B,
             random.effects = FALSE, #if FALSE, attributes any excess dispersion to an underestimation of the analytical uncertainties. This akin to a ‘model-1’ isochron regression.
             detect.outliers = TRUE, #logical flag indicating whether outliers should be detected and rejected using Chauvenet’s Criterion.
             plot = TRUE, #plot the data
             from = NA, 
             to = NA, 
             levels = NA, 
             clabel = "", 
             rect.col = c("#00FF0080", "#FF000080"), 
             outlier.col = "#00FFFF80", 
             sigdig = 2, 
             type = 2, #plot the 206/238 age because this smaple is younger than 1200
             #cutoff.76 = 1100, #only used if type = 4
             oerr = 2,  # Error format (2 = ±2σ)
             #cutoff.disc = dscf,
             exterr = FALSE, #do not propagate external uncertainties
             ranked = TRUE, #plot the aliquots in order of increasing age
             common.Pb = 0, #no common lead correction
             hide = c(1, 2, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 25), #
             omit = NA, 
             omit.col = NA)
dev.off()
```



### 91213-22A - Altered, aphanitic felsite

```{r}
a91213_22A <- IsoplotR::read.data("csv/91213-22A_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 3, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z])
                            ierr = 2) #2 sigma errors
```

```{r}
#91213-22A concordia all data
#pdf("pdf/91213-22A_concordia_all.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a91213_22A_plot <- concordia(
              x = a91213_22A,
              type = 1,  # Wetherill (207Pb/235U vs. 206Pb/238U)
              #xlim = c(0.025, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              cutoff.disc = discfilter()  # Discordance filter
)
#dev.off()
```

```{r}
#a91213_22A concordia reduced
pdf("pdf/91213_22A_concordia_reduced.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a91213_22A_plot_reduced <- concordia(
              x = a91213_22A,
              type = 2, # Wetherill (207Pb/235U vs. 206Pb/238U)
              #xlim = c(0.0, 0.2),       # Adjust x-axis limits
              #ylim = c(0.00, 0.03),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              anchor = 0, 
              #cutoff.disc = dscf,  # Discordance filter
             # hide = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 22, 23, 24, 25), #Need to change which ones to hide
              #hide = c(1, 19),
              omit = NULL, 
              omit.fill = NA, 
              omit.stroke = NA)
dev.off()
```

```{r}
#91213_22A WM
pdf("pdf/91213_22A_WM.pdf", width = 8, height = 8)
par(pty = "s")  # 's' means square plotting region
a91213_22A_WM <- weightedmean( x = a91213_22A,
             random.effects = FALSE, 
             detect.outliers = TRUE,
             plot = TRUE, 
             from = NA, 
             to = NA, 
             levels = NA, 
             clabel = "", 
             rect.col = c("#00FF0080", "#FF000080"), 
             outlier.col = "#00FFFF80", 
             sigdig = 1, 
             type = 2, 
             #cutoff.76 = 1100, 
             oerr = 2, 
             #cutoff.disc = discfilter(option = 0), 
             exterr = FALSE, 
             ranked = TRUE, 
             common.Pb = 0, 
            hide = c(2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 22, 23, 24, 25), 
            # hide = c(1, 19),
             omit = NA, 
             omit.col = NA)
dev.off()
```


### 52014-12A - Augen gneiss

```{r}
a52014_12A <- IsoplotR::read.data("csv/52014-12A_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 3, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z])
                            ierr = 2) #2 sigma errors
```

```{r}
#52014-12A concordia all data
#pdf("pdf/52014-12A_concordia_all.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a52014_12A_plot <- concordia(
              x = a52014_12A,
              type = 2,  #Tera-Wasserburg– 207Pb/206Pb vs. 238U/206Pb (because sampe is old)
              #xlim = c(0.025, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              cutoff.disc = discfilter()  # Discordance filter
)
#dev.off()
```

```{r}
#52014_12A concordia reduced

#filter data
dscf5 <- discfilter(option = '2', before = TRUE, cutoff = c(-5, 15))


#pdf("pdf/52014_12A_concordia_reduced.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a52014_12A_plot_reduced <- concordia(
              x = a52014_12A,
              type = 2,  # #Tera-Wasserburg– 207Pb/206Pb vs. 238U/206Pb (because sample is old )
              #xlim = c(0.020, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # plot the data without calcualting an age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              anchor = 0, #do not anchor the isochon
              cutoff.disc = dscf5, #cutoff samples with % discordance less than -5 or greater than 20
              hide = c(9, 11, 30), #Omit samples that are not concordant
              omit = NULL, 
              omit.fill = NA, 
              omit.stroke = NA)
#dev.off()
```

```{r}
#52014_12A WM
#pdf("pdf/52014_12A_WM.pdf", width = 8, height = 8)
par(pty = "s")  # 's' means square plotting region
a52014_12A_WM <- weightedmean( x = a52014_12A,
             random.effects = FALSE, #if FALSE, attributes any excess dispersion to an underestimation of the analytical uncertainties. This akin to a ‘model-1’ isochron regression.
             detect.outliers = TRUE, #logical flag indicating whether outliers should be detected and rejected using Chauvenet’s Criterion.
             plot = TRUE, #plot the data
             from = NA, 
             to = NA, 
             levels = NA, 
             clabel = "", 
             rect.col = c("#00FF0080", "#FF000080"), 
             outlier.col = "#00FFFF80", 
             sigdig = 2, 
             type = 3, #plot the 206/207 age because this sample is older than 1200 Ma
             #cutoff.76 = 1100, #only used if type = 4
             oerr = 2,  # Error format (2 = ±2σ)
             cutoff.disc = discfilter(option = '2', before = TRUE, cutoff = c(-5, 15)), #cutoff samples with % discordance less than -5 or greater than 15
             exterr = FALSE, #do not propagate external uncertainties
             ranked = TRUE, #plot the aliquots in order of increasing age
             common.Pb = 0, #no common lead correction
             #hide = c(9, 11, 30), #
             omit = NA, 
             omit.col = NA)
#dev.off()


#type = scalar indicating whether to plot the 207Pb/235U age (type=1), the 206Pb/238U age (type=2), the 207Pb/206Pb age (type=3), the 207Pb/206Pb-206Pb/238U age (type=4), the concordia_age (type=5), or the 208Pb/232Th age (type=6).
```


### 52015-6A - Aphanitic to porphyritic felsite, likely rhyolite

```{r}
a52015_6A <- IsoplotR::read.data("csv/52015-6A_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 3, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z])
                            ierr = 2) #2 sigma errors
```

```{r}
#52015-6A concordia all data
#pdf("pdf/52015-6A_concordia_all.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a52015_6A_plot <- concordia(
              x = a52015_6A,
              type = 1,  # Wetherill (207Pb/235U vs. 206Pb/238U)
              #xlim = c(0.025, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              #cutoff.disc = discfilter()  # Discordance filter
              #hide = c(1, 17)
)
#dev.off()
```

```{r}
#52015_6A concordia reduced
#pdf("pdf/52015_6A_concordia_reduced.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a52015_6A_plot_reduced <- concordia(
              x = a52015_6A,
              type = 2,  # Wetherill (207Pb/235U vs. 206Pb/238U)
              #xlim = c(230, 250),       # Adjust x-axis limits
              #ylim = c(0.04, 0.06),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              anchor = 0,
              #cutoff.disc = discfilter(),  # Discordance filter
              hide = c(1, 17)) #Need to change which ones to hide
              #hide = c(27),
              #omit = NULL, 
              #omit.fill = NA, 
             # omit.stroke = NA)
#dev.off()
```

```{r}
#52015_6A WM
#pdf("pdf/52015_6A_WM.pdf", width = 8, height = 8)
par(pty = "s")  # 's' means square plotting region
a52015_6A_WM <- weightedmean(x = a52015_6A,
             random.effects = FALSE, 
             detect.outliers = FALSE,
             plot = TRUE, 
             from = NA, 
             to = NA, 
             levels = NA, 
             clabel = "", 
             rect.col = c("#00FF0080", "#FF000080"), 
             outlier.col = "#00FFFF80", 
             sigdig = 1, 
             type = 1, 
             #cutoff.76 = 1100, 
             oerr = 2, 
             #cutoff.disc = discfilter(), 
             exterr = FALSE, 
             ranked = TRUE, 
             common.Pb = 0, 
             hide = c(1, 17, 20, 28))
#dev.off()
```


### 91917-3A - Altered, aphanitic felsite

```{r}
a91917_3A <- IsoplotR::read.data("csv/91917-3A_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 3, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z])
                            ierr = 2) #2 sigma errors
```

```{r}
#91917-3A concordia all data
pdf("pdf/91917-3A_concordia_all.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a91917_3A_plot <- concordia(
              x = a91917_3A,
              type = 2,  # Wetherill (207Pb/235U vs. 206Pb/238U)
              #xlim = c(0.025, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              #cutoff.disc = discfilter()  # Discordance filter
)
dev.off()
```
```{r}
#91917_3A concordia reduced
pdf("pdf/91917_3A_concordia_reduced.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a91917_3A_plot_reduced <- concordia(
              x = a91917_3A,
              type = 2,  # Wetherhill
              #xlim = c(0.020, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              anchor = 0, 
              #cutoff.disc = dscf,  # Discordance filter
              hide = c(3, 5, 15, 23, 17, 8, 26, 9, 18), #Removed discordant grains and grains that do not 
              omit = NULL, 
              omit.fill = NA, 
              omit.stroke = NA)
dev.off()
```

```{r}
#91917_3A WM
pdf("pdf/91917_3A_WM.pdf", width = 8, height = 8)
par(pty = "s")  # 's' means square plotting region
a91917_3A_WM <- weightedmean( x = a91917_3A,
             random.effects = FALSE, 
             detect.outliers = TRUE,
             plot = TRUE, 
             from = NA, 
             to = NA, 
             levels = NA, 
             clabel = "", 
             rect.col = c("#00FF0080", "#FF000080"), 
             outlier.col = "#00FFFF80", 
             sigdig = 1, 
             type = 3, 
             #cutoff.76 = 1100, 
             oerr = 2, 
             #cutoff.disc = dscf, 
             exterr = FALSE, 
             ranked = TRUE, 
             common.Pb = 0, 
             hide = c(3, 5, 15, 23, 17, 8, 26, 9, 18), #hide discordant grains and grains that do not match the main population
             omit = NA, 
             omit.col = NA)
dev.off()
```


### 91917-5B - Rhyolitic porphyry, altered

```{r}
a91917_5B <- IsoplotR::read.data("csv/91917-5B_data.csv", 
                            method = "U-Pb", #U-Pb zircon
                            format = 3, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z])
                            ierr = 2) #2 sigma errors
```

```{r}
#91917-5B concordia all data
#pdf("pdf/91917-5B_concordia_all.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a91917_5B_plot <- concordia(
              x = a91917_5B,
              type = 2,  #Tera-Wasserburg– 207Pb/206Pb vs. 238U/206Pb (because sampe is old)
              #xlim = c(0.025, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              cutoff.disc = discfilter()  # Discordance filter
)
#dev.off()
```

```{r}
#91917_5B concordia reduced
#pdf("pdf/91917_5B_concordia_reduced.pdf", width = 8, height = 8) #save as PDF
par(pty = "s")  # 's' means square plotting region
a91917_5B_plot_reduced <- concordia(
              x = a91917_5B,
              type = 2,  #Tera-Wasserburg– 207Pb/206Pb vs. 238U/206Pb (because sample is old)
              #xlim = c(0.020, 0.05),       # Adjust x-axis limits
              #ylim = c(0.0035, 0.0055),     # Adjust y-axis limits
              show.age = 0, # Only show concordia age
              show.numbers = FALSE,        # Suppress grain labels
              ellipse.fill = c("#00FF0040", "#FF000040"),  # Ellipse fill colors with transparency
              ellipse.stroke = "black",    # Ellipse outline color
              concordia.col = "darksalmon",# Concordia curve color
              exterr = FALSE,               # Include external uncertainties
              oerr = 2,                     # Error format (2 = ±2σ)
              sigdig = 2,                  # Number of significant digits
              common.Pb = 0,               # No common Pb correction
              ticks = 5,                   # Number of tick marks on concordia
              anchor = 0, 
              cutoff.disc = dscf,  # Discordance filter
              hide = c(5, 9, 10), #Need to change which ones to hide
              omit = NULL, 
              omit.fill = NA, 
              omit.stroke = NA)
#dev.off()
```

```{r}
#91917_5B WM
#pdf("pdf/91917_5B_WM.pdf", width = 8, height = 8)
par(pty = "s")  # 's' means square plotting region
a91917_5B_WM <- weightedmean( x = a91917_5B,
             random.effects = FALSE, 
             detect.outliers = TRUE,
             plot = TRUE, 
             from = NA, 
             to = NA, 
             levels = NA, 
             clabel = "", 
             rect.col = c("#00FF0080", "#FF000080"), 
             outlier.col = "#00FFFF80", 
             sigdig = 1, 
             type = 3, #plot the 206/207 age because this sample is older than 1200 Ma 
             #cutoff.76 = 1100, 
             oerr = 2, 
             cutoff.disc = dscf, 
             exterr = FALSE, 
             ranked = TRUE, 
             common.Pb = 0, 
             hide = c(5, 9, 10), #hide grains that do not match the mian population 
             omit = NA, 
             omit.col = NA)
#dev.off()
```


## ***Load in detrital zircon U-Pb data:***

```{r}
#a91113DZ KDE
a91113DZ <- IsoplotR::read.data("csv/91113-13A_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 3, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z]
                            ierr = 2) #2 sigma errors

```

```{r}
#a52014DZ KDE
a52014DZ <- IsoplotR::read.data("csv/52014-14A_data.csv",
                            method = "U-Pb", #U-Pb zircon
                            format = 3, #X=07/35, err[X], Y=06/38, err[Y], Z=07/06, err[Z] (, rho[X,Y]) (, rho[Y,Z]
                            ierr = 2) #2 sigma errors 

```

## ***Make KDE plots for detrital zircon data:***

```{r}
#KDE plot of 91113-13A

#filter data
dscf1 <- discfilter(option = '2', before = TRUE, cutoff = c(-5, 20))

pdf("pdf/dz91113-13A_KDE2.pdf", width = 8, height = 5)
#par(pty = "s")  # 's' means square plotting region
kde(x = a91113DZ, 
     from = 0, 
     to = 3500, 
     bw = NA, 
     adaptive = TRUE, 
     log = FALSE, 
     n = 512, #number of segments of the density estimate
     plot = TRUE, 
     rug = TRUE, 
     xlab = "age [Ma]", 
     ylab = "", 
     kde.col = rgb(1, 0, 1, 0.6), 
     hist.col = rgb(0, 1, 0, 0.2), 
     show.hist = FALSE, 
     bty = "n", 
     binwidth = NA, 
     type = 4, 
     cutoff.76 = 1300, 
     cutoff.disc=dscf1,
     common.Pb = 0, 
     hide = NULL, 
     nmodes = 0,
     sigdig = 2)
 
dev.off()

```



```{r}
#KDE plot of 52014-14A

#filter data (option 2 is for percent discordance, BEFORE = before common lead correction, removes values that are above 20% discordant, or beloe -5% discordant)
dscf2 <- discfilter(option = '2', before = TRUE, cutoff = c(-5, 20))

pdf("pdf/dz52014-14A_KDE2.pdf", width = 8, height = 5)
#par(pty = "s")  # 's' means square plotting region 
kde(x = a52014DZ, 
     from = 0, 
     to = 3500, 
     bw = NA, 
     adaptive = TRUE, 
     log = FALSE, 
     n = 512, #number of segments of the density estimate
     plot = TRUE, 
     rug = TRUE, 
     xlab = "age [Ma]", 
     ylab = "", 
     kde.col = rgb(0, 0, 1, 0.6), 
     hist.col = rgb(0, 1, 0, 0.2), 
     show.hist = FALSE, 
     bty = "n", 
     binwidth = NA, 
     type = 4, 
     cutoff.76 = 1300, 
     cutoff.disc=dscf2,
     common.Pb = 0, 
     hide = NULL, 
     nmodes = 0,
     sigdig = 2)

dev.off()


```




## ***Bonus! Make bar chart for kinematics:***

```{r}
# Load ggplot2
library(ggplot2)

# Create data frame with kinematics counts
kinematic_data <- data.frame(
  Kinematics=c("Top-SW","Top-NE","Mix","Coaxial") ,  
  Count=c(56,12,7,2)
  )

# Barplot
kinematics_plot <- ggplot(kinematic_data, aes(x=reorder(Kinematics, -Count), y=Count, fill=Kinematics)) + 
  geom_bar(stat = "identity", color="black") +
  scale_fill_manual(values = c("black", "green", "blue","red")) +
  labs(x = "Kinematics", 
       y = "Count")+
  theme_clean()+
  theme(legend.position="none")

kinematics_plot

ggsave(plot = kinematics_plot, filename = "kinematics_plot.pdf", width = 7.5, height = 6, units = c("in"), dpi = 400)

```



## ***Bonus! Make bar chart for deformation conditions:***

```{r}
# Load ggplot2
library(ggplot2)

# Create data frame with kinematics counts
deformation_data <- data.frame(
  Deformation=c("1","2","3") ,  
  Count=c(30,32,7)
  )

# Barplot
deformation_plot <- ggplot(deformation_data, aes(x=Deformation, y=Count, fill=Deformation)) + #aes(x=reorder(Deformation, -Count),
  geom_bar(stat = "identity", color="black") +
  scale_fill_manual(values = c("pink", "violet", "gold")) +
  labs(x = "Deformation Category", 
       y = "Count")+
  theme_clean()+
  theme(legend.position="none")

deformation_plot

ggsave(plot = deformation_plot, filename = "deformation_plot.pdf", width = 7.5, height = 6, units = c("in"), dpi = 400)

```
