In blongworth/HybridGIS: Data, Functions, and Analysis for the NOSAMS Hybrid Gas Ion Source

Figures for SNEAP talk

Setup

library(tidyverse)
library(here)
library(amstools)
library(HybridGIS)
library(patchwork)
library(broom)
library(gt)

Current vs flow

Good figure in first tests. Combine with data from dual bellows?

# factors for flow calc
r <- 2.5E-5
u <- 1.49E-5
x <- 8.2

data <- get_hgis_data(here("data/USAMS040317R.txt")) %>% 
  filter(Pos == 13)

# Assign Run numbers
counter <- 0
Run <- c(1:190)
for (i in 1:nrow(data)) {
  if (data$Meas[i] == 1) {
    counter <- counter + 1
  }
  Run[i] <- counter
} 
data$Run <- Run

# Load pressure data
pres <- read.csv(here("data/USAMS040317Pressure.csv"))
pres$Run <- as.numeric(row.names(pres))

# Join to run data
data <- inner_join(data, pres) %>%
  mutate(Source = Source * 1E-7, # I think the source pressures were to the -7, but unsure...
         flow = flowcalc(Bottle, r, u, x),
         eff = hgis_eff(flow, he12C, cs = 1)) # cs 1 to estimate LE current

# summarize
d.s <- data %>% 
  group_by(Run, Bottle, Source, flow) %>% 
  summarize(le12C = mean(le12C),
            he12C = mean(he12C),
            counts = sum(CntTotGT),
            interr = 1/sqrt(counts),
            C1412he = mean(X14.12he),
            C1412he_sd = sd(X14.12he),
            ferr = interr * C1412he,
            cor1412he = mean(cor1412he),
            eff = mean(eff))

# plot
filter(d.s, flow > 0) %>%
  ggplot(aes(flow, le12C)) +
  geom_smooth() +
  geom_point() +
  labs(x = expression(CO[2]~flow~rate~~mu~L %.% min^{-1}),
       y = bquote('12C-  ' *mu~ 'A')) +
  scale_y_reverse() +
  ylim(0,-10) +
  xlim(0,15.5)

Data from dual bellows tests

read.csv(here("data/14Nov2017LPtest.csv"), skip = 4) %>%
  mutate(flow = flowcalc(supply/10, 2E-5, 1.49E-5, 1.58)) %>% 
  ggplot(aes(flow, current)) +
  geom_smooth() + 
  geom_point() +
  labs(x = "Flow (ul/min)", y = "Source current (uA)")

Efficiency vs flow

Dual bellows

Breakseal results for OX-I, C-1, OX-II

data <- get_hgis_data(here("data/USAMS052318RHIGS.txt")) %>% 
                  mutate(Num = ifelse(Pos == 15, "S", "U"))
stdrat <- mean(data$cor1412he[data$Num == "S" & data$outlier == FALSE])
data_523 <- mutate(data, normFm = norm_gas(cor1412he, stdrat))

data <- get_hgis_data(here("data/USAMS053018R.txt"), as.Date("2018-05-30")) %>% 
           mutate(Pos = ifelse(Sample.Name == "OX-II", 6, Pos),
                  Num = ifelse(Pos == 5, "S", "U"))
stdrat <- mean(data$cor1412he[data$Num == "S"])
data_530 <- mutate(data, normFm = norm_gas(cor1412he, stdrat))

data <- get_hgis_data(here("data/USAMS060618R.txt")) %>% 
  filter(Pos != 1,                           #remove spurious pos1 runs
         !(Pos == 3 & X13.12he > 0.011)) %>%     #remove solid OX-I fliers
                  mutate(Num = ifelse(Pos == 5, "S", "U"))
stdrat <- mean(data$cor1412he[data$Num == "S" & data$outlier == FALSE])
data_606 <- mutate(data, normFm = norm_gas(cor1412he, stdrat))

dbd <- rbind(data_523, data_530, data_606)

std <- getStdTable()

dbds <- sum_hgis(dbd) %>% 
  mutate(rec_num = case_when(str_detect(Sample.Name, "Live") ~ 101730,
                             Sample.Name == "OX-I" ~ 34148,
                             Sample.Name == "OX-II" ~ 34149,
                             str_starts(Sample.Name, "C1") ~ 83028)) %>% 
 left_join(select(std, rec_num, fm_consensus), by = "rec_num") %>% 
 mutate(fm_consensus = case_when(rec_num == 101730 ~ 1.0398,
                                 rec_num == 72446 ~ 0.0013,
                                 TRUE ~ fm_consensus),
        fm_diff = mean - fm_consensus)

dbds %>% 
  filter(Cur < 10) %>% 
  ggplot(aes(fm_consensus, fm_diff)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(aes(ymin = fm_diff - merr, ymax = fm_diff + merr),
                 position = position_dodge2(width = .1),
                 size = .5) +
  labs(title = "Closed inlet")

Open split

Dillution curve

From helium_displacement.Rmd

Maybe also show long dillution runs plotted by CO2 flow?

Yes. Calculate flow at time using concCO2

data <- get_hgis_data(here("data/USAMS040121R.txt"), as.Date("2021-04-02")) %>% 
  filter(Pos == 5) %>% 
  mutate(time = cum_acqtime/60,
         co2flow = concCO2(time, r = 244) * 30 #30ul/min
         )

Sample current vs time

# flow vs time subplot
flow_time <- data.frame(x = 0:125) %>% 
  ggplot(aes(x)) +
  stat_function(fun=function(x) concCO2(x, r = 244) * 30, aes(color = "CO2")) +
  stat_function(fun=function(x) 30 - (concCO2(x, r = 244) * 30), aes(color = "Helium")) +
  scale_color_manual("Gas", values = c("green", "blue")) +
  labs(title = "Gas flows to source",
       subtitle = "7mL vial, 244uL/min helium, 30ul/min to source",
       y = "Gas flow (uL/min)") +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())

# current vs time subplot
cur_time <- ggplot(data, aes(time, he12C)) +
  geom_smooth(span = .4, se=FALSE, color = "green") +
  geom_point() +
  xlim(0, 125) +
  labs(title = "Ion current",
       x = "Time (min)",
       y = "12C current (uA)")

flow_time / cur_time

Expected CO2 ratio/flow

Plot volume of CO2 in sample vial during displacement by helium. It looks like optimal currents are maintained until about 2000s, which corresponds to about 2mL of CO2. Stable ratios until about 1.25 mL CO2. Can we back calculate flow to source if current falls off below ~2ul/min CO2?

Add flow at time to data and plot current vs flow. 30 ul/m is estimate of capillary flow based on past results.

ggplot(data, aes(co2flow, he12C)) +
  geom_smooth(span = .3, se=FALSE) +
  geom_point() +
  labs(title = "Vial dillution current",
       subtitle = "250uL/min displacement, 30ul/min delivery",
       x = "CO2 flow (ul/m)",
       y = "12C current (uA)")

Reproducibility of standards

Data from 1-08, 1-22, 1-29

Blanks - blank values

compare dead gas, vial dead gas, C1, and TIRI-F

Blanks - estimation of constant current contaminant

data<- get_hgis_data(here("data/USAMS101320R.txt"), as.Date("2020-11-17")) %>% 
  mutate(Num = ifelse(Pos %in% c(2, 4), "S", ifelse(Num == "S", "U", Num)))
stdrat <- mean(data$cor1412he[data$Num == "S" & !data$outlier])
data1013 <- data %>% mutate(normFm = norm_gas(cor1412he, stdrat))

data <- get_hgis_data(here("data/USAMS120320R.txt"), as.Date("2020-12-04")) %>% 
  mutate(Num = ifelse(Pos %in% c(2, 4), "S", ifelse(Num == "S", "U", Num)))
stdrat <- mean(data$cor1412he[data$Num == "S" & !data$outlier])
data1204 <- data %>% mutate(normFm = norm_gas(cor1412he, stdrat))

data1211 <- get_hgis_data(here("data/USAMS120320R.txt")) %>% 
  filter(Pos %in% 1:4 | as.Date(ts) == "2020-12-11")

data <- get_hgis_data(here("data/USAMS120320R.txt"), as.Date("2020-12-18")) %>% 
  filter(Pos != 0) %>% 
  mutate(Num = ifelse(Pos %in% c(2, 4), "S", ifelse(Num == "S", "U", Num)))
stdrat <- mean(data$cor1412he[data$Num == "S" & !data$outlier])
data1218 <- data %>% mutate(normFm = norm_gas(cor1412he, stdrat))

data0108 <- get_hgis_data(here("data/USAMS120320R.txt"), as.Date("2021-01-08"))

# no solid standards
# df <- get_hgis_data(here("data/USAMS120320R.txt"), as.Date("2021-01-15")) 
# stdrat <- mean(df$cor1412he[df$pos_name == "8 - LiveGasOS" & !df$outlier])
# data0115 <- mutate(df, normFm = norm_gas(cor1412he, stdrat))

data <- get_hgis_data(here("data/USAMS020521R.txt")) %>% 
  mutate(Num = ifelse(Pos %in% c(2, 4), "S", ifelse(Num == "S", "U", Num)))
stdrat <- mean(data$cor1412he[data$Num == "S" & !data$outlier])
data205 <- data %>% mutate(normFm = norm_gas(cor1412he, stdrat))

data <- get_hgis_data(here("data/USAMS030421R.txt"), as.Date("2021-03-05")) %>% 
  mutate(Num = ifelse(Pos %in% c(2, 4), "S", ifelse(Num == "S", "U", Num)))
stdrat <- mean(data$cor1412he[data$Num == "S" & !data$outlier])
data304 <- data %>% mutate(normFm = norm_gas(cor1412he, stdrat))

data <- get_hgis_data(here("data/USAMS040121R.txt"), as.Date("2021-04-09")) %>% 
  mutate(Num = ifelse(Pos %in% c(22, 24), "S", ifelse(Num == "S", "U", Num)))
stdrat <- mean(data$cor1412he[data$Num == "S" & !data$outlier])
data409 <- data %>% mutate(normFm = norm_gas(cor1412he, stdrat))

data <- get_hgis_data(here("data/USAMS041521R.txt"), as.Date("2021-04-16"))  %>% 
  mutate(Num = ifelse(Pos %in% c(22, 24), "S", ifelse(Num == "S", "U", Num)))
stdrat <- mean(data$cor1412he[data$Num == "S" & !data$outlier])
data415 <- data %>% mutate(normFm = norm_gas(cor1412he, stdrat))

std <- getStdTable()

data <- rbind(data1013, data1204, data1211, data1218, data205, data304, data409, data415) %>% 
  mutate(rec_num = case_when(str_detect(Sample.Name, "LiveGas") ~ 101730,
                             str_starts(Sample.Name, "C-1") ~ 83028,
                             str_starts(Sample.Name, "C1") ~ 83028,
                             str_starts(Sample.Name, "TIRI-F") ~ 2138,
                             str_starts(Sample.Name, "TIRI-I") ~ 17185,
                             str_starts(Sample.Name, "C-2") ~ 1082,
                             str_starts(Sample.Name, "C2") ~ 1082,
                             str_starts(Sample.Name, "NOSAMS") ~ 38809,
                             str_detect(Sample.Name, "DeadGas") ~ 72446)) %>% 
 left_join(select(std, rec_num, fm_consensus), by = "rec_num") %>% 
 mutate(fm_consensus = case_when(rec_num == 101730 ~ 1.0398,
                                 rec_num == 72446 ~ 0.0013,
                                 TRUE ~ fm_consensus))

xl <- 8

livefm_cur <- data %>% 
  filter(!(he12C < 2.5 & normFm > 1.05) & fm_consensus > 0.9) %>% 
  ggplot(aes(he12C, normFm)) + 
  geom_point() +
  xlim(0, xl) +
  ylim(.8, 1.2) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank())

deadfm_cur <- data %>% 
  filter(fm_consensus < 0.1 & normFm < 0.3 & Sample.Name != "DilDeadGas6") %>% 
  ggplot(aes(he12C, normFm)) + 
  geom_point() +
  xlim(0, xl) +
  ylim(0, .4) 

livefm_cur / deadfm_cur

Current dependent model

ds <- sum_hgis(data) %>% 
  filter((str_detect(Sample.Name, "LiveGas") & Cur > 1) |
         (str_detect(Sample.Name, "DeadGas") & Cur > 1) |
          str_detect(Sample.Name, "Blank")) %>% 
  mutate(Gas = case_when(str_detect(Sample.Name, "LiveGas") ~ "LiveGas",
                         str_detect(Sample.Name, "DeadGas") ~ "DeadGas",
                         TRUE ~ "Blank"),
         Cur_inv = 1/Cur)

# Fit linear model
fits <- ds %>% 
  ungroup() %>% 
  nest(data = -Gas) %>% 
  mutate(fit = map(data, ~lm(mean ~ Cur_inv, data = .x)),
         tidied = map(fit, tidy)) %>% 
  unnest(tidied) 

blankfit <- fits %>% 
  select(Gas, term, estimate) %>% 
  pivot_wider(names_from = c(Gas, term), values_from = estimate) %>% 
  mutate(inv_m_blank = -(`DeadGas_(Intercept)` - `LiveGas_(Intercept)`)/(DeadGas_Cur_inv - LiveGas_Cur_inv),
         Fm_blank = `DeadGas_(Intercept)` + DeadGas_Cur_inv * inv_m_blank,
         m_blank = 1/inv_m_blank) 

# Modeled Fm and mass of blank
blankfit %>% 
  select(Fm_blank, m_blank)

# Plot of data and model with measured blanks
ggplot(ds, aes(Cur_inv, mean, color = Gas)) +
  geom_abline(slope = blankfit$LiveGas_Cur_inv, 
              intercept = blankfit$`LiveGas_(Intercept)`) +
  geom_abline(slope = blankfit$DeadGas_Cur_inv, 
              intercept = blankfit$`DeadGas_(Intercept)`) +
  geom_pointrange(aes(ymin = mean - merr, ymax = mean + merr), size = .3) + 
  xlim(0, 10) +
  #ylim(0, 1.2)
  ggtitle("Fm of live and dead CO2 vs. inverse 12C current",
          subtitle = "Fit lines should cross at Fm and current of blank")

Results for carbonates

Combine all runs 2-05, 3-05, 4-09, 4-16

show as agreement with consensus vs consensus value with sample names

Need normalized values, BC values, Errors, and rec_num or consensus value

carb_data <- map_dfr(list.files(here("data_analysed"), pattern = "USAMS", full.names = TRUE), read_csv) %>% 
  filter(!is.na(rec_num),
         Cur > 2) %>% 
  mutate(fm_diff = mean - fm_consensus,
         fm_sigma = sigma(mean, fm_consensus, merr),
         fmbc_diff = fm_lbc - fm_consensus,
         fmbc_sigma = sigma(fm_lbc, fm_consensus, merr))

Current independent blank

Linear fit used to correct:

$$Fm_s = Fm_m - Fm_{blk}\frac{Fm_{std} - Fm_m}{Fm_{std}}$$

Mean Fm of blanks: r carb_data %>% filter(fm_consensus < 0.1 & Cur > 2) %>% pull(mean) %>% mean() %>% format(digits = 3)

carb_data %>% filter(fm_consensus < 0.1 &
              Cur > 2) %>% 
  ggplot(aes(Cur, mean)) +
  geom_errorbar(aes(ymin = mean - merr, ymax = mean + merr)) +
  geom_errorbarh(aes(xmin = Cur - Cur.sd, xmax = Cur + Cur.sd)) +
  geom_point(aes(color = wheel), size = 3) +
  scale_color_manual("wheel", values = c("#b7bf10", "#00a9e0", "#0069b1")) +
  labs(title = "Blanks",
       x = "12C Current (μA)",
       y = "Fraction modern") +
  theme_classic() +
  theme(legend.position = c(0.82, 0.76),
        legend.background = element_rect(fill = "white", color = "black"))

Secondary standards

Currents > 2 μA
Current dependence
Underestimation of error

Data summary

Mean difference: r format(mean(carb_data$fmbc_diff), digits = 1)
Standard deviation of difference: r format(sd(carb_data$fmbc_diff), digits = 2)

ps <- 0.3
dw <- 0.03
yl <- .1

secs <- ggplot(carb_data, aes(fm_consensus, fm_diff, color = Cur)) +
  geom_hline(yintercept = 0) +
  geom_smooth(method = "lm", se = FALSE) +
  geom_pointrange(aes(ymin = fm_diff - merr, ymax = fm_diff + merr), 
                  position = position_dodge2(width = dw),
                  size = ps,) +
  lims(x = c(-0.05, 1.1),
       y = c(-yl, yl)) +
  labs(title = "Difference from expected Fm",
       subtitle = "Uncorrected",
       colour = "Cur (μA)",
       y = "Fm difference") +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.direction = "horizontal",
        legend.position = c(.80, 1),
        legend.justification = "bottom")

secs_lbc <- ggplot(carb_data, aes(fm_consensus, fmbc_diff, color = Cur)) +
  geom_hline(yintercept = 0) +
  geom_smooth(method = "lm", se = FALSE) +
  geom_pointrange(aes(ymin = fmbc_diff - merr, ymax = fmbc_diff + merr), 
                  position = position_dodge2(width = dw),
                  size = ps) +
  lims(x = c(-0.05, 1.1),
       y = c(-yl, yl)) +
  labs(subtitle = "Blank corrected",
       x = "Fm expected",
       y = "Fm difference") +
  theme_classic() +
  theme(legend.position = "None")

secs / secs_lbc

same for sigma

ps <- 0.3
dw <- 0.03
yl <- .1

secs <- ggplot(carb_data, aes(fm_consensus, fm_sigma, color = Cur)) +
  geom_hline(yintercept = 0) +
  geom_smooth(method = "lm", se = FALSE) +
  geom_point(position = position_dodge2(width = dw)) +
  #lims(x = c(-0.05, 1.1),
  #     y = c(-yl, yl)) +
  labs(title = "Sigma vs. expected Fm",
       subtitle = "Uncorrected",
       colour = "Cur (μA)",
       y = "Sigma") +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.direction = "horizontal",
        legend.position = c(.80, 1),
        legend.justification = "bottom")

secs_lbc <- ggplot(carb_data, aes(fm_consensus, fmbc_sigma, color = Cur)) +
  geom_hline(yintercept = 0) +
  geom_smooth(method = "lm", se = FALSE) +
  geom_point(position = position_dodge2(width = dw)) +
  #lims(x = c(-0.05, 1.1),
  #     y = c(-yl, yl)) +
  labs(subtitle = "Blank corrected",
       x = "Fm expected",
       y = "Sigma") +
  theme_classic() +
  theme(legend.position = "None")

secs / secs_lbc

ggplot(carb_data, aes(fm_consensus, sigma)) +
  geom_hline(yintercept = 0) +
  geom_point() +
  ggtitle("Sigma of LBC-corrected vs consensus")

Data summary

Mean difference: r format(mean(carb_data$fmbc_diff), digits = 1)
Standard deviation of difference: r format(sd(carb_data$fmbc_diff), digits = 2)
Number of samples: r length(carb_data$fmbc_diff)

Data table

carb_data_sum <- carb_data %>% 
  select(-dil_factor, -sd, -exterr, -interr, -diff, -normFm, -sigma) %>% 
  rename(norm_ratio = mean,
         sig_norm_ratio = merr)
write_csv(carb_data_sum, here("data_analysed/carb_data_summary.csv"))

Carb data summary by type

carb_data_sum %>% 
  mutate(Name = case_when(rec_num == 101730 ~ "LiveGas",
                          rec_num == 83028 ~ "C-1",
                          rec_num == 2138 ~ "TIRI-F",
                          rec_num == 17185 ~ "TIRI-I",
                          rec_num == 1082 ~ "C-2",
                          rec_num == 38809 ~ "NOSAMS2",
                          rec_num == 72446 ~ "DeadGas")) %>% 
  select(rec_num, Name, norm_ratio, fm_lbc, sig_norm_ratio) %>% 
  group_by(Name) %>% 
  summarize(across(c(norm_ratio, fm_lbc, sig_norm_ratio), c(mean = mean, sd = sd)), N = n()) %>%
  arrange(norm_ratio_mean) %>% 
  gt() %>% 
  fmt_number(columns = 2:7,
             decimals = 4)

blongworth/HybridGIS documentation built on Dec. 19, 2021, 10:41 a.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

blongworth/HybridGIS
Data, Functions, and Analysis for the NOSAMS Hybrid Gas Ion Source

In blongworth/HybridGIS: Data, Functions, and Analysis for the NOSAMS Hybrid Gas Ion Source

Setup

Current vs flow

Efficiency vs flow

Dual bellows

Open split

Dillution curve

Expected CO2 ratio/flow

Reproducibility of standards

Blanks - blank values

Blanks - estimation of constant current contaminant

Current dependent model

Results for carbonates

Current independent blank

Linear fit used to correct:

Secondary standards

Data summary

same for sigma

Data summary

Data table

R Package Documentation

Browse R Packages

We want your feedback!

blongworth/HybridGIS Data, Functions, and Analysis for the NOSAMS Hybrid Gas Ion Source

In blongworth/HybridGIS: Data, Functions, and Analysis for the NOSAMS Hybrid Gas Ion Source

Setup

Current vs flow

Efficiency vs flow

Dual bellows

Open split

Dillution curve

Expected CO2 ratio/flow

Reproducibility of standards

Blanks - blank values

Blanks - estimation of constant current contaminant

Current dependent model

Results for carbonates

Current independent blank

Linear fit used to correct:

Secondary standards

Data summary

same for sigma

Data summary

Data table

R Package Documentation

Browse R Packages

We want your feedback!

blongworth/HybridGIS
Data, Functions, and Analysis for the NOSAMS Hybrid Gas Ion Source