diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index f9bc4effd..cbea44ffa 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -48,7 +48,9 @@ jobs:
 
       - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          extra-packages: any::rcmdcheck
+          extra-packages: 
+            any::rcmdcheck,
+            vdiffr@1.0.5
           needs: check
 
       - uses: r-lib/actions/check-r-package@v2
diff --git a/DESCRIPTION b/DESCRIPTION
index 1fc01775c..4cccac57c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -53,14 +53,14 @@ Imports:
     viridisLite
 Suggests:
     alphashape3d,
-    digest,
     imager,
     knitr,
     mapproj,
     rgl,
     rmarkdown,
     testthat(>= 2.99.0),
-    vdiffr
+    vdiffr,
+    digest
 VignetteBuilder: 
     knitr
 Config/Needs/website:
diff --git a/NAMESPACE b/NAMESPACE
index 2c336cea6..287ea6d0f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -9,6 +9,7 @@ S3method(plot,rimg)
 S3method(plot,rspec)
 S3method(plot,sensmod)
 S3method(points,colspace)
+S3method(stitch,rspec)
 S3method(subset,colspace)
 S3method(subset,rspec)
 S3method(subset,vismodel)
@@ -64,6 +65,7 @@ export(sensdata)
 export(sensmodel)
 export(simulate_spec)
 export(spec2rgb)
+export(stitch)
 export(tcspace)
 export(tcsplot)
 export(tcspoints)
diff --git a/NEWS.md b/NEWS.md
index 88f520005..ec9ed34e0 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,11 +2,14 @@
 
 ## NEW FEATURES AND SIGNIFICANT CHANGES
 
+- Added a new function `stitch.rspec()`, which facilitates the row-wise merging of spectra. For example, when combining spectra whose wavelength ranges do not overlap, or overlap only slightly (as is common when recording separate UV-VIS and NIR spectra). See `?stitch.rspec` for examples and information, as well as the handbook.
 - Added a new function `simulate_spec()`, which allows for the flexible simulation of naturalistic spectra (inc. reflectance, irradiance, radiance, absorbance). See `?simulate_spec` for examples and information, and the handbook for further discussion.
 - `plot.rspec()` now accepts a logical `labels` argument (and `labels.cex`), to control whether text labels identifying each spectrum should be added to the outer plot margins. This was previously only available, and was required, for 'stacked' plot types, but is now optional for both 'overlay' (the default) and 'stacked' spectral plots.
+- the `wlmin` and `wlmax` arguments in `summary.rspec()` are being deprecated in favour of a single `lim` argument, for consistency across functions.
 
 ## MINOR FEATURES AND BUG FIXES
 
+- `summary.rspec()` has been rewritten for efficiency, and now only calculates the required variables when `subset` is used. As a result, the function is also slightly slower (0.5 x) when calculating the full set of variables, but much faster (10 x) when calculating a subset.
 - Removed the start-up message.
 - Removed the previously-deprecated `margins` argument from various colourspace plots.
 - Replace `rgl.triangles` with `triangles3d()` internally to avoid a deprecation issue.
diff --git a/R/stitch.rspec.R b/R/stitch.rspec.R
new file mode 100644
index 000000000..ee26598dd
--- /dev/null
+++ b/R/stitch.rspec.R
@@ -0,0 +1,155 @@
+#' Stitch together two rspec objects
+#'
+#' Stitch (row-wise merge) two `rspec` objects of differing wavelength ranges into
+#' a single `rspec` object.
+#'
+#' @param rspec1,rspec2 (required) `rspec` objects of differing wavelength ranges
+#' to stitch together by row.
+#' @param overlap_method the method for modifying reflectance values if regions
+#' of the spectra overlap in their wavelength range. Defaults to `mean`.
+#' @param interp logical argument specifying whether reflectance values should be
+#' interpolated between the two sets of spectra if their wavelength ranges
+#' do not overlap. Defaults to `TRUE`.
+#'
+#' @export
+#'
+#' @examples
+#'
+#' # Simulate a UV-VIS and NIR reflectance spectrum whose wavelength regions
+#' # slightly overlap then stitch them together, with the overlapping
+#' # regions being averaged.
+#'
+#' # Simulate specs
+#' reflect1 <- simulate_spec(wl_peak = 550, xlim = c(300, 700))
+#' reflect2 <- simulate_spec(wl_inflect = 1100, xlim = c(650, 1200))
+#'
+#' # Ensure the names of the spectra match
+#' names(reflect1) <- names(reflect2) <- c('wl', 'sample_1')
+#'
+#' # Stitch the spectra together by their wavelength column
+#' full_spec <- stitch(reflect1, reflect2)
+#'
+#' # Plot the resulting spectrum
+#' plot(full_spec)
+#'
+#' # Simulate another set of UV-VIS and NIR spectra. Note two additional complexities,
+#' # both of which are handled without issue. First, the wavelength ranges are
+#' # non-overlapping (with a 100 nm gap). We'll keep the default interp = TRUE argument
+#' # to allow the missing reflectance region to be interpolated. Second, the names of
+#' # the spectra match, but are in a different order in the two rspec objects. This isn't
+#' # an issue, as the function can match up the spectra by name irrespective of their
+#' # ordering
+#'
+#' # Simulate UV-VIS and NIR spectra
+#' reflect_vis <- merge(simulate_spec(wl_peak = 550, xlim = c(300, 700)),
+#'                      simulate_spec(wl_peak = 550, xlim = c(300, 700)))
+#' reflect_nir <- merge(simulate_spec(wl_inflect = 1000, xlim = c(800, 1250)),
+#'                      simulate_spec(wl_inflect = 1100, xlim = c(800, 1250)))
+#'
+#' # Ensure the names of the spectra exist in each, albeit in a different order
+#' names(reflect_vis) <- c('wl', 'sample_1', 'sample_2')
+#' names(reflect_nir) <- c('wl', 'sample_2', 'sample_1')
+#'
+#' # Stitch together by their wavelength column, with missing regions being
+#' # interpolated
+#' reflect_vis_nir <- stitch(reflect_vis, reflect_nir)
+#'
+#' # Plot the resulting spectrum
+#' plot(reflect_vis_nir)
+#'
+#' @author Thomas White \email{thomas.white026@@gmail.com}
+#' @author Hugo Gruson \email{hugo.gruson+R@@normalesup.org}
+#'
+#' @seealso [as.rspec()], [merge.rspec()]
+
+stitch <- function(rspec1, rspec2, overlap_method, interp) {
+  UseMethod("stitch")
+}
+
+#' @rdname stitch
+#'
+#' @export
+stitch.rspec <- function(rspec1, rspec2,
+                                 overlap_method = c("mean", "minimum", "maximum"),
+                                 interp = TRUE) {
+
+  # Class check
+  if (!inherits(rspec1, "rspec") || !inherits(rspec2, "rspec")) {
+    stop("Both inputs must be of class 'rspec'")
+  }
+
+  # Validate overlap_method
+  overlap_method <- match.arg(overlap_method)
+
+  # Check that at least one spectrum has a matching name in both objects
+  common_cols <- intersect(names(rspec1), names(rspec2))
+  if (length(common_cols) <= 1) {
+    stop("At least one spectrum in both rspec objects must have a matching name")
+  }
+
+  # Warn if only subset is present across both rspec objects
+  if (length(common_cols) != ncol(rspec1) || length(common_cols) != ncol(rspec2)) {
+    warning("Not all spectra are present in both objects. Stitching only the common samples.")
+  }
+
+  # Identify unique spectra in both objects
+  unique_rspec1 <- setdiff(names(rspec1), common_cols)
+  unique_rspec2 <- setdiff(names(rspec2), common_cols)
+
+  # Create NA-filled columns in each rspec for unique spectra in the other
+  rspec1[, unique_rspec2] <- NA
+  rspec2[, unique_rspec1] <- NA
+
+  # Reorder columns of rspec2 to match rspec1
+  rspec2 <- rspec2[, names(rspec1)]
+
+  # Merge by wl
+  res <- rbind(rspec1, rspec2)
+
+  # Handle overlapping wl values
+  if (anyDuplicated(res$wl) > 0) {
+    overlap_wl <- unique(res$wl[duplicated(res$wl)])
+
+    for (wl in overlap_wl) {
+      idx <- which(res$wl == wl)
+
+      # Replace with a switch statement
+      switch(overlap_method,
+             mean = res[idx[1], -1] <- colMeans(as.matrix(res[idx, -1]), na.rm = TRUE),
+             minimum = res[idx[1], -1] <- apply(as.matrix(res[idx, -1]), 2, min, na.rm = TRUE),
+             maximum = res[idx[1], -1] <- apply(as.matrix(res[idx, -1]), 2, max, na.rm = TRUE)
+      )
+
+      # Remove extra rows
+      if (length(idx) > 1) {
+        res <- res[-idx[-1], ]
+      }
+    }
+  }
+
+  # Interpolate missing values
+  if (interp) {
+    full_wl_range <- min(res$wl):max(res$wl)
+    missing_wl <- setdiff(full_wl_range, res$wl)
+
+    if (length(missing_wl) > 0) {
+      new_rows <- data.frame(wl = missing_wl, matrix(rep(NA, length(names(res)) - 1), ncol = length(names(res)) - 1))
+      names(new_rows)[-1] <- names(res)[-1]
+
+      # Interpolate only common spectra
+      for (col in common_cols[-1]) {
+        new_values <- approx(res$wl, res[, col], xout = missing_wl)$y
+        new_rows[, col] <- new_values
+      }
+      res <- rbind(res, new_rows)
+    }
+  }
+
+  # Sort stitched spec by wl
+  res <- res[order(res$wl), ]
+
+  # Classes
+  class(res) <- c("rspec", "data.frame")
+
+  res
+}
diff --git a/R/summary.rspec.R b/R/summary.rspec.R
index f3f584a64..8abe287fa 100644
--- a/R/summary.rspec.R
+++ b/R/summary.rspec.R
@@ -10,8 +10,9 @@
 #' of the complete output (composed of B2, S8 and H1; the variables described in
 #' Andersson and Prager 2006) are returned. Finally, a user-specified string of variable
 #' names can be used in order to filter and show only those variables.
-#' @param wlmin,wlmax minimum and maximum used to define the range of wavelengths used in
-#' calculations (default is to use entire range in the `rspec` object).
+#' @param lim The range of wavelengths used in calculations. The default is to use
+#' the entire range in the `rspec` object (typically equivalent to `lim = c(300, 700)`).
+#' @param wlmin,wlmax Deprecated. Use the `lim` argument instead.
 #' @param ... class consistency (ignored)
 #'
 #' @return A data frame containing either 23 or 5 (`subset = TRUE`) variables described
@@ -91,6 +92,7 @@
 #'
 #' H5 (Hue): Wavelength at bmax. Sensitive to noise and may be variable if there is
 #' more than one maxima and minima. REF 5
+#'
 #' @note If minimum wavelength is over 400, UV chroma is not computed.
 #' @note Variables which compute bmax and bmaxneg should be used with caution, for they
 #' rely on smoothed curves to remove noise, which would otherwise result in spurious
@@ -103,11 +105,21 @@
 #' @export
 #'
 #' @examples
+#' # Load data
 #' data(sicalis)
+#'
+#' # Calculate and display all spectral summary variables
 #' summary(sicalis)
+#'
+#' # Calculate only subset of B2, S8 and H1 as per Andersson (1999)
 #' summary(sicalis, subset = TRUE)
+#'
+#' # Calculate user-specified subset of B1 and H4
 #' summary(sicalis, subset = c("B1", "H4"))
-#' @author Pierre-Paul Bitton \email{bittonp@@windsor.ca}, Rafael Maia \email{rm72@@zips.uakron.edu}
+#'
+#' @author Thomas E. White \email{thomas.white026@@gmail.com}
+#' @author Pierre-Paul Bitton \email{bittonp@@windsor.ca}
+#' @author Rafael Maia \email{rm72@@zips.uakron.edu}
 #'
 #' @references Montgomerie R. 2006. Analyzing colors. In Hill, G.E, and McGraw, K.J., eds.
 #' Bird Coloration. Volume 1 Mechanisms and measurements. Harvard University Press, Cambridge, Massachusetts.
@@ -158,230 +170,415 @@
 #' 13- Smiseth, P., J. Ornborg, S. Andersson, and T. Amundsen. 2001. Is male plumage reflectance
 #' correlated with paternal care in bluethroats? Behavioural Ecology 12:164-170.
 #'
-summary.rspec <- function(object, subset = FALSE, wlmin = NULL, wlmax = NULL, ...) {
-  chkDots(...)
+summary.rspec <- function(object, subset = FALSE, lim = NULL, wlmin = NULL, wlmax = NULL, ...) {
+    chkDots(...)
+
+    wl <- isolate_wl(object, keep = "wl")
+
+    # The control-flow here is for cases when users specify only of the deprecated
+    # wlmin or wlmax arguments. Can be removed once deprecated in next release.
+    if (!all(missing("wlmin"), missing("wlmin"))) {
+      warning("The 'wlmin' and 'wlmax' arguments are deprecated as of v2.9.0,
+      and will be removed in future releases. Use the lim() argument instead.")
+      if (!missing("wlmin") && missing("wlmax"))
+        lim  <- c(wlmin, max(wl))
+      if (missing("wlmin") && !missing("wlmax"))
+        lim  <- c(min(wl), wlmax)
+      if (all(missing("wlmin"), missing("wlmin")))
+        lim  <- c(wlmin, wlmax)
+    }
+
+    lambdamin <- max(lim[1], min(wl))
+    lambdamax <- min(lim[2], max(wl))
+
+    # wl-range checks
+    if (!missing(lim)) {
+      if (lambdamin > lim[1]) {
+        stop("Minimum specified wavelength is smaller than the range of spectral data. Check the lim argument.")
+      }
+      if (lambdamax < lim[2]) {
+        stop("Maximum specified wavelength is larger than the range of spectral data. Check the lim argument.")
+      }
+    }
 
-  wl <- isolate_wl(object, keep = "wl")
+    # Restrict to range of wlmin:wlmax
+    object <- object[which(wl == lambdamin):which(wl == lambdamax), ]
 
-  # Set WL min & max
-  lambdamin <- max(wlmin, min(wl))
-  lambdamax <- min(wlmax, max(wl))
+    wl <- isolate_wl(object, keep = "wl")
+    object <- isolate_wl(object, keep = "spec")
 
-  if (!is.null(wlmin) && lambdamin > wlmin) {
-    stop("wlmin is smaller than the range of spectral data")
-  }
-  if (!is.null(wlmax) && lambdamax < wlmax) {
-    stop("wlmax is larger than the range of spectral data")
-  }
+    # Establish variables and handle any subsetting
+    variable_names <- c("B1", "B2", "B3", "S1", "S2", "S3", "S4", "S5", "S6",
+                        "S7", "S8", "S9", "S10", "H1", "H2", "H3", "H4", "H5")
 
-  # Restrict to range of wlmin:wlmax
-  object <- object[which(wl == lambdamin):which(wl == lambdamax), ]
+    if (is.logical(subset)) {
+      if (subset) {
+        variable_names <- c("B2", "S8", "H1")
+      }
+    } else {
+      if (all(subset %in% variable_names)) {
+        variable_names <- subset
+      } else {
+        stop("Names in ", dQuote("subset"), " do not match color variable names")
+      }
+    }
+
+    # Preserve original variable names
+    variable_names_orig <- variable_names
 
-  wl <- isolate_wl(object, keep = "wl")
-  object <- isolate_wl(object, keep = "spec")
+    # Replace S1 with sub-variables
+    if ("S1" %in% variable_names) {
+      # Find index of 'S1'
+      s1_index <- which(variable_names == "S1")
 
-  output.mat <- matrix(nrow = ncol(object), ncol = 23)
+      # Remove 'S1'
+      variable_names <- variable_names[-s1_index]
 
-  # Three measures of brightness
-  B1 <- colSums(object)
+      # Full names
+      new_var_names <- c("S1U", "S1V", "S1B", "S1G", "S1Y", "S1R")
 
-  B2 <- colMeans(object)
+      # Append full names at the right place
+      variable_names <- append(variable_names, new_var_names, after = s1_index - 1)
+    }
 
-  B3 <- vapply(object, max, numeric(1))
+    # Define final summary output df
+    color.var <- data.frame(matrix(ncol = length(variable_names),
+                                   nrow = ncol(object)))
+    names(color.var) <- variable_names
+
+    # Dictionary to store computed variables
+    computed_vars <- list()
+
+    # Define a dependency structure among variables.
+    # This allows us to speed up computation by computing only the necessary
+    # variables (and their dependencies) when subset is not FALSE.
+    dependencies <- list(
+      B1 = list(fun = calc_B1, deps = NULL),
+      B2 = list(fun = calc_B2, deps = NULL),
+      B3 = list(fun = calc_B3, deps = NULL),
+      S1 = list(fun = calc_S1, deps = "B1"),
+      S2 = list(fun = calc_S2, deps = c("B3", "Rmin")),
+      S3 = list(fun = calc_S3, deps = c("H1", "B1")),
+      S4 = list(fun = calc_S4, deps = NULL),
+      S5 = list(fun = calc_S5, deps = NULL),
+      S6 = list(fun = calc_S6, deps = c("B3", "Rmin")),
+      S7 = list(fun = calc_S7, deps = c("B1", "Rmid")),
+      S8 = list(fun = calc_S8, deps = c("S6", "B2")),
+      S9 = list(fun = calc_S9, deps = NULL),
+      S10 = list(fun = calc_S10, deps = c("S8", "S4")),
+      H1 = list(fun = calc_H1, deps = NULL),
+      H2 = list(fun = calc_H2, deps = NULL),
+      H3 = list(fun = calc_H3, deps = c("B3", "Rmin", "Rmid")),
+      H4 = list(fun = calc_H4, deps = NULL),
+      H5 = list(fun = calc_H5, deps = NULL),
+      Rmin = list(fun = calc_Rmin, deps = NULL),
+      Rmid = list(fun = calc_Rmid, deps = NULL)
+    )
+
+    # Function to calculate a variable, which ensures each variable
+    # is calculated only once
+    compute_color_var <- function(var_name) {
+
+       # If the variable is already computed, return it
+      if (!is.null(computed_vars[[var_name]])) {
+        return(computed_vars[[var_name]])
+      }
+
+      # Get the function and its dependencies
+      fun <- dependencies[[var_name]]$fun
+      deps <- dependencies[[var_name]]$deps
+
+      # If the function has dependencies, calculate them first
+      if (!is.null(deps)) {
+        for (dep in deps) {
+          computed_vars[[dep]] <- compute_color_var(dep)
+        }
+      }
+
+      # Calculate and store variable
+      computed_vars[[var_name]] <- fun(object, wl, computed_vars)
+
+      # Special case handling for S1, which returns a list
+      if (var_name == "S1") {
+        for (sub_var_name in names(computed_vars[[var_name]])) {
+          color.var[, sub_var_name] <- computed_vars[[var_name]][[sub_var_name]]
+        }
+      } else {
+        color.var[, var_name] <- computed_vars[[var_name]]
+      }
+
+      # Return the computed variable
+      return(computed_vars[[var_name]])
+    }
 
-  Rmin <- vapply(object, min, numeric(1))
+    # Calculate all specified variables
+    for (var_name in variable_names_orig) {
+      if (var_name == "S1") {
+        color.var[, c("S1U", "S1V", "S1B", "S1G", "S1Y", "S1R")] <- compute_color_var("S1")
+      } else {
+        color.var[, var_name] <- compute_color_var(var_name)
+      }
+    }
 
-  Rmid <- (B3 + Rmin) / 2
+    row.names(color.var) <- names(object)
 
-  # Chromas
+    if ("Rmid" %in% names(color.var))
+      color.var <- color.var[, !names(color.var) %in% "Rmid"]
+    if ("Rmin" %in% names(color.var))
+      color.var <- color.var[, !names(color.var) %in% "Rmin"]
 
-  # Red
-  if (lambdamin <= 605 && lambdamax >= 700) {
-    Redchromamat <- object[wl >= 605 & wl <= 700, , drop = FALSE] # red 605-700nm inclusive
-    Redchroma <- colSums(Redchromamat) / B1 # S1 red
-    output.mat[, 9] <- Redchroma
-  } else {
-    warning("cannot calculate red chroma; wavelength range not between 605 and 700 nm", call. = FALSE)
+    # The double-conversion here is just so the attributes are in the same order
+    # as in the original formulation, for CI test consistency
+    as.data.frame(as.matrix(color.var))
   }
 
-  # Yellow
-  if (lambdamin <= 550 && lambdamax >= 625) {
-    Yellowchromamat <- object[wl >= 550 & wl <= 625, , drop = FALSE] # yellow 550-625nm
-    Yellowchroma <- colSums(Yellowchromamat) / B1 # S1 yellow
-    output.mat[, 8] <- Yellowchroma
-  } else {
-    warning("cannot calculate yellow chroma; wavelength range not between 550 and 625 nm", call. = FALSE)
-  }
 
-  # Green
-  if (lambdamin <= 510 && lambdamax >= 605) {
-    Greenchromamat <- object[wl >= 510 & wl <= 605, , drop = FALSE] # green 510-605nm inlusive
-    Greenchroma <- colSums(Greenchromamat) / B1 # S1 green
-    output.mat[, 7] <- Greenchroma
-  } else {
-    warning("cannot calculate green chroma; wavelength range not between 510 and 605 nm", call. = FALSE)
-  }
+## ----- Common intermediaries ----- ##
 
-  # Blue
-  if (lambdamin <= 400 && lambdamax >= 510) {
-    Bluechromamat <- object[wl >= 400 & wl <= 510, , drop = FALSE] # blue 400-510nm inclusive
-    Bluechroma <- colSums(Bluechromamat) / B1 # S1 blue
-    output.mat[, 6] <- Bluechroma
-  } else {
-    warning("cannot calculate blue chroma; wavelength range not between 400 and 510 nm", call. = FALSE)
-  }
+calc_Rmin <- function(object, wl, lambdamin, lambdamax) {
+  vapply(object, min, numeric(1))
+}
+
+calc_Rmid <- function(object, wl, lambdamin, lambdamax) {
+  (calc_B3(object) + calc_Rmin(object)) / 2
+}
+
+## ----- Brightness ----- ##
+
+calc_B1 <- function(object, wl, computed_vars, ...) {
+  colSums(object)
+}
 
-  # UV
-  if (lambdamin <= 400 && lambdamax >= 400) {
-    UVchromamat <- object[wl >= lambdamin & wl <= 400, , drop = FALSE]
-    UVchroma <- colSums(UVchromamat) / B1 # S1 UV
-    output.mat [, 4] <- UVchroma
-  } else {
-    warning("cannot calculate UV chroma; wavelength range not below 400 nm", call. = FALSE)
+calc_B2 <- function(object, wl, computed_vars, ...) {
+  colMeans(object)
+}
+
+calc_B3 <- function(object, wl, computed_vars, ...) {
+  vapply(object, max, numeric(1))
+}
+
+## ----- Saturation ----- ##
+
+calc_S1 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  B1 <- computed_vars[["B1"]]
+  lambdamin <- min(wl)
+  lambdamax <- max(wl)
+
+  # Define the chroma ranges
+  chroma_ranges <- list(
+    S1U = list(min = lambdamin, max = 400,
+      message = "cannot calculate UV chroma; wavelength range not below 400 nm"
+    ),
+    S1V = list(min = lambdamin, max = 415,
+      message = "cannot calculate violet chroma; wavelength below 415 nm"
+    ),
+    S1B = list(min = 400, max = 510,
+      message = "cannot calculate blue chroma; wavelength range not between 400 and 510 nm"
+    ),
+    S1G = list(min = 510, max = 605,
+      message = "cannot calculate green chroma; wavelength range not between 510 and 605 nm"
+    ),
+    S1Y = list(min = 550, max = 625,
+      message = "cannot calculate yellow chroma; wavelength range not between 550 and 625 nm"
+    ),
+    S1R = list(min = 605, max = 700,
+      message = "cannot calculate red chroma; wavelength range not between 605 and 700 nm"
+    )
+  )
+
+  output.mat <- matrix(ncol = length(chroma_ranges), nrow = ncol(object))
+
+  for (i in seq_along(chroma_ranges)) {
+    if (lambdamin <= chroma_ranges[[i]]$min && lambdamax >= chroma_ranges[[i]]$max) {
+      chromamat <- object[wl >= chroma_ranges[[i]]$min & wl <= chroma_ranges[[i]]$max, , drop = FALSE]
+      chroma <- colSums(chromamat) / B1
+      output.mat[, i] <- chroma
+    } else {
+      warning(chroma_ranges[[i]]$message, call. = FALSE)
+    }
   }
 
   if (lambdamin > 300 && lambdamin < 400) {
     warning("Minimum wavelength is ", lambdamin, "; UV-related variables may not be meaningful", call. = FALSE)
   }
 
-  # Violet
-  if (lambdamin <= 415 && lambdamax >= 415) {
-    Vchromamat <- object[wl >= lambdamin & wl <= 415, , drop = FALSE]
-    Vchroma <- colSums(Vchromamat) / B1 # S1 Violet
-    output.mat[, 5] <- Vchroma
-  } else {
-    warning("cannot calculate violet chroma; wavelength below 415 nm", call. = FALSE)
-  }
+  return(
+    list(
+      S1U = output.mat[, 1],
+      S1V = output.mat[, 2],
+      S1B = output.mat[, 3],
+      S1G = output.mat[, 4],
+      S1Y = output.mat[, 5],
+      S1R = output.mat[, 6]
+    )
+  )
+}
 
-  # lambda Rmax hue
-  H1 <- wl[max.col(t(object), ties.method = "first")]
+calc_S2 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  B3 <- computed_vars[["B3"]]
+  Rmin <- computed_vars[["Rmin"]]
 
-  # Segment-based variables
+  # Calc
+  B3 / Rmin
+}
 
-  segmts <- trunc(quantile(lambdamin:lambdamax, names = FALSE))
+calc_S3 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  H1 <- computed_vars[["H1"]]
+  B1 <- computed_vars[["B1"]]
 
+  # Calc
+  S3 <- vapply(seq_len(ncol(object)), function(col) {
+    spec <- object[, col]
+    H1_spec <- H1[col]
+    sum(spec[wl >= (H1_spec - 50) & wl <= (H1_spec + 50)])
+  }, numeric(1))
+  S3 / B1
+}
+
+calc_S4 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  diffsmooth <- apply(object, 2, diff)
+  incr <- apply(diffsmooth, 2, min) > 0
+
+  # Calc
+  bmaxneg <- abs(apply(diffsmooth, 2, min))
+  bmaxneg[incr] <- NA
+  bmaxneg
+}
+
+calc_S5 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  lambdamin <- min(wl)
+  lambdamax <- max(wl)
+
+  # Calc
+  segmts <- trunc(quantile(lambdamin:lambdamax, names = FALSE))
   Q1 <- wl >= segmts[1] & wl <= segmts[2]
   Q2 <- wl >= segmts[2] & wl <= segmts[3]
   Q3 <- wl >= segmts[3] & wl <= segmts[4]
   Q4 <- wl >= segmts[4] & wl <= segmts[5]
-
   S5R <- colSums(object[Q4, , drop = FALSE])
   S5Y <- colSums(object[Q3, , drop = FALSE])
   S5G <- colSums(object[Q2, , drop = FALSE])
   S5B <- colSums(object[Q1, , drop = FALSE])
+  sqrt((S5R - S5G) ^ 2 + (S5Y - S5B) ^ 2)
+}
 
-  S5 <- sqrt((S5R - S5G)^2 + (S5Y - S5B)^2)
+calc_S6 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  B3 <- computed_vars[["B3"]]
+  Rmin <- computed_vars[["Rmin"]]
 
-  H4 <- atan2(S5Y - S5B, S5R - S5G)
+  # Calc
+  B3 - Rmin
+}
 
-  # Carotenoid chroma
+calc_S7 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  B1 <- computed_vars[["B1"]]
+  Rmid <- computed_vars[["Rmid"]]
 
-  R450 <- as.numeric(object[which(wl == 450), , drop = FALSE])
-  R700 <- as.numeric(object[which(wl == 700), , drop = FALSE])
-  Carotchroma <- (R700 - R450) / R700
-
-  # H3
+  # Calc
   index_Rmid <- vapply(seq_len(ncol(object)), function(x) {
     which.min(abs(object[, x] - Rmid[x]))
   }, numeric(1))
-  H3 <- wl[index_Rmid]
-
-  # S7
-
   S7 <- vapply(seq_len(ncol(object)), function(col) {
     spec <- object[, col]
     index_Rmid_spec <- index_Rmid[col]
     spec_low <- spec[seq_len(index_Rmid_spec)]
     spec_high <- spec[index_Rmid_spec:length(spec)]
-
     return(sum(spec_low) - sum(spec_high))
   }, numeric(1))
+  S7 / B1
+}
 
-  S7 <- S7 / B1
+calc_S8 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  S6 <- computed_vars[["S6"]]
+  B2 <- computed_vars[["B2"]]
 
+  # Calc
+  S6 / B2
+}
 
-  # S3
-  S3 <- vapply(seq_len(ncol(object)), function(col) {
-    spec <- object[, col]
-    H1_spec <- H1[col]
-    sum(spec[wl >= (H1_spec - 50) & wl <= (H1_spec + 50)])
-  }, numeric(1))
-  S3 <- S3 / B1
+calc_S9 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  R450 <- as.numeric(object[which(wl == 450), , drop = FALSE])
+  R700 <- as.numeric(object[which(wl == 700), , drop = FALSE])
 
-  # Spectral saturation
-  S2 <- B3 / Rmin # S2
+  # Calc
+  (R700 - R450) / R700
+}
 
-  S6 <- B3 - Rmin # S6
+calc_S10 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  S8 <- computed_vars[["S8"]]
+  S4 <- computed_vars[["S4"]]
 
-  S8 <- S6 / B2 # S8
+  # Calc
+  S8 * S4
+}
 
-  # H2
-  diffsmooth <- apply(object, 2, diff)
+## ----- Hue ----- ##
+
+calc_H1 <- function(object, wl, computed_vars, ...) {
+  # Calc
+  wl[max.col(t(object), ties.method = "first")]
+}
 
-  # Spectra that are monotically increasing or decreasing
+calc_H2 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  diffsmooth <- apply(object, 2, diff)
   incr <- apply(diffsmooth, 2, min) > 0
   decr <- apply(diffsmooth, 2, max) < 0
 
-  lambdabmaxneg <- wl[apply(diffsmooth, 2, which.min)] # H2
+  # Calc
+  lambdabmaxneg <- wl[apply(diffsmooth, 2, which.min)]
   lambdabmaxneg[incr] <- NA
+  lambdabmaxneg
+}
 
-  # S4
-  bmaxneg <- abs(apply(diffsmooth, 2, min)) # S4
-  bmaxneg[incr] <- NA
-
-  # S10
-  S10 <- S8 * bmaxneg # S10
-
-  # H5
-  lambdabmax <- wl[apply(diffsmooth, 2, which.max)] # H5
-  lambdabmax[decr] <- NA
-
-
-  # Add remaining variables to output
-
-  output.mat[, 1] <- B1
-  output.mat[, 2] <- B2
-  output.mat[, 3] <- B3
-  output.mat[, 10] <- S2
-  output.mat[, 11] <- S3
-  output.mat[, 12] <- bmaxneg
-  output.mat[, 13] <- S5
-  output.mat[, 14] <- S6
-  output.mat[, 15] <- S7
-  output.mat[, 16] <- S8
-  output.mat[, 17] <- Carotchroma
-  output.mat[, 18] <- S10
-  output.mat[, 19] <- H1
-  output.mat[, 20] <- lambdabmaxneg
-  output.mat[, 21] <- H3 # Rmid
-  output.mat[, 22] <- H4
-  output.mat[, 23] <- lambdabmax
-
-  # PPB added S1v and S1Y
-
+calc_H3 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  B3 <- computed_vars[["B3"]]
+  Rmin <- computed_vars[["Rmin"]]
+  Rmid <- computed_vars[["Rmid"]]
+  index_Rmid <- vapply(seq_len(ncol(object)), function(x) {
+    which.min(abs(object[, x] - Rmid[x]))
+  }, numeric(1))
 
-  color.var <- data.frame(output.mat, row.names = names(object))
+  # Calc
+  wl[index_Rmid]
+}
 
-  colvarnames <- c(
-    "B1", "B2", "B3", "S1U", "S1V", "S1B", "S1G",
-    "S1Y", "S1R", "S2", "S3", "S4", "S5", "S6", "S7", "S8",
-    "S9", "S10", "H1", "H2", "H3", "H4", "H5"
-  )
+calc_H4 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  lambdamin <- min(wl)
+  lambdamax <- max(wl)
+  segmts <- trunc(quantile(lambdamin:lambdamax, names = FALSE))
+  Q1 <- wl >= segmts[1] & wl <= segmts[2]
+  Q2 <- wl >= segmts[2] & wl <= segmts[3]
+  Q3 <- wl >= segmts[3] & wl <= segmts[4]
+  Q4 <- wl >= segmts[4] & wl <= segmts[5]
+  S5R <- colSums(object[Q4, , drop = FALSE])
+  S5Y <- colSums(object[Q3, , drop = FALSE])
+  S5G <- colSums(object[Q2, , drop = FALSE])
+  S5B <- colSums(object[Q1, , drop = FALSE])
 
-  names(color.var) <- colvarnames
+  # Calc
+  atan2(S5Y - S5B, S5R - S5G)
+}
 
-  if (is.logical(subset)) {
-    if (subset) {
-      color.var <- color.var[c("B2", "S8", "H1")]
-    }
-  } else {
-    # check if any color variables selected don't exist
-    if (all(subset %in% colvarnames)) {
-      color.var <- color.var[subset]
-    } else {
-      stop("Names in ", dQuote("subset"), " do not match color variable names")
-    }
-  }
+calc_H5 <- function(object, wl, computed_vars, ...) {
+  # Dependencies
+  diffsmooth <- apply(object, 2, diff)
+  decr <- apply(diffsmooth, 2, max) < 0
 
-  color.var
+  # Calc
+  lambdabmax <- wl[apply(diffsmooth, 2, which.max)]
+  lambdabmax[decr] <- NA
+  lambdabmax
 }
diff --git a/cran-comments.md b/cran-comments.md
index 43cb5a2ae..1cac0c793 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,17 +1,12 @@
 ## Test environments
 
-* local macOS 12.4 install, R 4.2.1
-* macOS (on GitHub actions), R release
-* ubuntu 20.04 (on GitHub actions), R 4.1.3, R 4.2.1, R-devel
-* Windows (on GitHub actions), R 4.2.1
+* local macOS 13.4.1 install, R 4.3.1
+* macOS (on GitHub actions), R 4.3.1
+* ubuntu 22.04.3 (on GitHub actions), R 4.2.3, R 4.3.1, R-devel
+* Windows (on GitHub actions), R 4.3.1
 * win-builder (R release and devel)
-* rhub::check_for_cran()   
-* rhub::check_on_solaris()  
+* rhub::check_for_cran()     
 
 ## R CMD check results
 
-0 errors | 0 warnings | 0 note
-
-## Comments
-
-This release includes a fix for the math rendering problems identified by CRAN checks, as also kindly communicated by email from Kurt Hornik.
\ No newline at end of file
+0 errors | 0 warnings | 0 note
\ No newline at end of file
diff --git a/inst/CITATION b/inst/CITATION
index fa2ad5293..af9ea5398 100644
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -1,13 +1,12 @@
-citHeader("To cite package pavo in publications use:")
+citHeader("To cite the package pavo in publications please use:")
 
-citEntry(
-  entry = "Article",
+bibentry(
+  bibtype = "Article",
   title = "pavo 2: new tools for the spectral and spatial analysis of colour in R",
-  author = c("Rafael Maia", "Hugo Gruson", "John A. Endler", "Thomas E. White"),
+  author = "Rafael Maia and Hugo Gruson and John A. Endler and Thomas E. White",
   journal = "Methods in Ecology and Evolution",
   year = "2019",
   volume = "10",
   number = "7",
   page = "1097-1107",
-  doi = "10.1111/2041-210X.13174",
-  textVersion = "Maia R., Gruson H., Endler J. A., White T. E. 2019. pavo 2: new tools for the spectral and spatial analysis of colour in R. Methods in Ecology and Evolution, 10(7), 1097-1107")
+  doi = "10.1111/2041-210X.13174")
diff --git a/man/stitch.Rd b/man/stitch.Rd
new file mode 100644
index 000000000..f986a0a5d
--- /dev/null
+++ b/man/stitch.Rd
@@ -0,0 +1,84 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/stitch.rspec.R
+\name{stitch}
+\alias{stitch}
+\alias{stitch.rspec}
+\title{Stitch together two rspec objects}
+\usage{
+stitch(rspec1, rspec2, overlap_method, interp)
+
+\method{stitch}{rspec}(
+  rspec1,
+  rspec2,
+  overlap_method = c("mean", "minimum", "maximum"),
+  interp = TRUE
+)
+}
+\arguments{
+\item{rspec1, rspec2}{(required) \code{rspec} objects of differing wavelength ranges
+to stitch together by row.}
+
+\item{overlap_method}{the method for modifying reflectance values if regions
+of the spectra overlap in their wavelength range. Defaults to \code{mean}.}
+
+\item{interp}{logical argument specifying whether reflectance values should be
+interpolated between the two sets of spectra if their wavelength ranges
+do not overlap. Defaults to \code{TRUE}.}
+}
+\description{
+Stitch (row-wise merge) two \code{rspec} objects of differing wavelength ranges into
+a single \code{rspec} object.
+}
+\examples{
+
+# Simulate a UV-VIS and NIR reflectance spectrum whose wavelength regions
+# slightly overlap then stitch them together, with the overlapping
+# regions being averaged.
+
+# Simulate specs
+reflect1 <- simulate_spec(wl_peak = 550, xlim = c(300, 700))
+reflect2 <- simulate_spec(wl_inflect = 1100, xlim = c(650, 1200))
+
+# Ensure the names of the spectra match
+names(reflect1) <- names(reflect2) <- c('wl', 'sample_1')
+
+# Stitch the spectra together by their wavelength column
+full_spec <- stitch(reflect1, reflect2)
+
+# Plot the resulting spectrum
+plot(full_spec)
+
+# Simulate another set of UV-VIS and NIR spectra. Note two additional complexities,
+# both of which are handled without issue. First, the wavelength ranges are
+# non-overlapping (with a 100 nm gap). We'll keep the default interp = TRUE argument
+# to allow the missing reflectance region to be interpolated. Second, the names of
+# the spectra match, but are in a different order in the two rspec objects. This isn't
+# an issue, as the function can match up the spectra by name irrespective of their
+# ordering
+
+# Simulate UV-VIS and NIR spectra
+reflect_vis <- merge(simulate_spec(wl_peak = 550, xlim = c(300, 700)),
+                     simulate_spec(wl_peak = 550, xlim = c(300, 700)))
+reflect_nir <- merge(simulate_spec(wl_inflect = 1000, xlim = c(800, 1250)),
+                     simulate_spec(wl_inflect = 1100, xlim = c(800, 1250)))
+
+# Ensure the names of the spectra exist in each, albeit in a different order
+names(reflect_vis) <- c('wl', 'sample_1', 'sample_2')
+names(reflect_nir) <- c('wl', 'sample_2', 'sample_1')
+
+# Stitch together by their wavelength column, with missing regions being
+# interpolated
+reflect_vis_nir <- stitch(reflect_vis, reflect_nir)
+
+# Plot the resulting spectrum
+plot(reflect_vis_nir)
+
+}
+\seealso{
+\code{\link[=as.rspec]{as.rspec()}}, \code{\link[=merge.rspec]{merge.rspec()}}
+}
+\author{
+Thomas White \email{thomas.white026@gmail.com}
+
+Hugo Gruson \email{hugo.gruson+R@normalesup.org}
+}
diff --git a/man/summary.rspec.Rd b/man/summary.rspec.Rd
index 6e2b5f22c..1f7da63a2 100644
--- a/man/summary.rspec.Rd
+++ b/man/summary.rspec.Rd
@@ -4,7 +4,7 @@
 \alias{summary.rspec}
 \title{Colourimetric variables}
 \usage{
-\method{summary}{rspec}(object, subset = FALSE, wlmin = NULL, wlmax = NULL, ...)
+\method{summary}{rspec}(object, subset = FALSE, lim = NULL, wlmin = NULL, wlmax = NULL, ...)
 }
 \arguments{
 \item{object}{(required) a data frame, possibly an object of class \code{rspec},
@@ -17,8 +17,10 @@ of the complete output (composed of B2, S8 and H1; the variables described in
 Andersson and Prager 2006) are returned. Finally, a user-specified string of variable
 names can be used in order to filter and show only those variables.}
 
-\item{wlmin, wlmax}{minimum and maximum used to define the range of wavelengths used in
-calculations (default is to use entire range in the \code{rspec} object).}
+\item{lim}{The range of wavelengths used in calculations. The default is to use
+the entire range in the \code{rspec} object (typically equivalent to \code{lim = c(300, 700)}).}
+
+\item{wlmin, wlmax}{Deprecated. Use the \code{lim} argument instead.}
 
 \item{...}{class consistency (ignored)}
 }
@@ -115,10 +117,18 @@ Smoothing affects only B3, S2, S4, S6, S10, H2, and H5 calculation. All other
 variables can be reliably extracted using non-smoothed data.
 }
 \examples{
+# Load data
 data(sicalis)
+
+# Calculate and display all spectral summary variables
 summary(sicalis)
+
+# Calculate only subset of B2, S8 and H1 as per Andersson (1999)
 summary(sicalis, subset = TRUE)
+
+# Calculate user-specified subset of B1 and H4
 summary(sicalis, subset = c("B1", "H4"))
+
 }
 \references{
 Montgomerie R. 2006. Analyzing colors. In Hill, G.E, and McGraw, K.J., eds.
@@ -172,5 +182,9 @@ eastern bluebirds. Animal Behaviour 69:67-72.
 correlated with paternal care in bluethroats? Behavioural Ecology 12:164-170.
 }
 \author{
-Pierre-Paul Bitton \email{bittonp@windsor.ca}, Rafael Maia \email{rm72@zips.uakron.edu}
+Thomas E. White \email{thomas.white026@gmail.com}
+
+Pierre-Paul Bitton \email{bittonp@windsor.ca}
+
+Rafael Maia \email{rm72@zips.uakron.edu}
 }
diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml
index 108c79e51..df9b2a5b0 100644
--- a/pkgdown/_pkgdown.yml
+++ b/pkgdown/_pkgdown.yml
@@ -20,6 +20,7 @@ reference:
   - '`procspec`'
   - '`rimg2cimg`'
   - '`spec2rgb`'
+  - '`stitch.rspec`'
   - '`subset.rspec`'
   - '`simulate_spec`'
 
diff --git a/tests/testthat/test-S3rspec.R b/tests/testthat/test-S3rspec.R
index 813695949..5ee8a0d6b 100644
--- a/tests/testthat/test-S3rspec.R
+++ b/tests/testthat/test-S3rspec.R
@@ -86,25 +86,24 @@ test_that("summary.rspec", {
   expect_named(summary(sicalis, subset = c("B1", "H4")), c("B1", "H4"))
 
   # Different wl ranges
-  expect_warning(summary(sicalis, wlmin = 500), "wavelength range not between")
+  expect_warning(summary(sicalis, lim = c(500, 700)), "wavelength range not between")
   expect_warning(summary(sicalis[1:200, ]), "wavelength range not between")
-  expect_warning(summary(sicalis, wlmax = 600), "wavelength range not between")
-  expect_error(summary(sicalis, wlmin = 200), "wlmin is smaller")
-  expect_error(summary(sicalis, wlmax = 1000), "wlmax is larger")
+  expect_warning(summary(sicalis, lim = c(300, 600)), "wavelength range not between")
+  expect_error(summary(sicalis, lim = c(200, 700)), "smaller than the range")
+  expect_error(summary(sicalis, lim = c(300, 1000)), "larger than the range")
 
   # Test one spectrum rspec object
   one_spec <- sicalis[, c(1, 2)]
   expect_identical(dim(summary(one_spec)), c(1L, 23L))
   expect_warning(
-    expect_length(summary(one_spec, wlmin = 500), 23),
-    "blue chroma"
+    expect_length(summary(one_spec, lim = c(500, 700)), 23), "blue chroma"
   )
 
   # Error if subset vars do not exist
   expect_error(summary(sicalis, subset = "H9"), "do not match color variable names")
 
   # Warning about UV variables if full UV range is not included
-  expect_warning(summary(sicalis, wlmin = 350), "UV-related variables may not be meaningful")
+  expect_warning(summary(sicalis, lim = c(350, 700)), "UV-related variables may not be meaningful")
 })
 
 test_that("plot.rspec", {
diff --git a/tests/testthat/test-hashes.R b/tests/testthat/test-hashes.R
deleted file mode 100644
index 454bcbdd0..000000000
--- a/tests/testthat/test-hashes.R
+++ /dev/null
@@ -1,472 +0,0 @@
-local_edition(2)
-
-skip_on_cran()
-
-data(flowers)
-data(sicalis)
-data(teal)
-
-test_that("coldist", {
-  # JND transform
-  vis.flowers <- vismodel(flowers, visual = "apis")
-  cd.flowers <- coldist(vis.flowers, n = c(1, 1, 1))
-  jnd.flowers <- jnd2xyz(cd.flowers)
-  # expect_equal(digest::sha1(jndrot(jnd2xyz(coldist(vismodel(flowers, achromatic = "bt.dc", relative = FALSE), achromatic = TRUE))), digits = 4),
-  #              "07064d68561bad24d8f02c0413611b5ba49ec53a")
-
-  # Output
-  expect_identical(
-    digest::sha1(coldist(colspace(
-      vismodel(flowers, visual = "canis", achromatic = "ml")
-    ), achromatic = TRUE), digits = 4),
-    "bc460149b2263a857c9d573e77169556fa641f56"
-  )
-  # expect_equal(digest::sha1(coldist(vismodel(flowers, visual = 'canis', achromatic = 'ml'), achromatic = TRUE, n = c(1, 1)), digits = 4),
-  #              "7329a3c550fe1d2939423e4104066c868891914f")
-  expect_identical(
-    digest::sha1(
-      coldist(
-        colspace(vismodel(
-          flowers,
-          visual = "canis", achromatic = "all"
-        )),
-        n = c(1, 2),
-        achromatic = TRUE,
-        subset = "Hibbertia_acicularis"
-      ),
-      digits = 4
-    ),
-    "27ab9af8efe2b1651cd36f8506262f87e2b127a7"
-  )
-  expect_identical(
-    digest::sha1(
-      coldist(
-        colspace(
-          vismodel(
-            flowers,
-            visual = "apis",
-            achromatic = "all",
-            relative = FALSE,
-            vonkries = TRUE
-          ),
-          space = "hexagon"
-        ),
-        n = c(1, 2),
-        achromatic = TRUE,
-        subset = c("Hibbertia_acicularis", "Grevillea_buxifolia")
-      ),
-      digits = 4
-    ),
-    "754c01809100bdacc80d40db2359797f41180c23"
-  )
-  expect_identical(
-    digest::sha1(coldist(colspace(
-      vismodel(flowers, visual = "segment")
-    ), achromatic = TRUE), digits = 4),
-    "d65c018342664ae9c8dca35e715c57dde28de30a"
-  )
-  expect_identical(
-    digest::sha1(coldist(colspace(
-      vismodel(
-        flowers,
-        visual = "cie10",
-        illum = "D65",
-        vonkries = TRUE,
-        relative = FALSE
-      ),
-      "cielab"
-    )), digits = 4),
-    "ab8d1c2eac211561f68759137baa2b5d3005b199"
-  )
-})
-
-test_that("coldist_nolinux", {
-  skip_on_os("linux")
-
-  expect_identical(
-    digest::sha1(coldist(
-      as.matrix(vismodel(flowers, achro = "bt.dc")),
-      qcatch = "Qi",
-      achromatic = TRUE
-    ), digits = 3),
-    "1f797fe87a2e1502080e1c99251b3a768164e7c7"
-  )
-})
-
-test_that("bootcoldist", {
-  skip_on_os("linux")
-
-  # Empirical means
-  data(sicalis)
-  vm.sic <- vismodel(sicalis, visual = "apis", achromatic = "l")
-  gr.sic <- gsub("ind..", "", rownames(vm.sic))
-  bcd.sic <- suppressWarnings(bootcoldist(
-    vm.sic,
-    by = gr.sic,
-    n = c(1, 2, 3),
-    weber = 0.1,
-    weber.achro = 0.1
-  ))
-  expect_identical(
-    digest::sha1(c(bcd.sic[, 1], bcd.sic[, 4])),
-    "9d10fc3ed22e787974464b30a0c5254209522388"
-  )
-})
-
-test_that("special_colspace", {
-  skip_on_os("linux")
-
-  expect_identical(
-    digest::sha1(colspace(
-      vismodel(flowers, visual = "canis", achromatic = "all")
-    ), digits = 3),
-    "41e6dafe465f01a883abd395c00be168aec45b2d"
-  ) # dispace
-  expect_identical(
-    digest::sha1(colspace(
-      vismodel(flowers, visual = "apis", achromatic = "l")
-    ), digits = 3),
-    "f327404a1e417cef90d5c79a5cfd6ecb8c494de4"
-  ) # trispace
-  expect_identical(
-    digest::sha1(colspace(
-      vismodel(flowers, visual = "bluetit", achromatic = "ch.dc")
-    ), digits = 3),
-    "d6d882283d9f005f436066b7da0d5e9a4c4b1e15"
-  ) # tcs
-})
-
-test_that("special_colspace", {
-  expect_identical(
-    digest::sha1(colspace(
-      vismodel(flowers, visual = "musca", achro = "md.r1"),
-      space = "categorical"
-    ), digits = 4),
-    "b20853b3e52a60f2dd17b418a48b681d7f49e7d1"
-  ) # categorical
-  expect_identical(
-    digest::sha1(colspace(
-      vismodel(flowers, visual = "segment", achromatic = "bt.dc"),
-      space = "segment"
-    ), digits = 4),
-    "f47081fbc5f3f896fc50b2223937d91b6f61069e"
-  ) # segment
-  expect_identical(
-    digest::sha1(colspace(
-      vismodel(
-        flowers,
-        visual = "apis",
-        relative = FALSE,
-        qcatch = "Ei",
-        vonkries = TRUE,
-        achromatic = "l"
-      ),
-      space = "coc"
-    ), digits = 4),
-    "d6e5c22dd45d2604c0d2fc16509b8887cb7812d2"
-  ) # coc
-  expect_identical(
-    digest::sha1(
-      colspace(
-        vismodel(
-          flowers,
-          visual = "apis",
-          qcatch = "Ei",
-          vonkries = TRUE,
-          relative = FALSE,
-          achromatic = "l"
-        ),
-        space = "hexagon"
-      ),
-      digits = 4
-    ),
-    "a1fdd24e315413825c94d4caf1164b8be57c8156"
-  ) # hexagon
-  expect_identical(
-    digest::sha1(colspace(
-      vismodel(flowers, visual = "cie10"),
-      space = "ciexyz"
-    ), digits = 4),
-    "99684793a0db286562bff697354496ac3ef0abfb"
-  ) # ciexyz
-  expect_identical(
-    digest::sha1(colspace(
-      vismodel(flowers, visual = "cie10"),
-      space = "cielab"
-    ), digits = 4),
-    "55961f7e22403fba0c0c658868918722befb5f2c"
-  ) # cielab
-  expect_identical(
-    digest::sha1(colspace(
-      vismodel(flowers, visual = "cie10"),
-      space = "cielch"
-    ), digits = 4),
-    "e0ad250b695e97c9ffb53c0303f19d908e33d033"
-  ) # cielch
-
-  # sha1() has no method for the 'table' class
-  # expect_equal(
-  #   digest::sha1(summary(colspace(vismodel(flowers, visual = "cie10"), space = "cielch")), digits = 4),
-  #   "8d9c05ec7ae28b219c4c56edbce6a721bd68af82"
-  # )
-})
-
-test_that("voloverlap()", {
-  tcs.sicalis.C <- subset(colspace(vismodel(sicalis)), "C")
-  tcs.sicalis.T <- subset(colspace(vismodel(sicalis)), "T")
-  tcs.sicalis.B <- subset(colspace(vismodel(sicalis)), "B")
-
-  expect_identical(
-    digest::sha1(
-      voloverlap(tcs.sicalis.T, tcs.sicalis.B, type = "convex"),
-      digits = 4
-    ),
-    "3717422024683f1e3e1bd8dbfe832b177147afce"
-  )
-
-  expect_identical(
-    digest::sha1(
-      voloverlap(tcs.sicalis.T, tcs.sicalis.C, type = "convex"),
-      digits = 4
-    ),
-    "69b323778e83f2e43a91d60326f1e726eb2cd0e4"
-  )
-})
-
-test_that("processing & general", {
-  # Sensdata
-  expect_known_hash(
-    expect_silent(sensdata(
-      illum = "all",
-      bkg = "all",
-      trans = "all"
-    )),
-    "b084b37ec7"
-  )
-
-  # Peakshape
-  expect_known_hash(
-    expect_silent(peakshape(flowers, absolute.min = TRUE)),
-    "7fbaba1738"
-  )
-
-  # Simulate
-  # Ideal
-  expect_known_hash(
-    digest::sha1(summary(simulate_spec(ylim = c(
-      0, 50
-    ))), digits = 4),
-    "88a30943a99d1ea14bc42ac4e4f4d54d"
-  )
-  # Sigmoidd low-high
-  expect_known_hash(
-    digest::sha1(summary(simulate_spec(wl_inflect = 550)), digits = 4),
-    "8e3891ab1c6cf8a10bcab6410fa2f4a0"
-  )
-  # Sigmoid high-low
-  expect_known_hash(
-    digest::sha1(summary(
-      simulate_spec(wl_inflect = 550, ylim = c(100, 0))
-    ), digits = 4),
-    "2b9191b2b7107e2903702b77a80c432d"
-  )
-  # Gaussian
-  expect_known_hash(
-    digest::sha1(summary(simulate_spec(wl_peak = 400)), digits = 4),
-    "9e86606691410fa965031f7b6b180c48"
-  )
-
-  # Merge
-  teal1 <- teal[, c(1, 3:5)]
-  teal2 <- teal[, c(1, 2, 6:12)]
-  expect_known_hash(
-    expect_silent(merge(teal1, teal2, by = "wl")),
-    "02df3eedf3"
-  )
-
-  # Subset
-  vis.sicalis <- vismodel(sicalis)
-  tcs.sicalis <- colspace(vis.sicalis, space = "tcs")
-  expect_identical(
-    digest::sha1(subset(vis.sicalis, "C"), digits = 4),
-    "93ae671c250d2d4f0f5dcf9e714eb497d8baf74f"
-  )
-  expect_identical(
-    digest::sha1(subset(sicalis, "T", invert = TRUE), digits = 4),
-    "332a97ed1c25045b70d871a8686e268d09cefd76"
-  )
-
-  # Summary
-  expect_known_hash(
-    expect_silent(summary(teal)),
-    "c64e1fd403"
-  )
-  expect_known_hash(
-    expect_silent(summary(sicalis)),
-    "66129550f3"
-  )
-})
-
-test_that("images", {
-  suppressWarnings(RNGversion("3.5.0")) # back compatibility for now
-  set.seed(2231)
-
-  papilio <-
-    getimg(system.file("testdata", "images", "butterflies", "papilio.png", package = "pavo"))
-  snakes <-
-    getimg(system.file("testdata", "images", "snakes", package = "pavo"))
-
-  expect_identical(
-    digest::sha1(summary(papilio), digits = 4),
-    "aa1c46d4796c523f51c4e959ac90a692dd8ecfe4"
-  )
-  expect_identical(
-    digest::sha1(summary(snakes), digits = 4),
-    "dd7fc9fd7c41da84c641181a9a3701da74f3c41e"
-  )
-})
-
-test_that("vismodel", {
-  # Output
-  expect_identical(
-    digest::sha1(
-      vismodel(
-        flowers,
-        visual = "canis",
-        achromatic = "all",
-        illum = "bluesky"
-      ),
-      digits = 4
-    ),
-    "61879badc0cb518ebd8f62f9c8838c7b32cb51ff"
-  )
-  expect_identical(
-    digest::sha1(
-      vismodel(
-        flowers,
-        visual = "apis",
-        qcatch = "fi",
-        achromatic = "ml",
-        scale = 10000
-      ),
-      digits = 4
-    ),
-    "4a3539c87d1c672510df68992b9dc6954337a736"
-  )
-  expect_identical(
-    digest::sha1(
-      vismodel(
-        flowers,
-        visual = "bluetit",
-        achromatic = "ch.dc",
-        trans = "bluetit"
-      ),
-      digits = 4
-    ),
-    "fbd9f6b5368f2c81f11ec86a78322e7a14cc7f47"
-  )
-  expect_identical(
-    digest::sha1(
-      vismodel(
-        flowers,
-        visual = "musca",
-        achromatic = "md.r1",
-        relative = FALSE
-      ),
-      digits = 4
-    ),
-    "edcb721c2093c7af40efdae94837c4e01031c8ae"
-  )
-  expect_identical(
-    digest::sha1(
-      vismodel(
-        flowers,
-        visual = "apis",
-        relative = FALSE,
-        qcatch = "Ei",
-        bkg = "green",
-        vonkries = TRUE,
-        achromatic = "l"
-      ),
-      digits = 4
-    ),
-    "a6bc51f272c930a4ac9e69a1851eca16f5a3a1a0"
-  )
-  expect_identical(
-    digest::sha1(vismodel(flowers, visual = "cie10"), digits = 4),
-    "fc5f5f2f11fefdcff1bbdd264e28d520f0812712"
-  )
-
-  # Attributes
-  if (getRversion() < "4.0.0") {
-    expect_identical(
-      digest::sha1(attributes(
-        vismodel(
-          flowers,
-          visual = "canis",
-          achromatic = "all",
-          illum = "bluesky"
-        )
-      ), digits = 4),
-      "0f788526e4db68c9921e441066779146f8a4c377"
-    )
-    expect_identical(
-      digest::sha1(attributes(
-        vismodel(
-          flowers,
-          visual = "apis",
-          qcatch = "fi",
-          achromatic = "ml",
-          scale = 10000
-        )
-      ), digits = 4),
-      "b9d488a8e36bca04a66e4513e781c21b66c10ce9"
-    )
-    expect_identical(
-      digest::sha1(attributes(
-        vismodel(
-          flowers,
-          visual = "bluetit",
-          achromatic = "ch.dc",
-          trans = "bluetit"
-        )
-      ), digits = 4),
-      "0a3fb5b639d4c4224cf91045d5f8a13cc06f8550"
-    )
-  } else {
-    expect_identical(
-      digest::sha1(attributes(
-        vismodel(
-          flowers,
-          visual = "canis",
-          achromatic = "all",
-          illum = "bluesky"
-        )
-      ), digits = 4),
-      "397eda31c6948884e09ef58a1b8f0ec5d4f3401c"
-    )
-    expect_identical(
-      digest::sha1(attributes(
-        vismodel(
-          flowers,
-          visual = "apis",
-          qcatch = "fi",
-          achromatic = "ml",
-          scale = 10000
-        )
-      ), digits = 4),
-      "60744da46c20782fa8af0bf2a93316ed0d2e6e9d"
-    )
-    expect_identical(
-      digest::sha1(attributes(
-        vismodel(
-          flowers,
-          visual = "bluetit",
-          achromatic = "ch.dc",
-          trans = "bluetit"
-        )
-      ), digits = 4),
-      "7bca56baefdd2a42fcb59c26614fe0bbec326865"
-    )
-  }
-})
diff --git a/tests/testthat/test-processing.R b/tests/testthat/test-processing.R
index 00f622a03..dc8d8e0b1 100644
--- a/tests/testthat/test-processing.R
+++ b/tests/testthat/test-processing.R
@@ -86,6 +86,46 @@ test_that("Aggregation", {
   expect_error(aggspec(teal, by = 7), "by not a multiple")
 })
 
+test_that("Stitch", {
+  
+  # Overlapping ranges
+  r1 <- simulate_spec(wl_peak = 550, xlim = c(300, 700))
+  r2 <- simulate_spec(wl_inflect = 1100, xlim = c(650, 1200))
+  names(r1) <- names(r2) <- c('wl', 'sample_1')
+  r_stitch <- stitch(r1, r2)
+  
+  expect_equal(dim(r_stitch), c(901, 2))
+  expect_equal(r_stitch$wl, 300:1200)
+  expect_equal(summary(r_stitch)$H5, 520)
+  expect_equal(summary(r_stitch)$S6, 100)
+  expect_equal(summary(r_stitch)$B3, 100)
+
+  # Non-overlapping ranges
+  r_vis <- merge(simulate_spec(wl_peak = 550, xlim = c(300, 700)),
+                 simulate_spec(wl_peak = 550, xlim = c(300, 700)))
+  r_nir <- merge(simulate_spec(wl_inflect = 1000, xlim = c(800, 1250)),
+                 simulate_spec(wl_inflect = 1100, xlim = c(800, 1250)))
+  names(r_vis) <- c('wl', 'sample_1', 'sample_2')  # Names match, different order
+  names(r_nir) <- c('wl', 'sample_2', 'sample_1')
+  r_vis_nir <- stitch(r_vis, r_nir)
+  
+  expect_equal(dim(r_vis_nir), c(951, 3))
+  expect_equal(r_vis_nir$wl, 300:1250)
+  expect_equal(summary(r_vis_nir)$H5, c(520, 520))
+  expect_equal(summary(r_vis_nir)$S6, c(100, 100))
+  expect_equal(summary(r_vis_nir)$B3, c(100, 100))
+  
+  # Errors and warnings
+  names(r_vis) <- c('wl', 'sample_1', 'sample_2')
+  names(r_nir) <- c('wl', 'sample_3', 'sample_4')
+  expect_error(stitch(r_vis, r_nir), "matching name")
+  
+  names(r_vis) <- c('wl', 'sample_1', 'sample_2')
+  names(r_nir) <- c('wl', 'sample_1', 'sample_3')
+  expect_warning(stitch(r_vis, r_nir), "Not all spectra are present")
+  
+})
+
 test_that("Convert", {
   # Flux/irrad
   illum <- sensdata(illum = "forestshade")
diff --git a/tests/testthat/test-regression.R b/tests/testthat/test-regression.R
new file mode 100644
index 000000000..9096f85dc
--- /dev/null
+++ b/tests/testthat/test-regression.R
@@ -0,0 +1,219 @@
+data(flowers)
+data(sicalis)
+data(teal)
+
+test_that("coldist", {
+  # JND transform
+  vis.flowers <- vismodel(flowers, visual = "apis")
+  cd.flowers <- coldist(vis.flowers, n = c(1, 1, 1))
+  jnd.flowers <- jnd2xyz(cd.flowers)
+  result1 <- jndrot(jnd2xyz(coldist(vismodel(flowers, achromatic = "bt.dc", relative = FALSE), achromatic = TRUE)))
+  
+  expect_equal(dim(result1), c(36, 4))
+  expect_equal(round(mean(unlist(result1)^3), 3), -132.915)
+
+  # Output
+  result2 <- coldist(colspace(vismodel(flowers, visual = "canis", achromatic = "ml")), achromatic = TRUE)
+  expect_equal(dim(result2), c(630, 4))
+  expect_equal(round(mean(unlist(result2[,3:4])), 4), 0.8693)
+  
+  result3 <- coldist(colspace(vismodel(flowers, visual = "canis", achromatic = "all")),
+             n = c(1, 2), achromatic = TRUE, subset = "Hibbertia_acicularis")
+  expect_equal(dim(result3), c(35, 4))
+  expect_equal(round(mean(unlist(result3[, 3:4])), 4), 0.5388)
+  
+  result4 <- coldist(colspace(vismodel(flowers, visual = "apis", achromatic = "all",
+                    relative = FALSE, vonkries = TRUE), space = "hexagon"),
+             n = c(1, 2), achromatic = TRUE, subset = c("Hibbertia_acicularis", "Grevillea_buxifolia"))
+  expect_equal(dim(result4), c(1, 4))
+  expect_equal(round(mean(unlist(result4[, 3:4])), 4), 0.464)
+  
+  result5 <- coldist(colspace(vismodel(flowers, visual = "segment")), achromatic = TRUE)
+  expect_equal(dim(result5), c(630, 4))
+  expect_equal(round(mean(unlist(result5[, 3:4])), 4), 0.2685)
+  
+  result6 <- coldist(colspace(vismodel(flowers, 
+                                       visual = "cie10", 
+                                       illum = "D65", 
+                                       vonkries = TRUE, 
+                                       relative = FALSE), "cielab"))
+  expect_equal(dim(result6), c(630, 4))
+  expect_equal(round(mean(unlist(result6[, 3])), 4), 28.7164)
+  
+  result7 <- coldist(as.matrix(vismodel(flowers, achro = "bt.dc")), qcatch = "Qi", achromatic = TRUE)
+  expect_equal(dim(result7), c(630, 4))
+  expect_equal(round(mean(unlist(result7[, 3])), 4), 11.7606)
+  
+})
+
+test_that("bootcoldist", {
+
+  # Empirical means
+  data(sicalis)
+  vm.sic <- vismodel(sicalis, visual = "apis", achromatic = "l")
+  gr.sic <- gsub("ind..", "", rownames(vm.sic))
+  bcd.sic <- suppressWarnings(bootcoldist(
+    vm.sic,
+    by = gr.sic,
+    n = c(1, 2, 3),
+    weber = 0.1,
+    weber.achro = 0.1
+  ))
+  expect_equal(
+    unname(round(c(bcd.sic[, 1], bcd.sic[, 4]), 4)), 
+  c(4.5760, 1.1340, 5.6378, 8.2510, 0.0124, 8.2387))
+})
+
+test_that("special_colspace", {
+
+  # Dispace
+  result8 <- colspace(vismodel(flowers, visual = "canis", achromatic = "all"))
+  expect_equal(dim(result8), c(36, 5))
+  expect_equal(round(mean(unlist(result8)), 4), 0.4061)
+  
+  # Trispace
+  result9 <- colspace(vismodel(flowers, visual = "apis", achromatic = "l"))
+  expect_equal(dim(result9), c(36, 8))
+  expect_equal(round(mean(unlist(result9)), 4), 0.0876)
+  
+  # tcs
+  result10 <- colspace(vismodel(flowers, visual = "bluetit", achromatic = "ch.dc"))
+  expect_equal(dim(result10), c(36, 17))
+  expect_equal(round(mean(unlist(result10)), 4), 0.0923)
+  
+  # categorical
+  result11 <- colspace(vismodel(flowers, visual = "musca", achro = "md.r1"), space = "categorical")
+  expect_equal(dim(result11), c(36, 10))
+  expect_equal(round(mean(unlist(result11[, -7])), 4), -0.0949)
+  
+  # segment
+  result12 <- colspace(vismodel(flowers, visual = "segment", achromatic = "bt.dc"), space = "segment")
+  expect_equal(dim(result12), c(36, 9))
+  expect_equal(round(mean(unlist(result12)), 4), 20.4668)
+  
+  # coc
+  result13 <- colspace(vismodel(flowers, visual = "apis", relative = FALSE, qcatch = "Ei", 
+                       vonkries = TRUE, achromatic = "l"), space = "coc")
+  expect_equal(dim(result13), c(36, 7))
+  expect_equal(round(mean(unlist(result13)), 4), 0.5304)
+  
+  # hexagon
+  result14 <- colspace(vismodel(flowers, visual = "apis", qcatch = "Ei", 
+                                vonkries = TRUE, relative = FALSE, achromatic = "l"), space = "hexagon")
+  expect_equal(dim(result14), c(36, 10))
+  expect_equal(round(mean(unlist(result14[-9])), 4), 20.3489)
+  
+  # ciezyx
+  result15 <- colspace(vismodel(flowers, visual = "cie10"), space = "ciexyz")
+  expect_equal(dim(result15), c(36, 6))
+  expect_equal(round(mean(unlist(result15)), 4), 0.3596)
+  
+  # cielab
+  result16 <- colspace(vismodel(flowers, visual = "cie10"), space = "cielab")
+  expect_equal(dim(result16), c(36, 6))
+  expect_equal(round(mean(unlist(result16)), 4), 9.5446)
+  
+  # cielch
+  result17 <- colspace(vismodel(flowers, visual = "cie10"), space = "cielch")
+  expect_equal(dim(result17), c(36, 8))
+  expect_equal(round(mean(unlist(result17)), 4), 26.2162)
+  
+})
+
+
+test_that("voloverlap()", {
+  tcs.sicalis.C <- subset(colspace(vismodel(sicalis)), "C")
+  tcs.sicalis.T <- subset(colspace(vismodel(sicalis)), "T")
+  tcs.sicalis.B <- subset(colspace(vismodel(sicalis)), "B")
+
+  expect_equal(round(mean(unlist(voloverlap(tcs.sicalis.T, tcs.sicalis.B, type = "convex"))), 4), 0.0395)
+  expect_equal(round(mean(unlist(voloverlap(tcs.sicalis.T, tcs.sicalis.C, type = "convex"))[1:2]^-1), 4), 201959.93)
+  
+})
+
+test_that("processing & general", {
+  # Sensdata
+  result18 <- sensdata(illum = "all", bkg = "all", trans = "all")
+  expect_equal(dim(result18), c(401, 7))
+  expect_equal(round(mean(unlist(result18)), 4), 72.8731)
+
+  # Peakshape
+  result19 <- peakshape(flowers, absolute.min = TRUE)
+  expect_equal(dim(result19), c(36, 7))
+  expect_equal(round(mean(unlist(result19[, 2:6]), na.rm = TRUE), 4), 287.7377)
+
+  # Simulate
+  # Ideal
+  result20 <- summary(simulate_spec(ylim = c(0, 50)))
+  expect_equal(dim(result20), c(1, 23))
+  expect_equal(round(mean(unlist(result20)), 4), 928.3311)
+  
+  # Sigmoidd low-high
+  result21 <- summary(simulate_spec(wl_inflect = 550))
+  expect_equal(dim(result21), c(1, 23))
+  expect_equal(round(mean(unlist(subset(result21, select = -c(S2))), na.rm = TRUE), 4), 1479.6564)
+  
+  # Sigmoid high-low
+  result22 <- summary(simulate_spec(wl_inflect = 550, ylim = c(100, 0)))
+  expect_equal(dim(result22), c(1, 23))
+  expect_equal(round(mean(unlist(subset(result22, select = -c(S2, S9))), na.rm = TRUE), 4), 1885.9745)
+  
+  # Gaussian
+  result23 <- summary(simulate_spec(wl_peak = 400))
+  expect_equal(dim(result23), c(1, 23))
+  expect_equal(round(mean(unlist(subset(result23, select = -c(S2, S9))), na.rm = TRUE), 4), 694.4126)
+
+  # Merge
+  teal1 <- teal[, c(1, 3:5)]
+  teal2 <- teal[, c(1, 2, 6:12)]
+  expect_equal(round(mean(unlist(merge(teal1, teal2, by = "wl"))), 4), 51.1874)
+
+  # Subset
+  vis.sicalis <- vismodel(sicalis)
+  tcs.sicalis <- colspace(vis.sicalis, space = "tcs")
+  expect_equal(round(mean(unlist(subset(vis.sicalis, "C")), na.rm = TRUE), 4), 0.25)
+  expect_equal(round(mean(unlist(subset(sicalis, "T", invert = TRUE)), na.rm = TRUE), 4), 40.1018)
+
+  # Summary
+  expect_equal(round(mean(unlist(summary(teal))), 4), 287.9638)
+  expect_equal(round(mean(unlist(subset(summary(sicalis), select = -c(S2)))), 4), 313.7217)
+  
+})
+
+test_that("images", {
+
+  papilio <-
+    getimg(system.file("testdata", "images", "butterflies", "papilio.png", package = "pavo"))
+  snakes <-
+    getimg(system.file("testdata", "images", "snakes", package = "pavo"))
+
+  expect_equal(unlist(unname(summary(papilio))), 'papilio')
+  expect_equal(unlist(unname(summary(snakes))), c('snake_01', 'snake_02'))
+})
+
+test_that("vismodel", {
+  result24 <- vismodel(flowers, visual = "canis", achromatic = "all", illum = "bluesky")
+  expect_equal(dim(result24), c(36, 3))
+  expect_equal(round(mean(unlist(result24)), 4), 0.4572)
+  
+  result25 <- vismodel(flowers, visual = "apis", qcatch = "fi", achromatic = "ml", scale = 10000)
+  expect_equal(dim(result25), c(36, 4))
+  expect_equal(round(mean(unlist(result25)), 4), 2.3597)
+  
+  result26 <- vismodel(flowers, visual = "bluetit", achromatic = "ch.dc", trans = "bluetit")
+  expect_equal(dim(result26), c(36, 5))
+  expect_equal(round(mean(unlist(result26)), 4), 0.2832)
+  
+  result27 <- vismodel(flowers, visual = "musca", achromatic = "md.r1", relative = FALSE)
+  expect_equal(dim(result27), c(36, 5))
+  expect_equal(round(mean(unlist(result27)), 4), 0.1904)
+  
+  result28 <- vismodel(flowers, visual = "apis", relative = FALSE, qcatch = "Ei", bkg = "green",
+        vonkries = TRUE, achromatic = "l")
+  expect_equal(dim(result28), c(36, 4))
+  expect_equal(round(mean(unlist(result28)), 4), 0.5457)
+  
+  result29 <- vismodel(flowers, visual = "cie10")
+  expect_equal(dim(result29), c(36, 4))
+  expect_equal(round(mean(unlist(result29), na.rm = TRUE), 4), 0.3859)
+})