| Title: | Biochar Characterisation and Adsorption Data Analysis |
|---|---|
| Description: | A toolkit for analysing biochar characterisation and batch adsorption experiments. Provides functions to parse structured sample identifiers encoding pyrolysis conditions, read raw FTIR and XRD instrument output, compute adsorption capacity and removal efficiency, fit adsorption isotherms following Langmuir (1918) <doi:10.1021/ja02242a004> and Sips (1948) <doi:10.1063/1.1746922> among other models, fit adsorption kinetics following Ho and McKay (1999) <doi:10.1016/S0032-9592(98)00112-5> and Chien and Clayton (1980) <doi:10.2136/sssaj1980.03615995004400020013x> among other models, fit batches of samples at once, compute van't Hoff thermodynamic parameters, baseline-correct and pick peaks in FTIR spectra, deconvolve XRD patterns into a crystallinity index, compute BET surface area following Brunauer, Emmett, and Teller (1938) <doi:10.1021/ja01269a023>, compute proximate and ultimate analysis summaries including directly from a thermogravimetric analysis (TGA) curve, compute a smoothed derivative thermogravimetric (DTG) curve and pick its decomposition peaks, fit non-isothermal decomposition kinetics from multi-heating-rate TGA data following Kissinger (1957) <doi:10.1021/ac60131a045>, build correlation matrices with p-values, and produce publication-style base-graphics figures including 600 dpi TIFF export. Built on base R ('stats', 'graphics', 'grDevices') so it has no dependency on packages that require external CRAN network access to install. |
| Authors: | Sukamal Sarkar [aut, cre] |
| Maintainer: | Sukamal Sarkar <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.3.0 |
| Built: | 2026-07-21 14:46:03 UTC |
| Source: | https://github.com/cran/biocharkit |
Matches a vector of observed DTG peak temperatures (e.g. from
find_dtg_peaks()) against a stage-range reference table. Where a
temperature falls inside more than one range (lignin decomposition
overlaps the others, since it is broad and underlying), the first
matching row of reference wins, so the default reference lists the
narrower hemicellulose/cellulose ranges before the broad lignin range.
assign_dtg_peaks( peak_temperatures_C, reference = tga_decomposition_reference() )assign_dtg_peaks( peak_temperatures_C, reference = tga_decomposition_reference() )
peak_temperatures_C |
Numeric vector of observed DTG peak temperatures (degrees C). |
reference |
A reference table as returned by
|
A data frame with columns peak_temperature_C, stage,
description. Peaks matching no range get NA stage/description.
assign_dtg_peaks(c(80, 300, 350))assign_dtg_peaks(c(80, 300, 350))
Matches a vector of observed peak positions (e.g. peak-picked maxima) against a band-assignment reference table.
assign_ftir_peaks(peaks, reference = ftir_band_reference())assign_ftir_peaks(peaks, reference = ftir_band_reference())
peaks |
Numeric vector of observed peak wavenumbers (cm^-1). |
reference |
A reference table as returned by
|
A data.frame with columns peak_cm1, group, assignment.
Peaks matching no range in the reference get NA group/assignment.
assign_ftir_peaks(c(3420, 2920, 1710, 1600, 1030))assign_ftir_peaks(c(3420, 2920, 1710, 1600, 1030))
Removes a simple estimated baseline from a spectrum before further
analysis (e.g. functional_group_density() or find_ftir_peaks()).
Two methods are provided: "linear" subtracts the straight line
connecting the spectrum's first and last points (a common quick
correction for a roughly flat, tilted baseline); "rolling_min"
subtracts a rolling-minimum envelope, which follows a more uneven
baseline at the cost of also shaving the base of broad peaks.
baseline_correct(spectrum, method = c("linear", "rolling_min"), window = 25)baseline_correct(spectrum, method = c("linear", "rolling_min"), window = 25)
spectrum |
A data frame with columns |
method |
|
window |
For |
A data frame in the same shape as spectrum, with intensity
replaced by the baseline-corrected values (never negative; clipped
at 0).
spec <- data.frame(wavenumber_cm1 = seq(4000, 400, by = -4), intensity = 0.1 + runif(901, 0, 0.05)) corrected <- baseline_correct(spec)spec <- data.frame(wavenumber_cm1 = seq(4000, 400, by = -4), intensity = 0.1 + runif(901, 0, 0.05)) corrected <- baseline_correct(spec)
Fits the linearised multi-point BET (Brunauer-Emmett-Teller) equation
1 / (Q * (P0/P - 1)) = 1/(Qm*C) + ((C-1)/(Qm*C)) * (P/P0) by ordinary
least squares, then converts the BET monolayer capacity Qm to a
specific surface area using the adsorbate's molecular cross-sectional
area.
bet_surface_area( P_P0, Q, cross_sectional_area_nm2 = 0.162, molar_volume_cm3mol = 22414 )bet_surface_area( P_P0, Q, cross_sectional_area_nm2 = 0.162, molar_volume_cm3mol = 22414 )
P_P0 |
Numeric vector of relative pressures (P/P0, dimensionless, between 0 and 1). |
Q |
Numeric vector of quantity adsorbed at each |
cross_sectional_area_nm2 |
Molecular cross-sectional area of the
adsorbate, in nm^2. Default |
molar_volume_cm3mol |
Molar volume of an ideal gas at STP, in
cm^3/mol. Default |
Conventionally only points in the relative pressure range
P/P0 roughly 0.05-0.35 are used (where the BET model is valid); this
function does not enforce that range for you, so filter P_P0 and Q
before calling if your data extends beyond it.
A list with elements:
Qm_cm3g |
BET monolayer capacity (cm^3/g STP) |
C |
BET constant (related to adsorbate-adsorbent interaction energy) |
SBET_m2g |
Specific surface area (m^2/g) |
R2 |
Coefficient of determination of the linear BET fit |
model |
The underlying |
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)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)
Computes a pairwise correlation matrix (Spearman by default) across all
numeric columns of a data frame, along with a matching matrix of
p-values from stats::cor.test(). Implemented with base R only (no
Hmisc dependency).
correlation_matrix( df, method = c("spearman", "pearson", "kendall"), use = "pairwise.complete.obs" )correlation_matrix( df, method = c("spearman", "pearson", "kendall"), use = "pairwise.complete.obs" )
df |
A |
method |
Correlation method passed to |
use |
Passed to |
A list with elements estimate (correlation coefficient
matrix) and p_value (matching p-value matrix), both with row/column
names equal to the numeric column names of df.
df <- data.frame(temp = c(300, 450, 600, 300, 450, 600), EC = c(1.1, 1.8, 2.6, 1.0, 1.9, 2.5), pH = c(7.1, 8.0, 9.2, 7.0, 8.1, 9.0)) correlation_matrix(df)df <- data.frame(temp = c(300, 450, 600, 300, 450, 600), EC = c(1.1, 1.8, 2.6, 1.0, 1.9, 2.5), pH = c(7.1, 8.0, 9.2, 7.0, 8.1, 9.0)) correlation_matrix(df)
Finds local maxima in dtg$dtg (weight-loss-rate peaks): a point is a
candidate peak if it is the highest point within window indices on
either side, following the same local-maxima-plus-prominence approach
as find_ftir_peaks(). Peaks with non-positive dtg (baseline drift
or mass gain, e.g. oxidation) are dropped.
find_dtg_peaks(dtg, window = 5, min_prominence = 0)find_dtg_peaks(dtg, window = 5, min_prominence = 0)
dtg |
A data frame with columns |
window |
Half-width (in points) used to test local maxima.
Default |
min_prominence |
Minimum prominence (in %/degC) for a local
maximum to be kept. Default |
A data frame with columns peak_temperature_C, dtg_max,
prominence, sorted by decreasing dtg_max.
set.seed(1) T <- seq(25, 800, by = 5) W <- 100 - 5 * (T > 110) - 40 / (1 + exp(-(T - 340) / 20)) + rnorm(length(T), 0, 0.15) tga <- data.frame(temperature_C = T, weight_pct = pmax(W, 20)) find_dtg_peaks(tga_dtg(tga))set.seed(1) T <- seq(25, 800, by = 5) W <- 100 - 5 * (T > 110) - 40 / (1 + exp(-(T - 340) / 20)) + rnorm(length(T), 0, 0.15) tga <- data.frame(temperature_C = T, weight_pct = pmax(W, 20)) find_dtg_peaks(tga_dtg(tga))
Finds local maxima in spectrum$intensity: a point is a candidate peak
if it is the highest point within window indices on either side. Weak
shoulders are then filtered out by a minimum prominence requirement
(height above the higher of the two valleys flanking the peak).
Intended as a quick first pass, not a substitute for visual inspection
of the spectrum.
find_ftir_peaks(spectrum, window = 5, min_prominence = 0)find_ftir_peaks(spectrum, window = 5, min_prominence = 0)
spectrum |
A data frame with columns |
window |
Half-width (in points) used to test local maxima.
Default |
min_prominence |
Minimum prominence (in intensity units) for a
local maximum to be kept. Default |
A data frame with columns peak_cm1, intensity,
prominence, sorted by decreasing intensity.
wn <- seq(4000, 400, by = -2) spec <- data.frame( wavenumber_cm1 = wn, intensity = exp(-((wn - 1710)^2) / (2 * 20^2)) + exp(-((wn - 2920)^2) / (2 * 30^2)) ) find_ftir_peaks(spec)wn <- seq(4000, 400, by = -2) spec <- data.frame( wavenumber_cm1 = wn, intensity = exp(-((wn - 1710)^2) / (2 * 20^2)) + exp(-((wn - 2920)^2) / (2 * 30^2)) ) find_ftir_peaks(spec)
Wraps stats::confint() around the underlying nls/lm fit object
produced by this package's fit_*() functions, returning a tidy data
frame rather than a bare matrix.
fit_confint(fit, level = 0.95)fit_confint(fit, level = 0.95)
fit |
A fit object returned by any of |
level |
Confidence level. Default |
A data frame with columns parameter, estimate, lower,
upper.
Ce <- c(2, 5, 10, 20, 35, 50, 70) qe <- c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, 4.4) fit <- fit_langmuir(Ce, qe) fit_confint(fit)Ce <- c(2, 5, 10, 20, 35, 50, 70) qe <- c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, 4.4) fit <- fit_langmuir(Ce, qe) fit_confint(fit)
Fits qe = qs * exp(-Kad * epsilon^2), where
epsilon = R * temperature_K * log(1 + 1/Ce) is the Polanyi potential.
Also returns the mean free energy of adsorption
E = 1 / sqrt(2 * Kad) (kJ/mol), conventionally used to distinguish
physisorption (E < 8 kJ/mol) from chemisorption (8-16 kJ/mol).
fit_dr(Ce, qe, temperature_K = 298.15, start = NULL)fit_dr(Ce, qe, temperature_K = 298.15, start = NULL)
Ce |
Numeric vector of equilibrium concentrations (mg/L). |
qe |
Numeric vector of equilibrium adsorption capacities (mg/g). |
temperature_K |
Absolute temperature at which the isotherm was
measured, in Kelvin. Default |
start |
Optional named list |
A list with elements qs (mg/g), Kad (mol^2/kJ^2), E
(kJ/mol), R2, model.
Ce <- c(2, 5, 10, 20, 35, 50, 70) qe <- c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, 4.4) fit_dr(Ce, qe)Ce <- c(2, 5, 10, 20, 35, 50, 70) qe <- c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, 4.4) fit_dr(Ce, qe)
Fits qt = (1/beta) * log(1 + alpha * beta * t), where alpha is the
initial sorption rate (mg/g/min) and beta relates to the extent of
surface coverage / activation energy for chemisorption (g/mg).
fit_elovich(t, qt, start = NULL)fit_elovich(t, qt, start = NULL)
t |
Numeric vector of contact times. |
qt |
Numeric vector of adsorption capacity at time |
start |
Optional named list |
A list with elements alpha, beta, R2, model.
t <- c(5, 15, 30, 60, 120, 240) qt <- c(1.2, 2.1, 2.9, 3.6, 4.0, 4.2) fit_elovich(t, qt)t <- c(5, 15, 30, 60, 120, 240) qt <- c(1.2, 2.1, 2.9, 3.6, 4.0, 4.2) fit_elovich(t, qt)
Fits qe = Kf * Ce^(1/n) by nonlinear least squares, with a fallback to
the linearised form log(qe) = log(Kf) + (1/n) * log(Ce) fit by OLS if
the nonlinear fit fails to converge.
fit_freundlich(Ce, qe, start = NULL)fit_freundlich(Ce, qe, start = NULL)
Ce |
Numeric vector of equilibrium concentrations (mg/L). |
qe |
Numeric vector of equilibrium adsorption capacities (mg/g). |
start |
Optional named list |
A list with elements Kf, n, R2, method, model (see
fit_langmuir() for details of the structure).
Ce <- c(2, 5, 10, 20, 35, 50) qe <- c(0.8, 1.6, 2.6, 3.6, 4.2, 4.5) fit_freundlich(Ce, qe)Ce <- c(2, 5, 10, 20, 35, 50) qe <- c(0.8, 1.6, 2.6, 3.6, 4.2, 4.5) fit_freundlich(Ce, qe)
Fits the linear form qt = kid * sqrt(t) + C, where kid is the
intraparticle diffusion rate constant (mg/(g min^0.5)) and C is
related to boundary-layer thickness (C = 0 would indicate
intraparticle diffusion is the sole rate-limiting step).
fit_intraparticle(t, qt)fit_intraparticle(t, qt)
t |
Numeric vector of contact times. |
qt |
Numeric vector of adsorption capacity at time |
A list with elements kid, C, R2, model (an lm fit on
qt ~ sqrt(t)).
t <- c(5, 15, 30, 60, 120, 240) qt <- c(1.2, 2.1, 2.9, 3.6, 4.0, 4.2) fit_intraparticle(t, qt)t <- c(5, 15, 30, 60, 120, 240) qt <- c(1.2, 2.1, 2.9, 3.6, 4.0, 4.2) fit_intraparticle(t, qt)
Splits df by group_col (e.g. a sample ID column) and fits the
chosen isotherm model independently within each group, returning one
row of parameters per group. Groups that fail to fit (e.g. too few
points, non-convergence) are skipped with a warning rather than
aborting the whole batch.
fit_isotherm_batch( df, group_col, ce_col, qe_col, model = c("langmuir", "freundlich", "temkin", "dr", "sips"), temperature_K = 298.15 )fit_isotherm_batch( df, group_col, ce_col, qe_col, model = c("langmuir", "freundlich", "temkin", "dr", "sips"), temperature_K = 298.15 )
df |
A data frame containing at least |
group_col |
Name of the grouping column (e.g. sample ID). |
ce_col, qe_col
|
Column names for equilibrium concentration (mg/L) and adsorption capacity (mg/g). |
model |
One of |
temperature_K |
Only used when |
A data frame with one row per successfully fitted group,
columns group_col plus the model's parameter columns and R2
(column names match what the corresponding single-sample fit_*()
function returns).
df <- data.frame( id = rep(c("A", "B"), each = 6), Ce = rep(c(2, 5, 10, 20, 35, 50), 2), qe = c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, # sample A 0.5, 1.0, 1.6, 2.2, 2.7, 3.0) # sample B ) fit_isotherm_batch(df, "id", "Ce", "qe", model = "langmuir")df <- data.frame( id = rep(c("A", "B"), each = 6), Ce = rep(c(2, 5, 10, 20, 35, 50), 2), qe = c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, # sample A 0.5, 1.0, 1.6, 2.2, 2.7, 3.0) # sample B ) fit_isotherm_batch(df, "id", "Ce", "qe", model = "langmuir")
Splits df by group_col (e.g. a sample ID column) and fits the
chosen kinetics model independently within each group, returning one
row of parameters per group. Groups that fail to fit are skipped with
a warning rather than aborting the whole batch.
fit_kinetics_batch( df, group_col, t_col, qt_col, model = c("pfo", "pso", "elovich", "intraparticle") )fit_kinetics_batch( df, group_col, t_col, qt_col, model = c("pfo", "pso", "elovich", "intraparticle") )
df |
A data frame containing at least |
group_col |
Name of the grouping column (e.g. sample ID). |
t_col, qt_col
|
Column names for contact time and adsorption capacity at time t (mg/g). |
model |
One of |
A data frame with one row per successfully fitted group.
df <- data.frame( id = rep(c("A", "B"), each = 6), t = rep(c(5, 15, 30, 60, 120, 240), 2), qt = c(1.2, 2.1, 2.9, 3.6, 4.0, 4.2, 0.8, 1.5, 2.0, 2.5, 2.8, 3.0) ) fit_kinetics_batch(df, "id", "t", "qt", model = "pso")df <- data.frame( id = rep(c("A", "B"), each = 6), t = rep(c(5, 15, 30, 60, 120, 240), 2), qt = c(1.2, 2.1, 2.9, 3.6, 4.0, 4.2, 0.8, 1.5, 2.0, 2.5, 2.8, 3.0) ) fit_kinetics_batch(df, "id", "t", "qt", model = "pso")
Fits qe = (Qmax * KL * Ce) / (1 + KL * Ce) by nonlinear least squares.
If nonlinear fitting fails to converge (common with few points or noisy
data), falls back to the linearised form Ce/qe = Ce/Qmax + 1/(KL*Qmax)
fit by ordinary least squares, and flags the result accordingly.
fit_langmuir(Ce, qe, start = NULL)fit_langmuir(Ce, qe, start = NULL)
Ce |
Numeric vector of equilibrium concentrations (mg/L). |
qe |
Numeric vector of equilibrium adsorption capacities (mg/g). |
start |
Optional named list of starting values |
A list with elements:
Qmax |
Maximum adsorption capacity (mg/g) |
KL |
Langmuir affinity constant (L/mg) |
R2 |
Coefficient of determination |
method |
|
model |
The underlying |
Ce <- c(2, 5, 10, 20, 35, 50) qe <- c(0.8, 1.6, 2.6, 3.6, 4.2, 4.5) fit_langmuir(Ce, qe)Ce <- c(2, 5, 10, 20, 35, 50) qe <- c(0.8, 1.6, 2.6, 3.6, 4.2, 4.5) fit_langmuir(Ce, qe)
Fits qt = qe * (1 - exp(-k1 * t)) by nonlinear least squares.
fit_pfo(t, qt, start = NULL)fit_pfo(t, qt, start = NULL)
t |
Numeric vector of contact times. |
qt |
Numeric vector of adsorption capacity at time |
start |
Optional named list |
A list with elements qe, k1, R2, model (an nls fit).
t <- c(5, 15, 30, 60, 120, 240) qt <- c(1.2, 2.1, 2.9, 3.6, 4.0, 4.2) fit_pfo(t, qt)t <- c(5, 15, 30, 60, 120, 240) qt <- c(1.2, 2.1, 2.9, 3.6, 4.0, 4.2) fit_pfo(t, qt)
qt = (k2 * qe^2 * t) / (1 + k2 * qe * t) by nonlinear least
squares.Fits qt = (k2 * qe^2 * t) / (1 + k2 * qe * t) by nonlinear least
squares.
fit_pso(t, qt, start = NULL)fit_pso(t, qt, start = NULL)
t |
Numeric vector of contact times. |
qt |
Numeric vector of adsorption capacity at time |
start |
Optional named list |
A list with elements qe, k2, R2, model (an nls fit).
t <- c(5, 15, 30, 60, 120, 240) qt <- c(1.2, 2.1, 2.9, 3.6, 4.0, 4.2) fit_pso(t, qt)t <- c(5, 15, 30, 60, 120, 240) qt <- c(1.2, 2.1, 2.9, 3.6, 4.0, 4.2) fit_pso(t, qt)
Fits qe = (Qmax * Ks * Ce^(1/n)) / (1 + Ks * Ce^(1/n)), a three-
parameter hybrid that reduces to Langmuir as n -> 1 and to Freundlich
at low Ce.
fit_sips(Ce, qe, start = NULL)fit_sips(Ce, qe, start = NULL)
Ce |
Numeric vector of equilibrium concentrations (mg/L). |
qe |
Numeric vector of equilibrium adsorption capacities (mg/g). |
start |
Optional named list |
A list with elements Qmax (mg/g), Ks, n, R2, model.
Ce <- c(2, 5, 10, 20, 35, 50, 70) qe <- c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, 4.4) fit_sips(Ce, qe)Ce <- c(2, 5, 10, 20, 35, 50, 70) qe <- c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, 4.4) fit_sips(Ce, qe)
Fits qe = B * log(KT * Ce), where B = RT/b is the Temkin constant
related to the heat of sorption and KT is the equilibrium binding
constant (L/mg).
fit_temkin(Ce, qe, start = NULL)fit_temkin(Ce, qe, start = NULL)
Ce |
Numeric vector of equilibrium concentrations (mg/L). |
qe |
Numeric vector of equilibrium adsorption capacities (mg/g). |
start |
Optional named list |
A list with elements KT, B, R2, model (an nls fit).
Ce <- c(2, 5, 10, 20, 35, 50, 70) qe <- c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, 4.4) fit_temkin(Ce, qe)Ce <- c(2, 5, 10, 20, 35, 50, 70) qe <- c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, 4.4) fit_temkin(Ce, qe)
Computes standard enthalpy (deltaH) and entropy (deltaS) of
adsorption from the linear van't Hoff relationship
log(Kc) = deltaS/R - deltaH/(R * temperature_K), fit by ordinary
least squares across temperatures. Standard Gibbs free energy
deltaG = deltaH - temperature_K * deltaS is then computed at each
supplied temperature.
fit_vant_hoff(temperature_K, Kc)fit_vant_hoff(temperature_K, Kc)
temperature_K |
Numeric vector of absolute temperatures (Kelvin) at which adsorption was measured. |
Kc |
Numeric vector of dimensionless equilibrium distribution
coefficients (e.g. |
A list with elements:
deltaH_kJmol |
Standard enthalpy of adsorption (kJ/mol) |
deltaS_Jmolk |
Standard entropy of adsorption (J/(mol K)) |
table |
A data frame with one row per input temperature:
|
R2 |
Coefficient of determination of the linear van't Hoff fit |
model |
The underlying |
temperature_K <- c(298, 308, 318, 328) Kc <- c(1.8, 2.3, 2.9, 3.4) fit_vant_hoff(temperature_K, Kc)temperature_K <- c(298, 308, 318, 328) Kc <- c(1.8, 2.3, 2.9, 3.4) fit_vant_hoff(temperature_K, Kc)
Returns a data.frame of commonly cited FTIR wavenumber ranges and
their associated functional-group assignments for biochar/carbon
materials, based on standard published band assignments. Intended as a
starting reference for assign_ftir_peaks(); users should adjust
ranges to match their own literature basis where needed.
ftir_band_reference()ftir_band_reference()
A data.frame with columns group, range_low_cm1,
range_high_cm1, assignment.
Computes a simple proxy for functional-group density: the integrated (trapezoidal) absorbance/intensity area within each wavenumber region of a reference table, over a full spectrum (not just picked peaks).
functional_group_density(spectrum, reference = ftir_band_reference())functional_group_density(spectrum, reference = ftir_band_reference())
spectrum |
A |
reference |
A reference table as returned by
|
A data.frame with columns group, range_low_cm1,
range_high_cm1, area (trapezoidal integral of intensity over the
region) and mean_intensity.
Parses sample IDs of the form <PREFIX><TEMPERATURE>-<TIME>
(e.g. "SBC450-30") into their pyrolysis condition components. This
matches naming conventions used for factorial pyrolysis designs where a
feedstock prefix is followed by pyrolysis temperature (degrees C) and
residence time (minutes), separated by a hyphen.
parse_sbc_id(id, prefix_pattern = "[A-Za-z]+")parse_sbc_id(id, prefix_pattern = "[A-Za-z]+")
id |
Character vector of sample identifiers, e.g.
|
prefix_pattern |
Regex describing the non-numeric feedstock prefix.
Defaults to one or more letters ( |
A data.frame with columns id, feedstock, temperature_C,
residence_time_min. Rows that do not match the expected pattern have
NA in the parsed columns and trigger a warning listing the offending
IDs, rather than failing the whole batch.
parse_sbc_id(c("SBC300-15", "SBC450-30", "SBC600-60"))parse_sbc_id(c("SBC300-15", "SBC450-30", "SBC600-60"))
Simple gravimetric helper: 100 * (initial_mass - final_mass) / initial_mass.
Useful for computing moisture % or volatile matter % from raw balance
readings before passing them into proximate_analysis().
pct_mass_change(initial_mass, final_mass)pct_mass_change(initial_mass, final_mass)
initial_mass, final_mass
|
Numeric vectors of mass measurements (any consistent unit, e.g. g). |
Numeric vector of percentage mass loss.
pct_mass_change(initial_mass = 5.000, final_mass = 4.812) # e.g. moisture %pct_mass_change(initial_mass = 5.000, final_mass = 4.812) # e.g. moisture %
Line plot of an FTIR spectrum in conventional display orientation (wavenumber decreasing left to right).
plot_ftir_spectrum(spectrum, main = "FTIR spectrum", ...)plot_ftir_spectrum(spectrum, main = "FTIR spectrum", ...)
spectrum |
A |
main |
Plot title. |
... |
Further arguments passed to |
Invisibly, NULL. Called for its side effect of drawing a plot
on the current graphics device.
Base-graphics scatter of observed (Ce, qe) points with the fitted isotherm curve overlaid.
plot_isotherm(fit, Ce, qe, model = NULL, main = "Adsorption isotherm", ...)plot_isotherm(fit, Ce, qe, model = NULL, main = "Adsorption isotherm", ...)
fit |
A fit object from |
Ce, qe
|
Original data used for the fit. |
model |
Which model |
main |
Plot title. |
... |
Further arguments passed to |
Invisibly, NULL. Called for its side effect of drawing a plot
on the current graphics device.
Base-graphics scatter of observed (t, qt) points with the fitted kinetics curve overlaid.
plot_kinetics( fit, t, qt, model = c("pfo", "pso", "elovich", "intraparticle"), main = "Adsorption kinetics", ... )plot_kinetics( fit, t, qt, model = c("pfo", "pso", "elovich", "intraparticle"), main = "Adsorption kinetics", ... )
fit |
A fit object from |
t, qt
|
Original data used for the fit. |
model |
Which model |
main |
Plot title. |
... |
Further arguments passed to |
Invisibly, NULL. Called for its side effect of drawing a plot
on the current graphics device.
Base-graphics dual-axis plot: the thermogravimetric (TG) weight-loss
curve on the left axis, with the derivative (DTG) curve overlaid as a
dashed line on a right axis. If dtg is not supplied, it is computed
internally via tga_dtg().
plot_tga( tga, dtg = NULL, temperature_col = "temperature_C", weight_col = "weight_pct", main = "TGA/DTG curve", ... )plot_tga( tga, dtg = NULL, temperature_col = "temperature_C", weight_col = "weight_pct", main = "TGA/DTG curve", ... )
tga |
A data frame with columns |
dtg |
Optional pre-computed DTG data frame from |
temperature_col, weight_col
|
Column names in |
main |
Plot title. |
... |
Further arguments passed to the TG |
Invisibly, NULL. Called for its side effect of drawing a plot
on the current graphics device.
set.seed(1) T <- seq(25, 800, by = 5) W <- 100 - 5 * (T > 110) - 40 / (1 + exp(-(T - 340) / 20)) + rnorm(length(T), 0, 0.15) tga <- data.frame(temperature_C = T, weight_pct = pmax(W, 20)) plot_tga(tga)set.seed(1) T <- seq(25, 800, by = 5) W <- 100 - 5 * (T > 110) - 40 / (1 + exp(-(T - 340) / 20)) + rnorm(length(T), 0, 0.15) tga <- data.frame(temperature_C = T, weight_pct = pmax(W, 20)) plot_tga(tga)
Combines independently measured moisture %, volatile matter %, and ash
% (e.g. each computed via pct_mass_change() from its own sub-sample,
following standard gravimetric proximate analysis protocols) into a
full proximate analysis table, computing fixed carbon by difference:
FC% = 100 - moisture% - VM% - ash%.
proximate_analysis(moisture_pct, volatile_matter_pct, ash_pct)proximate_analysis(moisture_pct, volatile_matter_pct, ash_pct)
moisture_pct, volatile_matter_pct, ash_pct
|
Numeric vectors (recycled to a common length) of each percentage, on an as-received mass basis. |
A data frame with columns moisture_pct, VM_pct, ash_pct,
fixed_carbon_pct. A warning is issued (not an error) for any row
where the components sum to more than 100%, since that indicates a
likely measurement or transcription issue rather than a computable
negative fixed-carbon value.
proximate_analysis(moisture_pct = 3.2, volatile_matter_pct = 28.5, ash_pct = 12.1)proximate_analysis(moisture_pct = 3.2, volatile_matter_pct = 28.5, ash_pct = 12.1)
Computes the equilibrium adsorption capacity for a batch adsorption
experiment: qe = (C0 - Ce) * V / m.
qe_batch(C0, Ce, V, m)qe_batch(C0, Ce, V, m)
C0 |
Initial solute concentration (mg/L). |
Ce |
Equilibrium solute concentration (mg/L). |
V |
Solution volume (L). |
m |
Adsorbent (biochar) mass (g). |
Numeric vector of adsorption capacities in mg/g. Recycled
element-wise across C0, Ce, V, m.
qe_batch(C0 = 50, Ce = 12.5, V = 0.05, m = 0.1)qe_batch(C0 = 50, Ce = 12.5, V = 0.05, m = 0.1)
Reads plain-text FTIR exports (e.g. from PerkinElmer ATR-FTIR software)
consisting of two numeric columns (wavenumber and
absorbance/transmittance), possibly preceded by header/metadata lines and
using either comma, tab, or whitespace as a delimiter. Uses base
readLines()/strsplit() rather than a heavier parser so that malformed
or partially-exported files degrade gracefully (bad lines are skipped and
reported) instead of aborting the whole read.
read_ftir_txt( path, skip_pattern = "^\\s*-?[0-9]", col_names = c("wavenumber_cm1", "intensity") )read_ftir_txt( path, skip_pattern = "^\\s*-?[0-9]", col_names = c("wavenumber_cm1", "intensity") )
path |
Path to the text file. |
skip_pattern |
Regex used to identify non-data header lines to skip. Defaults to lines that don't start with an optional minus sign and a digit. |
col_names |
Column names to assign, default |
A data.frame with the two numeric columns, sorted by
wavenumber_cm1 descending (conventional FTIR display order), with a
attr(., "skipped_lines") attribute recording line numbers that could
not be parsed.
Reads plain-text TGA exports (e.g. from TA Instruments, Netzsch, or
Shimadzu software) consisting of two or more numeric columns, typically
temperature, time, and weight percent, possibly preceded by
header/metadata lines and using comma, tab, or whitespace as a
delimiter. Uses the same defensive readLines()/strsplit() approach
as read_ftir_txt() so malformed or partially-exported lines are
skipped and reported rather than aborting the whole read.
read_tga_txt( path, skip_pattern = "^\\s*-?[0-9]", col_names = c("temperature_C", "time_min", "weight_pct") )read_tga_txt( path, skip_pattern = "^\\s*-?[0-9]", col_names = c("temperature_C", "time_min", "weight_pct") )
path |
Path to the text file. |
skip_pattern |
Regex used to identify non-data header lines to skip. Defaults to lines that don't start with an optional minus sign and a digit. |
col_names |
Column names to assign, in the order they appear in
the file. Default |
A data.frame with one column per entry in col_names, sorted
by the first column (temperature) ascending, with a
attr(., "skipped_lines") attribute recording line numbers that
could not be parsed.
tmp <- tempfile(fileext = ".txt") writeLines(c("Temp(C)\tTime(min)\tWeight(%)", "25\t0\t100.0", "50\t2\t99.8", "100\t6\t98.1"), tmp) read_tga_txt(tmp) unlink(tmp)tmp <- tempfile(fileext = ".txt") writeLines(c("Temp(C)\tTime(min)\tWeight(%)", "25\t0\t100.0", "50\t2\t99.8", "100\t6\t98.1"), tmp) read_tga_txt(tmp) unlink(tmp)
Reads plain-text XRD exports (e.g. from Rigaku SmartLab, Cu K-alpha)
consisting of two-theta and intensity columns. Same defensive
readLines()/strsplit() approach as read_ftir_txt().
read_xrd_txt( path, skip_pattern = "^\\s*-?[0-9]", col_names = c("two_theta", "intensity") )read_xrd_txt( path, skip_pattern = "^\\s*-?[0-9]", col_names = c("two_theta", "intensity") )
path |
Path to the text file. |
skip_pattern |
Regex used to identify non-data header lines to skip. Defaults to lines that don't start with an optional minus sign and a digit. |
col_names |
Column names to assign, default |
A data.frame with columns two_theta, intensity, sorted by
two_theta ascending.
Computes percentage removal: 100 * (C0 - Ce) / C0.
removal_efficiency(C0, Ce)removal_efficiency(C0, Ce)
C0 |
Initial solute concentration (mg/L). |
Ce |
Equilibrium solute concentration (mg/L). |
Numeric vector of removal efficiencies (%).
removal_efficiency(C0 = 50, Ce = 12.5)removal_efficiency(C0 = 50, Ce = 12.5)
Wraps a plotting expression in a grDevices::tiff() device at
publication resolution (600 dpi by default) with LZW compression, matching
a common journal figure requirement. The font family defaults to
"serif", which maps to Liberation Serif on most Linux systems (a
metric-compatible substitute for Times New Roman).
save_tiff(expr, filename, width = 6, height = 5, dpi = 600, family = "serif")save_tiff(expr, filename, width = 6, height = 5, dpi = 600, family = "serif")
expr |
An expression (typically a call to one of the |
filename |
Output file path, should end in |
width, height
|
Figure size in inches. |
dpi |
Resolution in dots per inch. Default |
family |
Font family for the device. Default |
Invisibly, the filename.
Ce <- c(2, 5, 10, 20, 35, 50, 70) qe <- c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, 4.4) path <- tempfile(fileext = ".tif") save_tiff(plot_isotherm(fit_langmuir(Ce, qe), Ce, qe), filename = path) file.exists(path) unlink(path)Ce <- c(2, 5, 10, 20, 35, 50, 70) qe <- c(0.8, 1.6, 2.3, 3.0, 3.6, 4.0, 4.4) path <- tempfile(fileext = ".tif") save_tiff(plot_isotherm(fit_langmuir(Ce, qe), Ce, qe), filename = path) file.exists(path) unlink(path)
A small, randomly generated dataset in the shape of a 3x3 factorial pyrolysis design (temperature x residence time), used only to demonstrate and test package functions in examples and tests.
sbc_examplesbc_example
A data frame with 9 rows and 7 columns:
Synthetic sample identifier, e.g. "SBC450-30"
Pyrolysis temperature (degrees C)
Pyrolysis residence time (minutes)
Synthetic pH value
Synthetic electrical conductivity (dS/m)
Synthetic equilibrium solute concentration (mg/L)
Synthetic adsorption capacity (mg/g)
This is synthetic data, not experimental measurements. Values
were generated with a fixed random seed to follow a plausible but
entirely artificial trend, and must not be cited, plotted alongside
real results, or otherwise treated as real biochar characterisation
data. See data-raw/generate_example_data.R for the generating code.
Randomly generated; see data-raw/generate_example_data.R.
Returns a data.frame of commonly cited temperature ranges for the
main decomposition stages seen in biomass/biochar-feedstock TGA
curves, based on standard published ranges. Intended as a starting
reference for assign_dtg_peaks(); users should adjust ranges to
match their own literature basis and feedstock where needed, the same
way ftir_band_reference() is a starting point rather than a fixed
truth.
tga_decomposition_reference()tga_decomposition_reference()
A data.frame with columns stage, range_low_C,
range_high_C, description.
Fits a smoothing spline to the weight-loss curve and differentiates it
analytically (stats::smooth.spline() plus its deriv = 1 predict
method), rather than taking a raw finite difference of noisy
instrument data. Sign convention: dtg is positive during weight
loss (-dW/dT), so decomposition stages show up as positive peaks,
matching how DTG curves are conventionally plotted.
tga_dtg( tga, temperature_col = "temperature_C", weight_col = "weight_pct", spar = NULL )tga_dtg( tga, temperature_col = "temperature_C", weight_col = "weight_pct", spar = NULL )
tga |
A data frame with columns |
temperature_col, weight_col
|
Column names. Defaults
|
spar |
Smoothing parameter passed to |
A data frame with columns temperature_C, weight_pct,
weight_pct_smooth, dtg (percent per degree C), sorted by
temperature_C ascending.
set.seed(1) T <- seq(25, 800, by = 5) W <- 100 - 5 * (T > 110) - 40 / (1 + exp(-(T - 340) / 20)) + rnorm(length(T), 0, 0.15) tga <- data.frame(temperature_C = T, weight_pct = pmax(W, 20)) dtg <- tga_dtg(tga)set.seed(1) T <- seq(25, 800, by = 5) W <- 100 - 5 * (T > 110) - 40 / (1 + exp(-(T - 340) / 20)) + rnorm(length(T), 0, 0.15) tga <- data.frame(temperature_C = T, weight_pct = pmax(W, 20)) dtg <- tga_dtg(tga)
Fits the Kissinger (1957) method for a single decomposition step
observed at several heating rates: ln(beta/Tp^2) = ln(A*R/Ea) - Ea/(R*Tp),
a linear regression of ln(beta/Tp^2) against 1/Tp. Requires the DTG
peak temperature of the same decomposition stage (e.g. matched via
assign_dtg_peaks()) from runs at several different heating rates on
the same sample — not a single TGA curve.
tga_kinetics_kissinger(heating_rate_Kmin, T_peak_C)tga_kinetics_kissinger(heating_rate_Kmin, T_peak_C)
heating_rate_Kmin |
Numeric vector of heating rates (degrees C/min), one per run. |
T_peak_C |
Numeric vector of DTG peak temperatures (degrees C)
for the same decomposition stage, one per run, matched to
|
Because the fit is a plain stats::lm() with an element named
model, fit_confint() works on the result directly, the same as on
any isotherm/kinetics fit in this package.
A list with elements Ea_kJmol (activation energy),
A_min1 (pre-exponential factor, per minute, consistent with
heating_rate_Kmin in degC/min), R2, model (the underlying
lm fit).
# Synthetic: true Ea = 150 kJ/mol, illustrating the expected trend # (peak temperature rises with heating rate). beta <- c(5, 10, 15, 20) Tp_C <- c(320, 335, 344, 351) tga_kinetics_kissinger(beta, Tp_C)# Synthetic: true Ea = 150 kJ/mol, illustrating the expected trend # (peak temperature rises with heating rate). beta <- c(5, 10, 15, 20) Tp_C <- c(320, 335, 344, 351) tga_kinetics_kissinger(beta, Tp_C)
Some instrument exports report raw sample mass (e.g. mg) rather than
weight percent. This rescales a mass column to percent of the first
(initial) reading, the form expected by tga_dtg(), tga_stages(),
and plot_tga().
tga_normalize(tga, mass_col = "weight_mg", out_col = "weight_pct")tga_normalize(tga, mass_col = "weight_mg", out_col = "weight_pct")
tga |
A data frame containing |
mass_col |
Name of the raw mass column. Default |
out_col |
Name to give the new percent column. Default
|
tga with an added/overwritten out_col column.
tga <- data.frame(temperature_C = c(25, 100, 500), weight_mg = c(8.2, 8.1, 5.4)) tga_normalize(tga)tga <- data.frame(temperature_C = c(25, 100, 500), weight_mg = c(8.2, 8.1, 5.4)) tga_normalize(tga)
Reads moisture, volatile matter, and ash percentages directly off a
single TGA weight-loss curve at three temperature breakpoints, then
hands them to proximate_analysis() for fixed-carbon-by-difference
(so the 100%-sum check and warning behaviour is identical either way,
whether the three percentages came from separate gravimetric
sub-samples or from one continuous TGA run).
tga_stages( tga, temperature_col = "temperature_C", weight_col = "weight_pct", moisture_end = 110, vm_end = 650 )tga_stages( tga, temperature_col = "temperature_C", weight_col = "weight_pct", moisture_end = 110, vm_end = 650 )
tga |
A data frame with columns |
temperature_col, weight_col
|
Column names. Defaults
|
moisture_end |
Temperature (degrees C) marking the end of the
drying/moisture-loss stage. Default |
vm_end |
Temperature (degrees C) marking the end of the volatile
matter (devolatilization) stage. Default |
This assumes weight_col is expressed as percent of the original
as-received mass (starting near 100%) and that the run proceeded far
enough for the residual weight at the final temperature to represent
ash (i.e. an oxidative burn-off, per protocols such as ASTM E1131).
If your run stayed under inert atmosphere throughout, the residual at
the final temperature is fixed carbon + ash combined, not ash alone;
in that case treat ash_pct in the output with caution, or supply a
separately measured ash percentage to proximate_analysis() directly
instead.
A data frame as returned by proximate_analysis()
(moisture_pct, VM_pct, ash_pct, fixed_carbon_pct), with two
extra columns moisture_end_C, vm_end_C recording the breakpoints
used.
set.seed(1) T <- seq(25, 800, by = 5) W <- 100 - 5 * (T > 110) - 40 / (1 + exp(-(T - 340) / 20)) + rnorm(length(T), 0, 0.15) tga <- data.frame(temperature_C = T, weight_pct = pmax(W, 20)) tga_stages(tga)set.seed(1) T <- seq(25, 800, by = 5) W <- 100 - 5 * (T > 110) - 40 / (1 + exp(-(T - 340) / 20)) + rnorm(length(T), 0, 0.15) tga <- data.frame(temperature_C = T, weight_pct = pmax(W, 20)) tga_stages(tga)
Splits df by group_col (e.g. a sample ID column) and applies
tga_stages() independently within each group, returning one row per
group. Groups that fail (e.g. too few points, breakpoints outside the
temperature range) are skipped with a warning rather than aborting the
whole batch, mirroring fit_isotherm_batch() and
fit_kinetics_batch().
tga_stages_batch( df, group_col, temperature_col = "temperature_C", weight_col = "weight_pct", moisture_end = 110, vm_end = 650 )tga_stages_batch( df, group_col, temperature_col = "temperature_C", weight_col = "weight_pct", moisture_end = 110, vm_end = 650 )
df |
A long-format data frame containing at least |
group_col |
Name of the grouping column (e.g. sample ID). |
temperature_col, weight_col
|
Column names. Defaults
|
moisture_end |
Temperature (degrees C) marking the end of the
drying/moisture-loss stage. Default |
vm_end |
Temperature (degrees C) marking the end of the volatile
matter (devolatilization) stage. Default |
A data frame with one row per successfully processed group,
group_col plus the columns returned by tga_stages().
set.seed(1) mk <- function(id, vm_center) { T <- seq(25, 800, by = 10) W <- 100 - 5 * (T > 110) - 40 / (1 + exp(-(T - vm_center) / 20)) + rnorm(length(T), 0, 0.15) data.frame(id = id, temperature_C = T, weight_pct = pmax(W, 20)) } df <- rbind(mk("SBC300-30", 300), mk("SBC450-30", 340), mk("SBC600-30", 380)) tga_stages_batch(df, "id")set.seed(1) mk <- function(id, vm_center) { T <- seq(25, 800, by = 10) W <- 100 - 5 * (T > 110) - 40 / (1 + exp(-(T - vm_center) / 20)) + rnorm(length(T), 0, 0.15) data.frame(id = id, temperature_C = T, weight_pct = pmax(W, 20)) } df <- rbind(mk("SBC300-30", 300), mk("SBC450-30", 340), mk("SBC600-30", 380)) tga_stages_batch(df, "id")
Converts elemental mass percentages (from a CHNS/CHNSO analyser) to atomic ratios (H/C, O/C, and optionally N/C), the standard inputs for a Van Krevelen diagram. Uses standard atomic weights (C 12.011, H 1.008, O 15.999, N 14.007).
ultimate_ratios(C_pct, H_pct, O_pct, N_pct = NULL)ultimate_ratios(C_pct, H_pct, O_pct, N_pct = NULL)
C_pct, H_pct, O_pct
|
Numeric vectors of carbon, hydrogen, and oxygen mass percentages. |
N_pct |
Optional numeric vector of nitrogen mass percentage; if
supplied, |
Note: O_pct is very often obtained by difference
(100 - C - H - N - S - ash) rather than measured directly; supply
whichever value your own workflow produced.
A data frame with columns HC_ratio, OC_ratio, and (if
N_pct supplied) NC_ratio.
ultimate_ratios(C_pct = 65.2, H_pct = 2.8, O_pct = 18.4, N_pct = 1.1)ultimate_ratios(C_pct = 65.2, H_pct = 2.8, O_pct = 18.4, N_pct = 1.1)
Computes a simple crystallinity index as the ratio of crystalline peak
area to total (crystalline + amorphous) area, following the common
peak-area deconvolution approach used for carbon materials:
CI = A_crystalline / (A_crystalline + A_amorphous).
xrd_crystallinity_index(area_crystalline, area_amorphous)xrd_crystallinity_index(area_crystalline, area_amorphous)
area_crystalline |
Numeric: integrated area of crystalline peak(s). |
area_amorphous |
Numeric: integrated area of the amorphous "hump". |
This function does not perform peak fitting/deconvolution itself; it expects the user to have already separated crystalline vs. amorphous contributions (e.g. via baseline subtraction and peak integration in their preferred XRD software), and simply computes the ratio from the resulting areas. This keeps the function honest about not fabricating peak-fitting results it cannot verify.
Numeric crystallinity index between 0 and 1.
xrd_crystallinity_index(area_crystalline = 1200, area_amorphous = 800)xrd_crystallinity_index(area_crystalline = 1200, area_amorphous = 800)
Fits a sum of Gaussian peaks to an XRD pattern at user-supplied
approximate peak centers (2theta), by simultaneous nonlinear least
squares. From the fitted peak areas and the total pattern area, also
computes a crystallinity index (see xrd_crystallinity_index()) by
treating everything not accounted for by the fitted crystalline peaks
as the amorphous contribution.
xrd_deconvolve(two_theta, intensity, peak_centers, start_width = 1)xrd_deconvolve(two_theta, intensity, peak_centers, start_width = 1)
two_theta, intensity
|
Numeric vectors: the diffraction angle
(degrees 2theta) and intensity of the pattern. Should already be
baseline-corrected (e.g. via an approach analogous to
|
peak_centers |
Numeric vector of approximate 2theta positions for the crystalline peaks to fit. |
start_width |
Starting value (degrees 2theta) for each peak's
Gaussian width. Default |
This is a practical, simplified tool intended for a handful of
well-separated peaks on an already baseline-corrected pattern — it is
not a substitute for Rietveld refinement or dedicated diffraction
software, and may fail to converge for patterns with many strongly
overlapping peaks (in which case, try narrowing peak_centers to the
most prominent peaks, or refine starting values manually).
A list with elements:
peaks |
A data frame with one row per peak: |
crystallinity_index |
|
R2 |
Coefficient of determination of the overall multi-Gaussian fit |
model |
The underlying |
set.seed(1) 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_deconvolve(two_theta, intensity, peak_centers = c(22, 29))set.seed(1) 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_deconvolve(two_theta, intensity, peak_centers = c(22, 29))