---
title: "Introduction to biocharkit"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to biocharkit}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(biocharkit)
```

All data in this vignette is **synthetically generated** with a fixed
random seed, purely to demonstrate the package's functions end to end. It
is not real experimental data.

## 1. Sample IDs

Pyrolysis samples are often named to encode their production conditions,
e.g. `SBC450-30` for a spirulina-derived biochar pyrolysed at 450 degC for
30 minutes. `parse_sbc_id()` splits these back out:

```{r}
ids <- c("SBC300-15", "SBC300-30", "SBC300-60",
         "SBC450-15", "SBC450-30", "SBC450-60",
         "SBC600-15", "SBC600-30", "SBC600-60")
parse_sbc_id(ids)
```

## 2. Adsorption capacity and removal efficiency

```{r}
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)
```

## 3. Isotherms: fitting and comparing multiple models

```{r}
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)
)
```

Comparing R² across models like this is the usual way to decide which
isotherm best describes a given system — here, unsurprisingly, Langmuir
and Sips (which nests Langmuir) come out on top, since the data was
simulated from a Langmuir relationship.

```{r, fig.alt = "Langmuir isotherm fit with fitted curve overlaid on observed points"}
plot_isotherm(langmuir, Ce_iso, qe_iso)
```

Confidence intervals on the fitted parameters:

```{r}
fit_confint(langmuir)
```

### Batch fitting across many samples at once

In practice you'll often have an isotherm for every sample in a factorial
design, not just one. `fit_isotherm_batch()` fits a model separately per
group and returns one tidy row per sample:

```{r}
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")
```

## 4. Kinetics

```{r}
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)
)
```

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

## 5. Thermodynamics

Given a distribution coefficient `Kc` measured at several temperatures,
`fit_vant_hoff()` recovers standard enthalpy, entropy, and Gibbs free
energy of adsorption:

```{r}
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
```

## 6. FTIR: baseline correction, peak picking, and functional groups

```{r, 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)
```

```{r}
spectrum <- baseline_correct(spectrum_raw)
find_ftir_peaks(spectrum, window = 10, min_prominence = 0.05)
functional_group_density(spectrum)
```

## 7. XRD: crystallinity index via peak deconvolution

```{r}
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
```

If you already have crystalline/amorphous peak areas from your own XRD
software, [xrd_crystallinity_index()] computes the index directly without
needing to fit anything in R.

## 8. BET surface area

```{r}
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)
```

## 9. Proximate and ultimate analysis

```{r}
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)
```

## 10. TGA: DTG, decomposition stages, and Kissinger kinetics

A single synthetic TGA curve, with a moisture step and one main
devolatilization step:

```{r}
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()` reads moisture/volatile-matter/ash off the curve at three
temperature breakpoints and hands them to `proximate_analysis()`:

```{r}
tga_stages(tga)
```

```{r, fig.alt = "TG and DTG curves against temperature"}
plot_tga(tga)
```

Kissinger kinetics needs the same decomposition step observed at several
heating rates, not a single curve; here from DTG peak temperatures
recorded across four heating rates:

```{r}
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)
```

## 11. Correlation matrix

```{r}
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")
```