1. Fitting emission lines

All four use forageLIBS: 368 spectra of dried forage samples, with laboratory concentrations of 12 elements.

Each row of spectra is one spectrum, the mean of 8 laser shots, and the column names are the wavelengths (nm):

Identifying the lines

libs_lines() lists the strongest lines of given species in a wavelength range, from the NIST Atomic Spectra Database, and plot_lines() marks them on a spectrum. line_finder() opens an interactive app that does both. These functions need an internet connection:

candidates <- libs_lines(c("K I", "Ca I", "Fe I", "Mn I"), wavelength = c(400, 410))
plot_lines(colMeans(spectra), candidates)
line_finder(spectra)

Correcting the wavelength axis

The measured lines are shifted from their tabulated wavelengths. wavelength_calibration() locates known lines in the mean spectrum and fits a correction of the axis. The instrument joins several detectors (the wavelengths step back at 766 nm), and each detector has its own calibration: we correct the first one, from 199 to 766 nm, which holds the lines used below.

first_detector <- spectra[seq_len(which(diff(wl) < 0)[1])]

reference <- c(`C I` = 247.856, `Mg I` = 285.213, `Ca II` = 317.933, `Ca II` = 370.603,
               `Ca II` = 373.690, `K I` = 404.414, `Ca I` = 430.253, `Ca I` = 445.478,
               `Mg I` = 516.732, `Mg I` = 517.268, `Mg I` = 518.360, `Na I` = 588.995,
               `Na I` = 589.592, `Ca I` = 616.217, `Ca I` = 643.907)
wavecal <- wavelength_calibration(first_detector, reference)
wavecal$segments
#> # A tibble: 1 × 6
#>   segment  from    to degree n_lines   rmse
#>     <int> <dbl> <dbl>  <int>   <int>  <dbl>
#> 1       1  199.  766.      1      15 0.0136

plot_wavelength_calibration(wavecal)

The offset grows from 0.05 nm in the UV to 0.12 nm in the red, about one channel, and the residuals of the linear correction are close to 0.01 nm. apply_calibration() relabels the columns with the corrected wavelengths. We also remove the continuum with baseline_arpls(), and average the spectra for a clean example:

calibrated <- apply_calibration(first_detector, wavecal)
corrected <- baseline_arpls(calibrated, lambda = 1e5, max.iter = 20)$correction
wl_corrected <- as.numeric(names(corrected))
mean_spectrum <- colMeans(corrected)

# The channels between two wavelengths, as a one-row table
window <- function(from, to, x = mean_spectrum) {
  keep <- wl_corrected > from & wl_corrected < to
  if (is.null(dim(x))) as_tibble(as.list(x[keep])) else x[keep]
}

Line profiles

An emission line of a laser-induced plasma is broadened by several mechanisms: the instrument and the Doppler effect give a Gaussian profile, the Stark effect a Lorentzian one, and both together a Voigt profile. specProc provides the four profiles, parametrized by their full widths at half maximum:

Function	Profile	Widths
`gaussian_profile()`	Gaussian	`wG`
`lorentzian_profile()`	Lorentzian	`wL`
`voigt_profile()`	exact Voigt (Faddeeva function, C++)	`wG`, `wL`
`pseudo_voigt_profile()`	Thompson–Cox–Hastings pseudo-Voigt	`wG`, `wL`

Fitting a single line

peak_fit() fits \(y = y_0 + A \cdot f(x; x_c, w)\) to each spectrum (one per row, wavelengths as column names) by nonlinear least squares, with starting values estimated from the data. The hydrogen H\(\alpha\) line at 656.28 nm comes from the water of the samples. Its width is dominated by the Stark effect, and the choice of profile matters:

halpha <- window(653.5, 659.5)
fits <- c("gaussian", "lorentzian", "pseudo_voigt", "voigt") |>
  set_names() |>
  map(\(p) peak_fit(halpha, profile = p))

tibble(
  profile = names(fits),
  AIC = map_dbl(fits, \(f) AIC(f$fit[[1]])),
  area = map_dbl(fits, \(f) coef(f$fit[[1]])[["A"]])
) |>
  mutate(across(c(AIC, area), round))
#> # A tibble: 4 × 3
#>   profile        AIC  area
#>   <chr>        <dbl> <dbl>
#> 1 gaussian      1081 14650
#> 2 lorentzian    1002 25884
#> 3 pseudo_voigt  1004 25886
#> 4 voigt         1004 25884

The Gaussian misses the wings of the line: its AIC is 80 units higher, and its area 40% smaller. The Lorentzian, pseudo-Voigt and Voigt profiles fit equally well. tidied holds the parameters, with their standard errors:

fits$voigt$tidied[[1]]
#> # A tibble: 5 × 5
#>   term      estimate  std.error statistic    p.value
#>   <chr>        <dbl>      <dbl>     <dbl>      <dbl>
#> 1 y0      260.        141.        1.84e+0 7.05 e-  2
#> 2 xc      657.          0.00843   7.79e+4 1.60 e-264
#> 3 wG        0.000841  115.        7.30e-6 1.000e+  0
#> 4 wL        1.69        0.138     1.23e+1 9.35 e- 19
#> 5 A     25884.       1376.        1.88e+1 3.27 e- 28

The Gaussian width of the Voigt fit is close to zero, with a standard error far larger than its value: the line is purely Lorentzian, and the Voigt profile has one parameter too many. plot_fit() shows the Lorentzian fit and its residuals:

plot_fit(fits$lorentzian, title = "H-alpha 656.28 nm, Lorentzian profile")

The width of H\(\alpha\) gives the electron density of the plasma (see vignette("plasma-diagnostics")); for a Voigt fit, voigt_fwhm() combines the two widths into the total width.

Overlapping lines

The potassium doublet at 404.41 and 404.72 nm is weak and its lines overlap. multipeak_fit() fits them together, on a common baseline; the parameters of line \(i\) carry the suffix _i:

doublet <- window(404.1, 405.0)
k_voigt <- multipeak_fit(doublet, peaks = c(404.414, 404.721), profiles = "voigt")
k_voigt$tidied[[1]]
#> # A tibble: 9 × 5
#>   term     estimate std.error statistic   p.value
#>   <chr>       <dbl>     <dbl>     <dbl>     <dbl>
#> 1 y0    -137.        428.      -3.19e-1 0.803    
#> 2 xc_1   404.          0.0230   1.76e+4 0.0000362
#> 3 wG_1     0.0712      0.224    3.17e-1 0.804    
#> 4 wL_1     0.0657      0.223    2.94e-1 0.818    
#> 5 A_1    285.        194.       1.47e+0 0.380    
#> 6 xc_2   405.          0.0155   2.61e+4 0.0000244
#> 7 wG_2     0.000837  107.       7.82e-6 1.000    
#> 8 wL_2     0.145       0.503    2.88e-1 0.821    
#> 9 A_2    343.        682.       5.03e-1 0.703

With about ten channels for 9 parameters, the Gaussian and Lorentzian widths of each line cannot be separated: their standard errors are larger than the estimates, and the areas are poorly determined. A Gaussian profile, with one width per line, is enough for lines only a few channels wide:

k_gauss <- multipeak_fit(doublet, peaks = c(404.414, 404.721), profiles = "gaussian")
k_gauss$tidied[[1]]
#> # A tibble: 7 × 5
#>   term  estimate std.error statistic  p.value
#>   <chr>    <dbl>     <dbl>     <dbl>    <dbl>
#> 1 y0     -14.4    97.4        -0.148 8.92e- 1
#> 2 xc_1   404.      0.00743 54465.    1.36e-14
#> 3 wG_1     0.119   0.0126      9.43  2.53e- 3
#> 4 A_1    232.     33.5         6.94  6.12e- 3
#> 5 xc_2   405.      0.0112  36253.    4.63e-14
#> 6 wG_2     0.183   0.0330      5.54  1.16e- 2
#> 7 A_2    234.     48.0         4.86  1.66e- 2
c(voigt = AIC(k_voigt$fit[[1]]), gaussian = AIC(k_gauss$fit[[1]]))
#>    voigt gaussian 
#> 135.1920 133.2817

plot_fit(k_gauss, title = "K I 404.41 and 404.72 nm, Gaussian profiles")

The areas are now determined, and the AIC is lower: prefer the simplest profile that describes the data.

Fitting every spectrum

With an id column, peak_fit() and multipeak_fit() fit each row, so the same line is measured in every spectrum. A fit that fails (here on a few spectra where the doublet is lost in the noise) gives NULL with a warning, and the others continue:

doublets <- bind_cols(
  select(forageLIBS, Measurement),
  window(404.1, 405.0, x = corrected)
)
k_fits <- multipeak_fit(doublets, peaks = c(404.414, 404.721), profiles = "gaussian",
                        id = "Measurement")

k_areas <- k_fits |>
  mutate(area = map_dbl(tidied, \(t) if (is.null(t)) NA else t$estimate[t$term == "A_1"])) |>
  select(Measurement, area) |>
  left_join(select(forageLIBS, Measurement, K), by = "Measurement")
sum(is.na(k_areas$area))
#> [1] 8

ggplot(k_areas, aes(K, area)) +
  geom_point(alpha = 0.5) +
  labs(x = "K (%)", y = "Area of K I 404.41 nm") +
  theme_bw()

The area follows the potassium content, with much scatter: the intensity of a line also depends on the ablated mass and the plasma, which the next vignettes correct by normalization and calibration.

Integrating instead of fitting

line_intensities() measures many lines at once: it finds the peak of each line near its tabulated wavelength (search) and integrates it over a window (half_width), or fits a Voigt profile with method = "voigt". It needs no starting values and is fast, but it does not separate overlapping lines:

lines <- c(`K I 404.41` = 404.414, `Mg I 518.36` = 518.360, `Fe II 259.94` = 259.940)
line_intensities(corrected[1:2, ], lines, half_width = 0.1)
#> # A tibble: 6 × 10
#>   spectrum line       wavelength peak_wavelength    shift height intensity   snr
#>      <int> <chr>           <dbl>           <dbl>    <dbl>  <dbl>     <dbl> <dbl>
#> 1        1 K I 404.41       404.            404.  0.00970  2908.     321.   66.9
#> 2        1 Mg I 518.…       518.            518.  0.0247   8692.    1006.  200. 
#> 3        1 Fe II 259…       260.            260. -0.00470   614.      65.3  14.1
#> 4        2 K I 404.41       404.            404.  0.00970  2092.     209.   47.0
#> 5        2 Mg I 518.…       518.            518.  0.0247   5547.     639.  125. 
#> 6        2 Fe II 259…       260.            260. -0.00470   596.      64.4  13.4
#> # ℹ 2 more variables: detected <lgl>, saturated <lgl>

The recipe step step_line_intensities() does the same within a preprocessing pipeline (vignette("preprocessing")).

Summary

Correct the wavelength axis first with wavelength_calibration(), so that lines can be found and identified at their tabulated wavelengths.
Choose the profile from the data: compare the AIC and the residuals. Strongly Stark-broadened lines such as H\(\alpha\) need a Lorentzian or Voigt profile; lines a few channels wide are often better fitted by a Gaussian.
Fit overlapping lines together with multipeak_fit(), and check that each parameter is determined (its standard error well below its value).
Keep the same profile and window for every spectrum of a study, and fit them all with id.