R Under development (unstable) (2025-07-25 r88451 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(chapensk)
Loading required package: Bessel
> 
> # graphics parameters
> width <- 10
> height <- 10
> mar <- c(5.1, 4.6, 0.2, 0.6)
> mai <- c(1.02, 0.82, 0.04, 0.04)
> family <- "Helvetica"
> pointsize <- 12
> lwd <- 1
> 
> Diffusion <- function(T, D0, b) {
+   return <- D0 * (T / pkg.env$T0)^b
+ }
> 
> helium <- Gas("helium")
> nitrogen <- Gas("nitrogen")
> argon <- Gas("argon")
> methane <- Gas("methane")
> ethane <- Gas("ethane")
> propane <- Gas("propane")
> butane <- Gas("n-butane")
> fluoromethane <- Gas("fluoromethane")
> difluoromethane <- Gas("difluoromethane")
> trifluoromethane <- Gas("trifluoromethane")
> 
> if (interactive()) {
+   if (exists("outputDevice")) {
+     if (outputDevice == "pdf") pdf(file = "man/figures/binary_diffusion.pdf", width = width, height = height, family = family, pointsize = pointsize)
+     if (outputDevice == "svg") svg(file = "man/figures/binary_diffusion.svg", width = width, height = height, family = family, pointsize = pointsize)
+   }
+   split.screen(figs = c(2, 2))
+ }
> 
> nitrogen_in_helium <- subset(binary_diffusion, (gas == "nitrogen" & bath_gas == "helium"))
> argon_in_helium <- subset(binary_diffusion, (gas == "argon" & bath_gas == "helium"))
> 
> if (interactive()) {
+   screen(1)
+ }
> par(mai = mai, mar = mar)
> 
> plot(nitrogen_in_helium$T, nitrogen_in_helium$D,
+   log = "xy",
+   xlab = expression(paste(italic(T), " / [K]")),
+   ylab = expression(paste(italic(D), "/ [cm2/s]")),
+   pch = 0, col = "blue"
+ )
> Dcalc <- nitrogen$binary_diffusion(T = nitrogen_in_helium$T, bathGas = helium)
> nitrogen_in_helium$Dcalc <- 1E4 * Dcalc
> lines(nitrogen_in_helium$T, 1E4 * Dcalc, col = "blue")
> nls_n2_he <- nls(D ~ Diffusion(T, D0, b), data = nitrogen_in_helium, weights = 1 / U_D, start = list(D0 = 0.62, b = 1.7), trace = TRUE)
3.363581    (1.61e+01): par = (0.62 1.7)
0.01371689  (2.96e-01): par = (0.6048772 1.666074)
0.01261376  (3.89e-04): par = (0.6049878 1.664479)
0.01261376  (6.81e-07): par = (0.604989 1.664475)
> print(summary(nls_n2_he))

Formula: D ~ Diffusion(T, D0, b)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
D0 0.604989   0.001488   406.6   <2e-16 ***
b  1.664475   0.003849   432.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03971 on 8 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 6.805e-07

> points(argon_in_helium$T, argon_in_helium$D, pch = 1, col = "red")
> Dcalc <- argon$binary_diffusion(T = argon_in_helium$T, bathGas = helium)
> argon_in_helium$Dcalc <- 1E4 * Dcalc
> lines(argon_in_helium$T, 1E4 * Dcalc, col = "red")
> nls_ar_he <- nls(D ~ Diffusion(T, D0, b), data = argon_in_helium, weights = 1 / U_D, start = list(D0 = 0.62, b = 1.7), trace = TRUE)
0.1761939   (4.29e+00): par = (0.62 1.7)
0.009025495 (3.42e-02): par = (0.6298177 1.665571)
0.009014922 (5.05e-06): par = (0.6299077 1.66545)
> print(summary(nls_ar_he))

Formula: D ~ Diffusion(T, D0, b)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
D0 0.629908   0.001231   511.6   <2e-16 ***
b  1.665450   0.003093   538.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03357 on 8 degrees of freedom

Number of iterations to convergence: 2 
Achieved convergence tolerance: 5.054e-06

> legend("bottomright",
+   legend = c("nitrogen", "argon"),
+   col = c("blue", "red"),
+   pch = c(0, 1)
+ )
> if (interactive()) {
+   legend("top", bty = "n", legend = "(a)")
+   close.screen()
+ }
> 
> methane_in_helium <- subset(binary_diffusion, (gas == "methane" & bath_gas == "helium"))
> ethane_in_helium <- subset(binary_diffusion, (gas == "ethane" & bath_gas == "helium"))
> propane_in_helium <- subset(binary_diffusion, (gas == "propane" & bath_gas == "helium"))
> butane_in_helium <- subset(binary_diffusion, (gas == "butane" & bath_gas == "helium"))
> 
> if (interactive()) {
+   screen(2)
+ }
> par(mai = mai, mar = mar)
> 
> plot(methane_in_helium$T, methane_in_helium$D,
+   log = "xy",
+   xlab = expression(paste(italic(T), " / [K]")),
+   ylab = expression(paste(italic(D), "/ [cm2/s]")),
+   pch = 0, col = "blue",
+   ylim = c(0.3, 3)
+ )
> Dcalc <- methane$binary_diffusion(T = nitrogen_in_helium$T, bathGas = helium)
> lines(methane_in_helium$T, 1E4 * Dcalc, col = "blue")
> methane_in_helium$Dcalc <- 1E4 * Dcalc
> nls_ch4_he <- nls(D ~ Diffusion(T, D0, b), data = methane_in_helium, weights = 1 / U_D, start = list(D0 = 0.62, b = 1.7), trace = TRUE)
14.44066    (3.77e+01): par = (0.62 1.7)
0.01271827  (5.18e-01): par = (0.5745618 1.677166)
0.01002813  (7.95e-04): par = (0.5746216 1.67497)
0.01002813  (1.54e-06): par = (0.5746235 1.674962)
> print(summary(nls_ch4_he))

Formula: D ~ Diffusion(T, D0, b)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
D0 0.574623   0.001273   451.4   <2e-16 ***
b  1.674962   0.003500   478.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03541 on 8 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 1.544e-06

> points(ethane_in_helium$T, ethane_in_helium$D, pch = 1, col = "red")
> Dcalc <- ethane$binary_diffusion(T = nitrogen_in_helium$T, bathGas = helium)
> lines(ethane_in_helium$T, 1E4 * Dcalc, col = "red")
> ethane_in_helium$Dcalc <- 1E4 * Dcalc
> nls_c2h6_he <- nls(D ~ Diffusion(T, D0, b), data = ethane_in_helium, start = list(D0 = 0.62, b = 1.7), trace = TRUE)
4.932062    (1.17e+02): par = (0.62 1.7)
0.0008203367 (1.13e+00): par = (0.4205432 1.688201)
0.0003604450 (3.41e-03): par = (0.4205818 1.682375)
0.0003604408 (8.84e-06): par = (0.4205928 1.682329)
> print(summary(nls_c2h6_he))

Formula: D ~ Diffusion(T, D0, b)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
D0 0.420593   0.002712   155.1 3.34e-15 ***
b  1.682329   0.007748   217.1 2.26e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.006712 on 8 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 8.839e-06

> points(propane_in_helium$T, propane_in_helium$D, pch = 2, col = "dark green")
> Dcalc <- propane$binary_diffusion(T = propane_in_helium$T, bathGas = helium)
> lines(propane_in_helium$T, 1E4 * Dcalc, col = "dark green")
> propane_in_helium$Dcalc <- 1E4 * Dcalc
> nls_c3h8_he <- nls(D ~ Diffusion(T, D0, b), data = propane_in_helium, start = list(D0 = 0.62, b = 1.7), trace = TRUE)
9.650276    (1.35e+02): par = (0.62 1.7)
0.002134078 (1.74e+00): par = (0.3406273 1.684231)
0.0005308258 (1.10e-02): par = (0.340691 1.670784)
0.0005307615 (2.66e-05): par = (0.3407205 1.670601)
0.0005307615 (1.41e-07): par = (0.3407207 1.6706)
> print(summary(nls_c3h8_he))

Formula: D ~ Diffusion(T, D0, b)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
D0 0.340721   0.003308   103.0 8.81e-14 ***
b  1.670600   0.011676   143.1 6.37e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.008145 on 8 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 1.411e-07

> points(butane_in_helium$T, butane_in_helium$D, pch = 3, col = "orange")
> Dcalc <- butane$binary_diffusion(T = butane_in_helium$T, bathGas = helium)
> lines(butane_in_helium$T, 1E4 * Dcalc, col = "orange")
> butane_in_helium$Dcalc <- 1E4 * Dcalc
> nls_c4h10_he <- nls(D ~ Diffusion(T, D0, b), data = butane_in_helium, start = list(D0 = 0.62, b = 1.7), trace = TRUE)
13.31892    (1.77e+02): par = (0.62 1.7)
0.004852906 (3.23e+00): par = (0.293759 1.677561)
0.0004258618 (3.68e-02): par = (0.2938575 1.651525)
0.0004252864 (1.90e-05): par = (0.2939035 1.651046)
0.0004252864 (2.08e-08): par = (0.2939036 1.651046)
> print(summary(nls_c4h10_he))

Formula: D ~ Diffusion(T, D0, b)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
D0 0.293904   0.002986   98.41 1.27e-13 ***
b  1.651046   0.012238  134.91 1.02e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.007291 on 8 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 2.08e-08

> legend("bottomright",
+   legend = c("methane", "ethane", "propane", "butane"),
+   col = c("blue", "red", "dark green", "orange"),
+   pch = c(0, 1, 2, 3)
+ )
> if (interactive()) {
+   legend("top", bty = "n", legend = "(b)")
+   close.screen()
+ }
> 
> if (interactive()) {
+   screen(3)
+ }
> par(mai = mai, mar = mar)
> methane_in_nitrogen <- subset(binary_diffusion, (gas == "methane" & bath_gas == "nitrogen"))
> ethane_in_nitrogen <- subset(binary_diffusion, (gas == "ethane" & bath_gas == "nitrogen"))
> propane_in_nitrogen <- subset(binary_diffusion, (gas == "propane" & bath_gas == "nitrogen"))
> butane_in_nitrogen <- subset(binary_diffusion, (gas == "butane" & bath_gas == "nitrogen"))
> plot(methane_in_nitrogen$T, methane_in_nitrogen$D,
+   log = "xy",
+   xlab = expression(paste(italic(T), " / [K]")),
+   ylab = expression(paste(italic(D), "/ [cm2/s]")),
+   pch = 0, col = "blue",
+   ylim = c(0.1, 1)
+ )
> Dcalc <- methane$binary_diffusion(T = methane_in_nitrogen$T, bathGas = nitrogen)
> methane_in_nitrogen$Dcalc <- 1E4 * Dcalc
> lines(methane_in_nitrogen$T, 1E4 * Dcalc, col = "blue")
> nls_ch4_n2 <- nls(D ~ Diffusion(T, D0, b), data = methane_in_nitrogen, weights = 1 / U_D, start = list(D0 = 0.62, b = 1.7), trace = TRUE)
1139.267    (3.65e+02): par = (0.62 1.7)
0.04987031  (2.19e+00): par = (0.2005089 1.712088)
0.008642811 (1.87e-02): par = (0.2005155 1.737083)
0.008639780 (3.77e-05): par = (0.2005248 1.736796)
0.008639780 (1.05e-07): par = (0.2005247 1.736797)
> print(summary(nls_ch4_n2))

Formula: D ~ Diffusion(T, D0, b)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
D0 0.2005247  0.0009176   218.5 6.19e-13 ***
b  1.7367970  0.0091826   189.1 1.47e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03795 on 6 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 1.047e-07

> points(ethane_in_nitrogen$T, ethane_in_nitrogen$D, pch = 1, col = "red")
> Dcalc <- ethane$binary_diffusion(T = ethane_in_nitrogen$T, bathGas = nitrogen)
> lines(ethane_in_nitrogen$T, 1E4 * Dcalc, col = "red")
> ethane_in_nitrogen$Dcalc <- 1E4 * Dcalc
> nls_c2h6_n2 <- nls(D ~ Diffusion(T, D0, b), data = ethane_in_nitrogen, weights = 1 / U_D, start = list(D0 = 0.62, b = 1.7), trace = TRUE)
9376.359    (5.66e+02): par = (0.62 1.7)
0.03179131  (2.95e-01): par = (0.1361981 1.700665)
0.02924257  (6.68e-04): par = (0.1361985 1.703023)
0.02924255  (2.92e-06): par = (0.1361998 1.703008)
> print(summary(nls_c2h6_n2))

Formula: D ~ Diffusion(T, D0, b)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
D0 0.1361998  0.0009141   149.0 4.60e-15 ***
b  1.7030078  0.0085603   198.9 4.56e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.06046 on 8 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.923e-06

> points(propane_in_nitrogen$T, propane_in_nitrogen$D, pch = 2, col = "dark green")
> Dcalc <- propane$binary_diffusion(T = propane_in_nitrogen$T, bathGas = nitrogen)
> lines(propane_in_nitrogen$T, 1E4 * Dcalc, col = "dark green")
> propane_in_nitrogen$Dcalc <- 1E4 * Dcalc
> nls_c3h8_n2 <- nls(D ~ Diffusion(T, D0, b), data = propane_in_nitrogen, weights = 1 / U_D, start = list(D0 = 0.62, b = 1.7), trace = TRUE)
5717.460    (5.03e+02): par = (0.62 1.7)
0.05137649  (1.13e+00): par = (0.1055947 1.70353)
0.02265381  (1.05e-02): par = (0.1055985 1.720616)
0.02265132  (7.63e-05): par = (0.1056118 1.720299)
0.02265132  (6.97e-07): par = (0.1056116 1.720302)
> print(summary(nls_c3h8_n2))

Formula: D ~ Diffusion(T, D0, b)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
D0 0.1056116  0.0008852   119.3 2.72e-14 ***
b  1.7203018  0.0131642   130.7 1.31e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05321 on 8 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 6.966e-07

> points(butane_in_nitrogen$T, butane_in_nitrogen$D, pch = 3, col = "orange")
> Dcalc <- butane$binary_diffusion(T = butane_in_nitrogen$T, bathGas = nitrogen)
> lines(butane_in_nitrogen$T, 1E4 * Dcalc, col = "orange")
> butane_in_nitrogen$Dcalc <- 1E4 * Dcalc
> nls_c4h10_n2 <- nls(D ~ Diffusion(T, D0, b), data = butane_in_nitrogen, weights = 1 / U_D, start = list(D0 = 0.62, b = 1.7), trace = TRUE)
6403.839    (9.03e+02): par = (0.62 1.7)
0.02645333  (1.54e+00): par = (0.09003597 1.702571)
0.007874967 (1.06e-02): par = (0.0900379 1.717643)
0.007874089 (3.28e-05): par = (0.09004512 1.717435)
0.007874089 (1.56e-07): par = (0.09004507 1.717436)
> print(summary(nls_c4h10_n2))

Formula: D ~ Diffusion(T, D0, b)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
D0 0.0900451  0.0006509   138.3 2.72e-13 ***
b  1.7174356  0.0109012   157.5 1.10e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03354 on 7 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 1.559e-07

> legend("topleft",
+   legend = c("methane", "ethane", "propane", "butane"),
+   col = c("blue", "red", "dark green", "orange"),
+   pch = c(0, 1, 2, 3)
+ )
> if (interactive()) {
+   legend("top", bty = "n", legend = "(c)")
+   close.screen()
+ }
> 
> if (interactive()) {
+   screen(4)
+ }
> par(mai = mai, mar = mar)
> fluoromethane_in_nitrogen <- subset(binary_diffusion, (gas == "fluoromethane" & bath_gas == "nitrogen"))
> difluoromethane_in_nitrogen <- subset(binary_diffusion, (gas == "difluoromethane" & bath_gas == "nitrogen"))
> trifluoromethane_in_nitrogen <- subset(binary_diffusion, (gas == "trifluoromethane" & bath_gas == "nitrogen"))
> plot(fluoromethane_in_nitrogen$T, fluoromethane_in_nitrogen$D,
+   log = "xy",
+   xlab = expression(paste(italic(T), " / [K]")),
+   ylab = expression(paste(italic(D), "/ [cm2/s]")),
+   pch = 0, col = "blue",
+   ylim = c(0.1, 0.6)
+ )
> fluoromethane_chisq <- function(p) {
+   fluoromethane$sigma <- p[1]
+   fluoromethane$zeta <- 0
+   fluoromethane$epsk <- p[2]
+   Dcalc <- 1E4 * fluoromethane$binary_diffusion(T = fluoromethane_in_nitrogen$T, bathGas = nitrogen)
+   return <- sum(((fluoromethane_in_nitrogen$D - Dcalc) / fluoromethane_in_nitrogen$U_D)^2)
+ }
> optim_fluoromethane <- optim(par = c(4, 180), fluoromethane_chisq, method = "L-BFGS-B", control = list(trace = 1, parscale = c(1, 100)))
iter   10 value 0.255421
final  value 0.248795 
converged
> print(optim_fluoromethane$par[1])
[1] 3.515299
> print(optim_fluoromethane$par[2])
[1] 173.8092
> fluoromethane$sigma <- optim_fluoromethane$par[1]
> fluoromethane$zeta <- 0
> fluoromethane$epsk <- optim_fluoromethane$par[2]
> Dcalc <- fluoromethane$binary_diffusion(T = fluoromethane_in_nitrogen$T, bathGas = nitrogen)
> lines(fluoromethane_in_nitrogen$T, 1E4 * Dcalc, col = "blue")
> fluoromethane_in_nitrogen$Dcalc <- 1E4 * Dcalc
> 
> points(difluoromethane_in_nitrogen$T, difluoromethane_in_nitrogen$D, pch = 1, col = "red")
> difluoromethane_chisq <- function(p) {
+   difluoromethane$sigma <- p[1]
+   difluoromethane$zeta <- 0
+   difluoromethane$epsk <- p[2]
+   Dcalc <- 1E4 * difluoromethane$binary_diffusion(T = difluoromethane_in_nitrogen$T, bathGas = nitrogen)
+   return <- sum(((difluoromethane_in_nitrogen$D - Dcalc) / difluoromethane_in_nitrogen$U_D)^2)
+ }
> optim_difluoromethane <- optim(par = c(4, 180), difluoromethane_chisq, method = "L-BFGS-B", control = list(trace = 1, parscale = c(1, 100)))
iter   10 value 1.120562
final  value 1.120319 
converged
> print(optim_difluoromethane$par[1])
[1] 3.931231
> print(optim_difluoromethane$par[2])
[1] 153.4328
> difluoromethane$sigma <- optim_difluoromethane$par[1]
> difluoromethane$zeta <- 0
> difluoromethane$epsk <- optim_difluoromethane$par[2]
> Dcalc <- difluoromethane$binary_diffusion(T = difluoromethane_in_nitrogen$T, bathGas = nitrogen)
> lines(difluoromethane_in_nitrogen$T, 1E4 * Dcalc, col = "red")
> difluoromethane_in_nitrogen$Dcalc <- 1E4 * Dcalc
> 
> points(trifluoromethane_in_nitrogen$T, trifluoromethane_in_nitrogen$D, pch = 2, col = "dark green")
> trifluoromethane_chisq <- function(p) {
+   trifluoromethane$sigma <- p[1]
+   trifluoromethane$zeta <- 0
+   trifluoromethane$epsk <- p[2]
+   Dcalc <- 1E4 * trifluoromethane$binary_diffusion(T = trifluoromethane_in_nitrogen$T, bathGas = nitrogen)
+   return <- sum(((trifluoromethane_in_nitrogen$D - Dcalc) / trifluoromethane_in_nitrogen$U_D)^2)
+ }
> optim_trifluoromethane <- optim(par = c(4, 180), trifluoromethane_chisq, method = "L-BFGS-B", control = list(trace = 1, parscale = c(1, 100)))
iter   10 value 3.403804
iter   20 value 2.378332
final  value 2.378332 
converged
> print(optim_trifluoromethane$par[1])
[1] 4.502888
> print(optim_trifluoromethane$par[2])
[1] 63.83145
> trifluoromethane$sigma <- optim_trifluoromethane$par[1]
> trifluoromethane$zeta <- 0
> trifluoromethane$epsk <- optim_trifluoromethane$par[2]
> Dcalc <- trifluoromethane$binary_diffusion(T = trifluoromethane_in_nitrogen$T, bathGas = nitrogen)
> lines(trifluoromethane_in_nitrogen$T, 1E4 * Dcalc, col = "dark green")
> trifluoromethane_in_nitrogen$Dcalc <- 1E4 * Dcalc
> legend("bottomright",
+   legend = c("fluoromethane", "difluoromethane", "trifluoromethane"),
+   col = c("blue", "red", "dark green"),
+   pch = c(0, 1, 2)
+ )
> 
> if (interactive()) {
+   legend("top", bty = "n", legend = "(d)")
+   close.screen()
+ }
> 
> dev.off()
null device 
          1 
> 
> proc.time()
   user  system elapsed 
   1.40    0.09    1.48