The package currently supports three varieties of stomatal conductance models (Ball et al. 1987; Leuning 1995; Medlyn et al. 2011).

library(dplyr)
library(photosynthesis)

round_to_nearest = function(x, values) {
  sapply(x, function(y, values) {
    values[which.min(abs(y - values))]
  }, values = values)
}

# Read in data
dat = system.file("extdata", "A_Ci_Q_data_1.csv", package = "photosynthesis") |>
  read.csv() |>
  # Convert RH to a proportion
  mutate(
    RH = RHcham / 100,
    PPFD = round_to_nearest(Qin, c(25, 50, 75, 100, 125, 375, 750, 1500))
  ) |>
  rename(A_net = A, C_air = Ca, g_sw = gsw, VPD = VPDleaf)

# Fit the Ball-Berry stomatal conductance models
fit = fit_gs_model(
  data = filter(dat, PPFD == 750),
  model = c("BallBerry")
  )

# Look at BallBerry model summary:
summary(fit[["BallBerry"]][["Model"]])

## 
## Call:
## lm(formula = g_sw ~ gs_mod_ballberry(A_net = A_net, C_air = C_air, 
##     RH = RH), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.05185 -0.03045  0.01009  0.02768  0.04035 
## 
## Coefficients:
##                                                          Estimate Std. Error
## (Intercept)                                             3.659e-01  1.303e-02
## gs_mod_ballberry(A_net = A_net, C_air = C_air, RH = RH) 3.148e-06  1.423e-06
##                                                         t value Pr(>|t|)    
## (Intercept)                                              28.076 7.63e-11 ***
## gs_mod_ballberry(A_net = A_net, C_air = C_air, RH = RH)   2.213   0.0513 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03363 on 10 degrees of freedom
## Multiple R-squared:  0.3287, Adjusted R-squared:  0.2615 
## F-statistic: 4.896 on 1 and 10 DF,  p-value: 0.05133

# Look at BallBerry parameters
fit[["BallBerry"]][["Parameters"]]

##          g0           g1
## 1 0.3659298 3.148369e-06

# Look at BallBerry plot
# fit[["BallBerry"]][["Graph"]]

# Fit many g_sw models
fits = fit_many(dat, funct = fit_gs_model, group = "PPFD", progress = FALSE)

## Error in nlsModel(formula, mf, start, wts) : 
##   singular gradient matrix at initial parameter estimates

# Look at the Medlyn_partial outputs at 750 PAR
# Model summary
summary(fits[["750"]][["Medlyn_partial"]][["Model"]])

## 
## Formula: g_sw ~ gs_mod_opti(A_net = A_net, C_air = C_air, VPD = VPD, g0, 
##     g1)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## g0  0.38778    0.03317  11.692 3.73e-07 ***
## g1 -1.09754    0.83022  -1.322    0.216    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04375 on 10 degrees of freedom
## 
## Number of iterations to convergence: 2 
## Achieved convergence tolerance: 1.49e-08

# Model parameters
fits[["750"]][["Medlyn_partial"]][["Parameters"]]

##          g0        g1
## 1 0.3877773 -1.097544

# Graph
# fits[["750"]][["Medlyn_partial"]][["Graph"]]

# Compile parameter outputs for BallBerry model
# Note that it's the first element for each PAR value
# First compile list of BallBerry fits
bbmods = compile_data(data = fits, output_type = "list", list_element = 1)

# Now compile the parameters (2nd element) into a dataframe
bbpars = compile_data(data = bbmods, output_type = "dataframe", list_element = 2)

#Compile graphs
graphs = compile_data(data = bbmods, output_type = "list", list_element = 3)

# Look at 3rd graph
# graphs[[3]]