## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4.5
)

## ----setup--------------------------------------------------------------------
library(biocharkit)

## -----------------------------------------------------------------------------
ids <- c("SBC300-15", "SBC300-30", "SBC300-60",
         "SBC450-15", "SBC450-30", "SBC450-60",
         "SBC600-15", "SBC600-30", "SBC600-60")
parse_sbc_id(ids)

## -----------------------------------------------------------------------------
set.seed(1)
C0 <- rep(50, 9)
temps <- rep(c(300, 450, 600), each = 3)
Ce <- pmax(1, 42 - 0.045 * temps + rnorm(9, 0, 2))

qe_batch(C0, Ce, V = 0.05, m = 0.1)
removal_efficiency(C0, Ce)

## -----------------------------------------------------------------------------
set.seed(2)
Ce_iso <- c(2, 5, 10, 20, 35, 50, 70, 90)
qe_iso <- (5.2 * 0.12 * Ce_iso) / (1 + 0.12 * Ce_iso) + rnorm(length(Ce_iso), 0, 0.08)

langmuir   <- fit_langmuir(Ce_iso, qe_iso)
freundlich <- fit_freundlich(Ce_iso, qe_iso)
temkin     <- fit_temkin(Ce_iso, qe_iso)
dr         <- fit_dr(Ce_iso, qe_iso)
sips       <- fit_sips(Ce_iso, qe_iso)

data.frame(
  model = c("Langmuir", "Freundlich", "Temkin", "Dubinin-Radushkevich", "Sips"),
  R2 = round(c(langmuir$R2, freundlich$R2, temkin$R2, dr$R2, sips$R2), 4)
)

## ----fig.alt = "Langmuir isotherm fit with fitted curve overlaid on observed points"----
plot_isotherm(langmuir, Ce_iso, qe_iso)

## -----------------------------------------------------------------------------
fit_confint(langmuir)

## -----------------------------------------------------------------------------
set.seed(3)
batch_df <- do.call(rbind, lapply(c("SBC300-30", "SBC450-30", "SBC600-30"), function(id) {
  Qmax <- runif(1, 3, 6); KL <- runif(1, 0.05, 0.2)
  Ce <- c(2, 5, 10, 20, 35, 50, 70)
  data.frame(id = id, Ce = Ce,
             qe = (Qmax * KL * Ce) / (1 + KL * Ce) + rnorm(length(Ce), 0, 0.05))
}))

fit_isotherm_batch(batch_df, "id", "Ce", "qe", model = "langmuir")

## -----------------------------------------------------------------------------
set.seed(4)
t_kin <- c(5, 15, 30, 60, 120, 180, 240, 360)
qt_kin <- (0.02 * 4.3^2 * t_kin) / (1 + 0.02 * 4.3 * t_kin) + rnorm(length(t_kin), 0, 0.05)

pfo           <- fit_pfo(t_kin, qt_kin)
pso           <- fit_pso(t_kin, qt_kin)
elovich       <- fit_elovich(t_kin, qt_kin)
intraparticle <- fit_intraparticle(t_kin, qt_kin)

data.frame(
  model = c("Pseudo-first-order", "Pseudo-second-order", "Elovich", "Intraparticle diffusion"),
  R2 = round(c(pfo$R2, pso$R2, elovich$R2, intraparticle$R2), 4)
)

## ----fig.alt = "Pseudo-second-order kinetics fit with fitted curve overlaid on observed points"----
plot_kinetics(pso, t_kin, qt_kin, model = "pso")

## -----------------------------------------------------------------------------
temperature_K <- c(298, 308, 318, 328)
Kc <- c(1.8, 2.3, 2.9, 3.4)
vh <- fit_vant_hoff(temperature_K, Kc)
vh[c("deltaH_kJmol", "deltaS_Jmolk", "R2")]
vh$table

## ----fig.alt = "Synthetic FTIR spectrum with five characteristic biochar functional-group peaks"----
wn <- seq(4000, 400, by = -4)
gaussian_peak <- function(x, center, height, width) height * exp(-((x - center)^2) / (2 * width^2))
baseline <- 0.08 + 0.00003 * (4000 - wn)
spectrum_raw <- data.frame(
  wavenumber_cm1 = wn,
  intensity = baseline +
    gaussian_peak(wn, 3400, 0.35, 120) + gaussian_peak(wn, 2920, 0.15, 40) +
    gaussian_peak(wn, 1710, 0.22, 25)  + gaussian_peak(wn, 1600, 0.28, 20) +
    gaussian_peak(wn, 1030, 0.30, 35)
)
plot_ftir_spectrum(spectrum_raw)

## -----------------------------------------------------------------------------
spectrum <- baseline_correct(spectrum_raw)
find_ftir_peaks(spectrum, window = 10, min_prominence = 0.05)
functional_group_density(spectrum)

## -----------------------------------------------------------------------------
set.seed(5)
two_theta <- seq(10, 40, by = 0.1)
intensity <- 300 * exp(-((two_theta - 22)^2) / (2 * 0.8^2)) +
             150 * exp(-((two_theta - 29)^2) / (2 * 0.6^2)) +
             rnorm(length(two_theta), 0, 3)
intensity <- pmax(0, intensity)

xrd <- xrd_deconvolve(two_theta, intensity, peak_centers = c(22, 29))
xrd$peaks
xrd$crystallinity_index

## -----------------------------------------------------------------------------
P_P0 <- c(0.05, 0.10, 0.15, 0.20, 0.25, 0.30)
Q <- c(12.1, 15.8, 18.9, 21.6, 24.1, 26.8)
bet_surface_area(P_P0, Q)

## -----------------------------------------------------------------------------
proximate_analysis(moisture_pct = 3.2, volatile_matter_pct = 28.5, ash_pct = 12.1)
ultimate_ratios(C_pct = 65.2, H_pct = 2.8, O_pct = 18.4, N_pct = 1.1)

## -----------------------------------------------------------------------------
Temp <- seq(25, 800, by = 5)
Weight <- 100 - 5 * (Temp > 110) -
  40 / (1 + exp(-(Temp - 340) / 20)) + rnorm(length(Temp), 0, 0.15)
tga <- data.frame(temperature_C = Temp, weight_pct = pmax(Weight, 20))

dtg <- tga_dtg(tga)
peaks <- find_dtg_peaks(dtg, min_prominence = 0.05)
peaks
assign_dtg_peaks(peaks$peak_temperature_C)

## -----------------------------------------------------------------------------
tga_stages(tga)

## ----fig.alt = "TG and DTG curves against temperature"------------------------
plot_tga(tga)

## -----------------------------------------------------------------------------
beta <- c(5, 10, 15, 20)          # deg C / min
Tp_C <- c(320, 335, 344, 351)     # DTG peak temperature at each beta
kis <- tga_kinetics_kissinger(beta, Tp_C)
kis[c("Ea_kJmol", "A_min1", "R2")]
fit_confint(kis)

## -----------------------------------------------------------------------------
pH <- 6.5 + 0.004 * temps + rnorm(9, 0, 0.15)
EC <- 0.8 + 0.003 * temps + rnorm(9, 0, 0.08)
df <- data.frame(temperature_C = temps, pH = pH, EC_dSm = EC)
correlation_matrix(df, method = "spearman")