diff --git a/NEWS.md b/NEWS.md
index b5f3fd03c..12e940079 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,23 @@ stability of steady states.
 
 ## Bug fixes
 
+- `plotCDF()`, `plotCDF2()`, `plotSpectra2()`, and `plotSpectraRelative()` now
+  support `size_axis = "l"` together with `total = TRUE`. The total spectrum and
+  cumulative distributions are summed across species on the length axis.
+
+- `ArrayResourceBySize` and `ArrayTimeByResourceBySize` plotting methods
+  (`plot2()`, `plotRelative()`, `addPlot()`) now support `size_axis = "l"` and
+  `per_log_size`, providing symmetry with species size-spectrum plotting methods.
+
+- `plotFeedingLevel()` and `plotlyFeedingLevel()` have been refactored to
+  delegate directly to `plot()` / `plotHover()` on `ArraySpeciesBySize` while
+  maintaining full support for `include_critical = TRUE` and non-clipping
+  proportion coordinate scaling.
+
+- `plotYield()` for a single simulation now delegates directly to
+  `plot(getYield(object))`. The `sim2` argument of `plotYield()` is deprecated in
+  favour of `plot2(getYield(sim1), getYield(sim2))`.
+
 - `plot(getFluxGradient(params), size_axis = "l")` now converts its values to a
   density with respect to length, and labels them `cm^-1/year`. The flux
   gradient is a rate of change of a number density, but it was not recognised as
diff --git a/R/ArrayResourceBySize-class.R b/R/ArrayResourceBySize-class.R
index 3b548ab3b..77ca51a6f 100644
--- a/R/ArrayResourceBySize-class.R
+++ b/R/ArrayResourceBySize-class.R
@@ -239,28 +239,41 @@ plot2.ArrayResourceBySize <- function(x, y, name1 = "First", name2 = "Second",
                                       ylim = c(NA, NA),
                                       total = FALSE, background = TRUE,
                                       y_ticks = 6,
-                                      wlim = c(NA, NA), ...) {
+                                      wlim = c(NA, NA), llim = c(NA, NA),
+                                      size_axis = c("w", "l"),
+                                      per_log_size = NULL, ...) {
     check_plot2_compatible(x, y, "ArrayResourceBySize")
     compare_array_metadata(x, y)
     warn_unused_resource_args(species, total, background)
+    size_axis <- plot_size_axis(size_axis)
+    check_per_log_size(x, per_log_size)
+    log_y <- array_log_y(x, log_y, log, !missing(log_y))
     log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     log_x <- log_axes$log_x
     log_y <- log_axes$log_y
     assert_that(length(wlim) == 2,
+                length(llim) == 2,
                 length(ylim) == 2)
 
     params <- attr(x, "params")
-    y_label <- array_y_label(x, default = "Value")
+    y_label <- array_y_label(x, default = "Value", size_axis = size_axis,
+                             per_log_size = per_log_size)
     plot_dat1 <- prepare_ArrayResourceBySize_plot_data(x, wlim = wlim)
     plot_dat2 <- prepare_ArrayResourceBySize_plot_data(y, wlim = wlim)
 
+    ylim <- array_ylim(x, ylim, log_y, c(plot_dat1[[2]], plot_dat2[[2]]))
+
     plotComparisonDataFrame(plot_dat1, plot_dat2, params,
                             name1 = name1, name2 = name2,
-                            xlab = "Weight (g)", ylab = y_label,
+                            xlab = plot_size_xlab(size_axis), ylab = y_label,
                             xtrans = if (log_x) "log10" else "identity",
                             ytrans = if (log_y) "log10" else "identity",
-                            xlim = wlim, ylim = ylim,
-                            y_ticks = y_ticks, legend_var = "Legend")
+                            xlim = plot_size_xlim(wlim, size_axis, llim),
+                            ylim = ylim,
+                            y_ticks = y_ticks, legend_var = "Legend",
+                            size_axis = size_axis,
+                            density_wrt = array_density_wrt(x),
+                            per_log_size = per_log_size)
 }
 
 #' @rdname plotRelative
@@ -271,11 +284,15 @@ plotRelative.ArrayResourceBySize <- function(x, y, species = NULL,
                                              ylim = c(NA, NA),
                                              total = FALSE,
                                              background = TRUE,
-                                             wlim = c(NA, NA), ...) {
+                                             wlim = c(NA, NA),
+                                             llim = c(NA, NA),
+                                             size_axis = c("w", "l"), ...) {
     check_plot2_compatible(x, y, "ArrayResourceBySize")
     compare_array_metadata(x, y)
     warn_unused_resource_args(species, total, background)
+    size_axis <- plot_size_axis(size_axis)
     assert_that(length(wlim) == 2,
+                length(llim) == 2,
                 length(ylim) == 2)
 
     params <- attr(x, "params")
@@ -283,10 +300,11 @@ plotRelative.ArrayResourceBySize <- function(x, y, species = NULL,
     plot_dat2 <- prepare_ArrayResourceBySize_plot_data(y, wlim = wlim)
 
     plotRelativeDataFrame(plot_dat1, plot_dat2, params,
-                          xlab = "Weight (g)",
+                          xlab = plot_size_xlab(size_axis),
                           xtrans = if (log_x) "log10" else "identity",
-                          xlim = wlim, ylim = ylim,
-                          legend_var = "Legend")
+                          xlim = plot_size_xlim(wlim, size_axis, llim),
+                          ylim = ylim,
+                          legend_var = "Legend", size_axis = size_axis)
 }
 
 #' @rdname addPlot
@@ -299,7 +317,10 @@ addPlot.ArrayResourceBySize <- function(plot, x, species = NULL,
                                         linetype = "dashed",
                                         linewidth = 0.8,
                                         alpha = 1,
-                                        wlim = c(NA, NA), ...) {
+                                        wlim = c(NA, NA),
+                                        llim = c(NA, NA),
+                                        size_axis = c("w", "l"),
+                                        per_log_size = NULL, ...) {
     if (!inherits(plot, "ggplot")) {
         stop("The `plot` argument must be a ggplot object.")
     }
@@ -307,14 +328,25 @@ addPlot.ArrayResourceBySize <- function(plot, x, species = NULL,
                 is.number(alpha),
                 alpha >= 0,
                 alpha <= 1,
-                length(wlim) == 2)
+                length(wlim) == 2,
+                length(llim) == 2)
     warn_unused_resource_args(species, total, background)
+    size_axis <- plot_size_axis(size_axis)
+    check_per_log_size(x, per_log_size)
 
     plot <- deep_copy(plot)
     plot_dat <- prepare_ArrayResourceBySize_plot_data(x, wlim = wlim)
+    params <- attr(x, "params")
+    plot_dat <- convert_plot_density_axis(plot_dat, params, size_axis,
+                                          density_wrt = array_density_wrt(x),
+                                          per_log_size = per_log_size)
+    if (identical(size_axis, "l")) {
+        plot_dat <- filter_plot_length_limits(plot_dat, llim)
+    }
+    x_var <- plot_size_x_var(size_axis)
     y_var <- names(plot_dat)[2]
-    check_addPlot_compatible(plot, x_var = "w", y_var = y_var,
-                             units = attr(x, "units"))
+    check_addPlot_compatible(plot, x_var = x_var, y_var = y_var,
+                             units = array_units(x, size_axis, per_log_size))
 
     # A resource array is a single line, so there is nothing to distinguish by
     # colour. Mapping colour to the "Resource" legend level would rely on the
@@ -330,7 +362,7 @@ addPlot.ArrayResourceBySize <- function(plot, x, species = NULL,
         }
     }
 
-    mapping <- aes(x = .data[["w"]], y = .data[[y_var]],
+    mapping <- aes(x = .data[[x_var]], y = .data[[y_var]],
                    group = .data[["Species"]])
     if (is.null(linetype)) {
         mapping$linetype <- rlang::quo(.data[["Legend"]])
@@ -656,7 +688,10 @@ plot2.ArrayTimeByResourceBySize <- function(x, y, name1 = "First",
                                             total = FALSE, background = TRUE,
                                             y_ticks = 6,
                                             time = NULL,
-                                            wlim = c(NA, NA), ...) {
+                                            wlim = c(NA, NA),
+                                            llim = c(NA, NA),
+                                            size_axis = c("w", "l"),
+                                            per_log_size = NULL, ...) {
     check_plot2_compatible(x, y, "ArrayTimeByResourceBySize")
     slice1 <- ArrayTimeByResourceBySize_slice(x, time = time)
     slice2 <- ArrayTimeByResourceBySize_slice(y, time = time)
@@ -665,7 +700,8 @@ plot2.ArrayTimeByResourceBySize <- function(x, y, name1 = "First",
                               species = species, log_x = log_x, log_y = log_y,
                               log = log, ylim = ylim, total = total,
                               background = background, y_ticks = y_ticks,
-                              wlim = wlim, ...)
+                              wlim = wlim, llim = llim, size_axis = size_axis,
+                              per_log_size = per_log_size, ...)
 }
 
 #' @rdname plotRelative
@@ -677,14 +713,17 @@ plotRelative.ArrayTimeByResourceBySize <- function(x, y, species = NULL,
                                                    total = FALSE,
                                                    background = TRUE,
                                                    time = NULL,
-                                                   wlim = c(NA, NA), ...) {
+                                                   wlim = c(NA, NA),
+                                                   llim = c(NA, NA),
+                                                   size_axis = c("w", "l"), ...) {
     check_plot2_compatible(x, y, "ArrayTimeByResourceBySize")
     slice1 <- ArrayTimeByResourceBySize_slice(x, time = time)
     slice2 <- ArrayTimeByResourceBySize_slice(y, time = time)
 
     plotRelative.ArrayResourceBySize(slice1, slice2, species = species,
                                      log_x = log_x, ylim = ylim, total = total,
-                                     background = background, wlim = wlim, ...)
+                                     background = background, wlim = wlim,
+                                     llim = llim, size_axis = size_axis, ...)
 }
 
 #' @rdname addPlot
@@ -698,12 +737,17 @@ addPlot.ArrayTimeByResourceBySize <- function(plot, x, species = NULL,
                                               linewidth = 0.8,
                                               alpha = 1,
                                               time = NULL,
-                                              wlim = c(NA, NA), ...) {
+                                              wlim = c(NA, NA),
+                                              llim = c(NA, NA),
+                                              size_axis = c("w", "l"),
+                                              per_log_size = NULL, ...) {
     slice <- ArrayTimeByResourceBySize_slice(x, time = time)
     addPlot.ArrayResourceBySize(plot, slice, species = species, total = total,
                                 background = background, colour = colour,
                                 linetype = linetype, linewidth = linewidth,
-                                alpha = alpha, wlim = wlim, ...)
+                                alpha = alpha, wlim = wlim, llim = llim,
+                                size_axis = size_axis,
+                                per_log_size = per_log_size, ...)
 }
 
 #' @rdname animate
diff --git a/R/ArraySpeciesBySize-class.R b/R/ArraySpeciesBySize-class.R
index 8d3264c0f..0f1a3d39f 100644
--- a/R/ArraySpeciesBySize-class.R
+++ b/R/ArraySpeciesBySize-class.R
@@ -534,6 +534,7 @@ plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second",
     compare_array_metadata(x, y)
     size_axis <- plot_size_axis(size_axis)
     check_per_log_size(x, per_log_size)
+    log_y <- array_log_y(x, log_y, log, !missing(log_y))
     log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     log_x <- log_axes$log_x
     log_y <- log_axes$log_y
@@ -551,6 +552,8 @@ plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second",
         y, species = species, all.sizes = all.sizes, wlim = wlim,
         total = total, background = background)
 
+    ylim <- array_ylim(x, ylim, log_y, c(plot_dat1[[2]], plot_dat2[[2]]))
+
     plotComparisonDataFrame(plot_dat1, plot_dat2, params,
                             name1 = name1, name2 = name2,
                             xlab = plot_size_xlab(size_axis), ylab = y_label,
diff --git a/R/ArrayTimeBySpecies-class.R b/R/ArrayTimeBySpecies-class.R
index c38ff75d7..7110bd94a 100644
--- a/R/ArrayTimeBySpecies-class.R
+++ b/R/ArrayTimeBySpecies-class.R
@@ -281,6 +281,7 @@ plot2.ArrayTimeBySpecies <- function(x, y, name1 = "First", name2 = "Second",
                                      tlim = c(NA, NA), ...) {
     check_plot2_compatible(x, y, "ArrayTimeBySpecies")
     compare_array_metadata(x, y)
+    log_y <- array_log_y(x, log_y, log, !missing(log_y))
     log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     log_x <- log_axes$log_x
     log_y <- log_axes$log_y
@@ -294,6 +295,8 @@ plot2.ArrayTimeBySpecies <- function(x, y, name1 = "First", name2 = "Second",
         y, species = species, tlim = tlim,
         ylim = ylim, total = total, background = background)
 
+    ylim <- array_ylim(x, ylim, log_y, c(plot_dat1[[2]], plot_dat2[[2]]))
+
     plotComparisonDataFrame(plot_dat1, plot_dat2, params,
                             name1 = name1, name2 = name2,
                             xlab = "Year", ylab = y_label,
diff --git a/R/ArrayTimeBySpeciesBySize-class.R b/R/ArrayTimeBySpeciesBySize-class.R
index 50b124780..a7a876a3c 100644
--- a/R/ArrayTimeBySpeciesBySize-class.R
+++ b/R/ArrayTimeBySpeciesBySize-class.R
@@ -393,6 +393,7 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
                 length(wlim) == 2, length(llim) == 2, length(ylim) == 2)
     size_axis <- plot_size_axis(size_axis)
     check_per_log_size(x, per_log_size)
+    log_y <- array_log_y(x, log_y, log, !missing(log_y))
     log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     log_x <- log_axes$log_x
     log_y <- log_axes$log_y
diff --git a/R/plots.R b/R/plots.R
index 210b6d163..d10df4a3d 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -218,9 +218,8 @@ plotDataFrame <- function(frame, params, style = "line", xlab = waiver(),
     if (ytrans == "log10") ybreaks <- log_breaks(n = y_ticks)
 
     # Set up axis limits. NA values mean auto-scale to data range.
-    # The reason why below `group = species` is included in `ggplot()`
-    # rather than in `geom_line` is because that puts it first in the
-    # plotly tooltips, due to a bug in plotly.
+    # The reason why below `group = species` is included first in `ggplot()`
+    # is because that puts it first in the plotly tooltips, due to a bug in plotly.
     p <- ggplot(frame, aes(group = .data[[group_var]])) +
         scale_y_continuous(trans = ytrans, breaks = ybreaks,
                            labels = prettyNum, name = ylab,
@@ -1295,7 +1294,7 @@ plotBiomass.MizerSim <- function(object, species = NULL,
         lifecycle::deprecate_warn("2.6.0", "plotBiomass(end_time)", "plotBiomass(tlim)")
         tlim[[2]] <- end_time
     }
-    log_axes <- parseTimePlotLog(log, log_x = log_x, log_y = log_y)
+    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     bm <- getBiomass(object, use_cutoff = use_cutoff,
                      min_w = min_w, max_w = max_w,
                      min_l = min_l, max_l = max_l)
@@ -1414,83 +1413,59 @@ plotYield.MizerSim <- function(object, sim2 = NULL,
                       ylim = c(NA, NA), tlim = c(NA, NA),
                       highlight = NULL, return_data = FALSE,
                       ...) {
-    log_axes <- parseTimePlotLog(log, log_x = log_x, log_y = log_y)
+    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     assert_that(is(object, "MizerSim"),
                 is.flag(total),
                 is.flag(log_axes$log_x),
                 is.flag(log_axes$log_y),
                 is.flag(return_data))
-    params <- object@params
-    species <- valid_species_arg(object, species, error_on_empty = TRUE)
     if (is.null(sim2)) {
         y <- getYield(object, ...)
-        times <- as.numeric(rownames(y))
-        if (!is.na(tlim[1])) {
-            y <- y[times >= tlim[1], , drop = FALSE]
-            times <- as.numeric(rownames(y))
-        }
-        if (!is.na(tlim[2])) {
-            y <- y[times <= tlim[2], , drop = FALSE]
-        }
-        y_total <- rowSums(y)
-        y <- y[, (as.character(dimnames(y)[[2]]) %in% species),
-               drop = FALSE]
-        if (total) {
-            # Include total
-            y <- cbind(y, "Total" = y_total)
-        }
-        plot_dat <- reshape2::melt(y, varnames = c("Year", "Species"),
-                                   value.name = "Yield")
-        plot_dat <- subset(plot_dat, plot_dat$Yield > 0)
-        # plotDataFrame() needs the columns in a particular order
-        plot_dat <- plot_dat[, c(1, 3, 2)]
-
-        if (nrow(plot_dat) == 0) {
-            warning("There is no yield to include.")
-        }
-        if (return_data) return(plot_dat)
-
-        plotDataFrame(plot_dat, params,
-                      ylab = "Yield [g/year]",
-                      xtrans = ifelse(log_axes$log_x, "log10", "identity"),
-                      ytrans = ifelse(log_axes$log_y, "log10", "identity"),
-                      ylim = ylim,
-                      highlight = highlight)
-    } else {
-        # We need to combine two plots
-        if (!all(dimnames(object@n)$time == dimnames(sim2@n)$time)) {
-            stop("The two simulations do not have the same times")
-        }
-        ym <- plotYield(object, species = species,
-                            tlim = tlim, total = total,
-                            log_x = log_axes$log_x, log_y = log_axes$log_y,
-                            ylim = ylim,
-                            highlight = highlight, return_data = TRUE, ...)
-        ym2 <- plotYield(sim2, species = species,
-                            tlim = tlim, total = total,
-                            log_x = log_axes$log_x, log_y = log_axes$log_y,
-                            ylim = ylim,
-                            highlight = highlight, return_data = TRUE, ...)
-        ym$Simulation <- rep(1, nrow(ym)) # We don't use recycling because that
-                                          # fails when there are zero rows.
-        ym2$Simulation <- rep(2, nrow(ym2))
-        ym <- rbind(ym, ym2)
-
-        if (return_data) return(ym)
-
-        plotDataFrame(ym, params,
-                      ylab = "Yield [g/year]",
-                      xtrans = ifelse(log_axes$log_x, "log10", "identity"),
-                      ytrans = ifelse(log_axes$log_y, "log10", "identity"),
-                      ylim = ylim,
-                      highlight = highlight, wrap_var = "Simulation")
-    }
+        attr(y, "value_name") <- "Yield"
+        return(plot(y, species = species,
+                    tlim = tlim, ylim = ylim,
+                    total = total, highlight = highlight,
+                    log_x = log_axes$log_x, log_y = log_axes$log_y,
+                    return_data = return_data))
+    }
+    lifecycle::deprecate_warn("2.6.0", "plotYield(sim2)", "plot2()",
+                              details = "Use plot2(getYield(sim1), getYield(sim2)) instead.")
+    params <- object@params
+    # We need to combine two plots
+    if (!all(dimnames(object@n)$time == dimnames(sim2@n)$time)) {
+        stop("The two simulations do not have the same times")
+    }
+    ym <- plotYield(object, species = species,
+                    tlim = tlim, total = total,
+                    log_x = log_axes$log_x, log_y = log_axes$log_y,
+                    ylim = ylim,
+                    highlight = highlight, return_data = TRUE, ...)
+    ym2 <- plotYield(sim2, species = species,
+                     tlim = tlim, total = total,
+                     log_x = log_axes$log_x, log_y = log_axes$log_y,
+                     ylim = ylim,
+                     highlight = highlight, return_data = TRUE, ...)
+    ym$Simulation <- rep(1, nrow(ym))
+    ym2$Simulation <- rep(2, nrow(ym2))
+    ym <- rbind(ym, ym2)
+
+    if (return_data) {
+        ym$Legend <- NULL
+        return(ym)
+    }
+
+    plotDataFrame(ym, params,
+                  ylab = "Yield [g/year]",
+                  xtrans = ifelse(log_axes$log_x, "log10", "identity"),
+                  ytrans = ifelse(log_axes$log_y, "log10", "identity"),
+                  ylim = ylim,
+                  highlight = highlight, wrap_var = "Simulation")
 }
 
 #' @rdname plotYield
 #' @usage NULL
 #' @export
-plotlyYield <- function(object, sim2,
+plotlyYield <- function(object, sim2 = NULL,
                         species = NULL,
                         total = FALSE,
                         log_x = FALSE, log_y = TRUE, log = NULL,
@@ -1499,39 +1474,7 @@ plotlyYield <- function(object, sim2,
                         ...) {
     argg <- as.list(environment())
     plotHover(do.call("plotYield", argg),
-             tooltip = c("Species", "Year", "Yield"))
-}
-
-# For time-series plots, keep backward compatibility with the historical
-# logical `log` argument while supporting the newer `log_x` / `log_y` form.
-#' Parse the `log` argument for time-series plots
-#'
-#' Translates the `log` argument of the time-series plotting functions into a
-#' list of logical flags for the x and y axes, falling back to `log_x` and
-#' `log_y` when `log` is `NULL`. For backwards compatibility a single logical
-#' `log` value sets only `log_y`.
-#'
-#' @param log `NULL`, a single logical value, or a character string containing
-#'   only `"x"` and/or `"y"`.
-#' @param log_x,log_y Default logical flags used when `log` is `NULL`.
-#' @return A list with logical elements `log_x` and `log_y`.
-#' @keywords internal
-parseTimePlotLog <- function(log, log_x, log_y) {
-    if (is.null(log)) {
-        return(list(log_x = log_x, log_y = log_y))
-    }
-    if (is.logical(log)) {
-        if (length(log) != 1 || is.na(log)) {
-            stop("`log` must be `NULL`, a single logical value, or a ",
-                 "character string containing only \"x\" and/or \"y\".")
-        }
-        return(list(log_x = FALSE, log_y = isTRUE(log)))
-    }
-    if (!is.character(log)) {
-        stop("`log` must be `NULL`, a single logical value, or a ",
-             "character string containing only \"x\" and/or \"y\".")
-    }
-    parsePlotLog(log, log_x = log_x, log_y = log_y)
+              tooltip = c("Species", "Year", "Yield"))
 }
 
 
@@ -1609,7 +1552,7 @@ plotYieldGear.MizerSim <- function(object,
                           ylim = c(NA, NA), tlim = c(NA, NA),
                           highlight = NULL, return_data = FALSE,
                           ...) {
-    log_axes <- parseTimePlotLog(log, log_x = log_x, log_y = log_y)
+    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     assert_that(is(object, "MizerSim"),
                 is.flag(total),
                 is.flag(log_axes$log_x),
@@ -2261,7 +2204,9 @@ plotCDF.MizerParams <- function(object, species = NULL,
 #' @keywords internal
 plot_cdf <- function(plot_dat, params, power, normalise, log_x, log_y, wlim, llim,
                      ylim, highlight, size_axis, return_data) {
+    params <- validParams(params)
     size_axis <- plot_size_axis(size_axis)
+    had_total <- "Total" %in% plot_dat$Species
     # A CDF value is cumulative *up to a size* — a boundary quantity — so it
     # belongs on the bin edges, not the bin centres. Under second-order
     # bin-averaging `plot_spectra()` returns its density on the geometric bin
@@ -2282,16 +2227,29 @@ plot_cdf <- function(plot_dat, params, power, normalise, log_x, log_y, wlim, lli
                                              drop_w = FALSE)
         plot_dat_l <- filter_plot_length_limits(plot_dat_l, llim)
         plot_dat <- plot_dat_l[, setdiff(names(plot_dat_l), "l"),
-                               drop = FALSE]
+                                drop = FALSE]
+        cdf_dat <- prepare_spectra_cdf_data(plot_dat, params,
+                                            normalise = FALSE)
+        cdf_dat <- convert_plot_size_axis(cdf_dat, params, size_axis)
+        if (had_total) {
+            cdf_dat <- add_total_line(cdf_dat, x_var = "l",
+                                      value_col = names(cdf_dat)[2])
+        }
+        if (normalise) {
+            totals <- stats::ave(cdf_dat[[names(cdf_dat)[2]]],
+                                 cdf_dat$Species, FUN = max)
+            cdf_dat[[names(cdf_dat)[2]]] <- cdf_dat[[names(cdf_dat)[2]]] / totals
+        }
+    } else {
+        cdf_dat <- prepare_spectra_cdf_data(plot_dat, params,
+                                            normalise = normalise)
+        cdf_dat <- convert_plot_size_axis(cdf_dat, params, size_axis)
     }
-    cdf_dat <- prepare_spectra_cdf_data(plot_dat, params,
-                                        normalise = normalise)
-    cdf_dat <- convert_plot_size_axis(cdf_dat, params, size_axis)
     cdf_y <- cdf_y_label(power, normalise)
     names(cdf_dat)[2] <- cdf_y
     if (return_data) return(cdf_dat)
 
-    plotDataFrame(cdf_dat, validParams(params),
+    plotDataFrame(cdf_dat, params,
                   xlab = plot_size_xlab(size_axis),
                   ylab = cdf_y,
                   xtrans = if (log_x) "log10" else "identity",
@@ -2437,14 +2395,14 @@ plotCDF2 <- function(object1, object2, name1 = "First", name2 = "Second",
     cf1 <- plotCDF(object1, species = species, wlim = wlim,
                    power = power, total = total, resource = resource,
                    background = background, normalise = normalise,
-                   size_axis = "w", return_data = TRUE, ...)
+                   size_axis = size_axis, return_data = TRUE, ...)
     cf2 <- plotCDF(object2, species = species, wlim = wlim,
                    power = power, total = total, resource = resource,
                    background = background, normalise = normalise,
-                   size_axis = "w", return_data = TRUE, ...)
+                   size_axis = size_axis, return_data = TRUE, ...)
     params <- if (is(object1, "MizerSim")) object1@params else object1
 
-    plotComparisonDataFrame(cf1, cf2, validParams(params),
+    plotComparisonDataFrame(cf1, cf2, params,
                             name1 = name1, name2 = name2,
                             xlab = plot_size_xlab(size_axis),
                             ylab = cdf_y_label(power, normalise),
@@ -2452,8 +2410,7 @@ plotCDF2 <- function(object1, object2, name1 = "First", name2 = "Second",
                             ytrans = if (log_axes$log_y) "log10" else "identity",
                             xlim = plot_size_xlim(wlim, size_axis, llim),
                             ylim = ylim, highlight = highlight,
-                            legend_var = "Legend",
-                            size_axis = size_axis)
+                            legend_var = "Legend")
 }
 
 #' Compare abundance and biomass spectra from two objects
@@ -2517,7 +2474,7 @@ plotSpectra2 <- function(object1, object2, name1 = "First", name2 = "Second",
                        return_data = TRUE, ...)
     params <- if (is(object1, "MizerSim")) object1@params else object1
 
-    plotComparisonDataFrame(sf1, sf2, validParams(params),
+    plotComparisonDataFrame(sf1, sf2, params,
                             name1 = name1, name2 = name2,
                             xlab = plot_size_xlab(size_axis),
                             ylab = spectra_y_label(
@@ -2737,8 +2694,9 @@ plotSpectraRelative <- function(object1, object2,
                        background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
     params <- if (is(object1, "MizerSim")) object1@params else object1
+    params <- validParams(params)
 
-    plotRelativeDataFrame(sf1, sf2, validParams(params),
+    plotRelativeDataFrame(sf1, sf2, params,
                           xlab = plot_size_xlab(size_axis),
                           xtrans = if (log_x) "log10" else "identity",
                           xlim = plot_size_xlim(wlim, size_axis, llim),
@@ -2936,29 +2894,24 @@ plotFeedingLevel.MizerSim <- function(object, species = NULL,
             return_data = FALSE,
             log_x = TRUE, log_y = FALSE, log = NULL,
             time_range, ...) {
-    size_axis <- plot_size_axis(size_axis)
-    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
-    assert_that(is.flag(all.sizes),
-                is.flag(include_critical),
-                length(wlim) == 2,
-                length(llim) == 2,
-                is.flag(return_data))
     if (missing(time_range)) {
         time_range  <- max(as.numeric(dimnames(object@n)$time))
     }
     params <- validParams(object@params)
     feed <- getFeedingLevel(object, time_range = time_range, drop = FALSE)
-    # If a time range was returned, average over it
     if (length(dim(feed)) == 3) {
         feed <- apply(feed, c(2, 3), mean)
     }
+    feed <- ArraySpeciesBySize(feed, value_name = "Feeding level",
+                               units = "", params = params,
+                               type = "proportion", representation = "average")
     plot_feeding_level(params, feed, species = species,
                        highlight = highlight, all.sizes = all.sizes,
                        include_critical = include_critical,
-                       log_x = log_axes$log_x, log_y = log_axes$log_y,
                        wlim = wlim, llim = llim,
                        size_axis = size_axis,
-                       return_data = return_data)
+                       return_data = return_data,
+                       log_x = log_x, log_y = log_y, log = log, ...)
 }
 
 #' @rdname plotFeedingLevel
@@ -2971,22 +2924,15 @@ plotFeedingLevel.MizerParams <- function(object, species = NULL,
             size_axis = c("w", "l"),
             return_data = FALSE,
             log_x = TRUE, log_y = FALSE, log = NULL, ...) {
-    size_axis <- plot_size_axis(size_axis)
-    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
-    assert_that(is.flag(all.sizes),
-                is.flag(include_critical),
-                length(wlim) == 2,
-                length(llim) == 2,
-                is.flag(return_data))
     params <- validParams(object)
-    feed <- getFeedingLevel(params, drop = FALSE)
+    feed <- getFeedingLevel(params)
     plot_feeding_level(params, feed, species = species,
                        highlight = highlight, all.sizes = all.sizes,
                        include_critical = include_critical,
-                       log_x = log_axes$log_x, log_y = log_axes$log_y,
                        wlim = wlim, llim = llim,
                        size_axis = size_axis,
-                       return_data = return_data)
+                       return_data = return_data,
+                       log_x = log_x, log_y = log_y, log = log, ...)
 }
 
 #' Build the feeding-level plot
@@ -3004,189 +2950,74 @@ plotFeedingLevel.MizerParams <- function(object, species = NULL,
 #'   are removed.
 #' @param include_critical Whether to also plot the critical feeding level.
 #' @param log_x,log_y Logical flags for log10 axes.
+#' @param log Optional base-R log argument string or boolean.
 #' @param wlim,llim Numeric vectors of length two giving the weight and length
 #'   limits.
 #' @param size_axis Either `"w"` (weight) or `"l"` (length).
 #' @param return_data If `TRUE`, return the plotting data frame instead of the
 #'   plot.
+#' @param ... Additional arguments passed to [plot.ArraySpeciesBySize()].
 #' @return A `mizer_plot` (ggplot2) object, or the plotting data frame if
 #'   `return_data = TRUE`.
 #' @keywords internal
 plot_feeding_level <- function(params, feed, species, highlight,
                                all.sizes, include_critical,
-                               log_x, log_y,
-                               wlim, llim, size_axis, return_data) {
+                               wlim, llim, size_axis, return_data,
+                               log_x = TRUE, log_y = FALSE, log = NULL, ...) {
     size_axis <- plot_size_axis(size_axis)
+    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
+    log_x <- log_axes$log_x
+    log_y <- log_axes$log_y
 
-    # selector for desired species
-    sel_sp <- valid_species_arg(params, species, return.logical = TRUE,
-                                error_on_empty = TRUE)
-    species <- dimnames(params@initial_n)$sp[sel_sp]
-    feed <- feed[sel_sp, , drop = FALSE]
-
-    plot_dat <- data.frame(w = rep(params@w, each = length(species)),
-                           value = c(feed),
-                           Species = species)
-
-    if (include_critical) {
-        feed_crit <- getCriticalFeedingLevel(params)[sel_sp, , drop = FALSE]
-        plot_dat_crit <- data.frame(
-            w = rep(params@w, each = length(species)),
-            value = c(feed_crit),
-            Species = species)
-        plot_dat$Type <- "actual"
-        plot_dat_crit$Type <- "critical"
-        plot_dat <- rbind(plot_dat, plot_dat_crit)
-    }
-
-    if (!all.sizes) {
-        # Remove feeding level for sizes outside a species' size range
-        for (sp in species) {
-            plot_dat$value[plot_dat$Species == sp &
-                               (plot_dat$w < params@species_params[sp, "w_min"] |
-                                    plot_dat$w > params@species_params[sp, "w_max"])] <- NA
+    if (!include_critical) {
+        if (return_data) {
+            return(plot(feed, species = species, all.sizes = all.sizes,
+                        highlight = highlight, wlim = wlim, llim = llim,
+                        size_axis = size_axis, return_data = TRUE,
+                        log_x = log_x, log_y = log_y, ...))
         }
-        plot_dat <- plot_dat[complete.cases(plot_dat), ]
-    }
-    if (!is.na(wlim[1])) plot_dat <- plot_dat[plot_dat$w >= wlim[1], ]
-    if (!is.na(wlim[2])) plot_dat <- plot_dat[plot_dat$w <= wlim[2], ]
-    plot_dat <- convert_plot_size_axis(plot_dat, params, size_axis)
-    if (identical(size_axis, "l")) {
-        plot_dat <- filter_plot_length_limits(plot_dat, llim)
-    }
-    if (log_y) {
-        # Remove non-positive values because log scales cannot represent them.
-        plot_dat <- subset(plot_dat, value > 0)
-    }
-
-    names(plot_dat)[2] <- "Feeding level"
-    if (return_data) return(plot_dat)
-    x_var <- plot_size_x_var(size_axis)
-
-    # Need to keep species in order for legend
-    legend_levels <-
-        intersect(names(params@linecolour), plot_dat$Species)
-    plot_dat$Legend <- factor(plot_dat$Species, levels = legend_levels)
-    linesize <- make_linesize(legend_levels, highlight)
-
-    # The feeding level array declares itself a proportion, so show the whole
-    # of [0, 1] — widened if the critical feeding level rises above it, which
-    # it can. A logarithmic axis is left to the data.
-    feeding_level_ylim <- if (log_y) NULL else
-        array_ylim(feed, c(NA, NA), log_y, plot_dat[["Feeding level"]])
-
-    # We do not use `plotDataFrame()` to create the plot because it would not
-    # handle the alpha transparency for the critical feeding level.
-
-    # The reason why below `group = species` is included in `ggplot()`
-    # rather than in `geom_line` is because that puts it first in the
-    # plotly tooltips, due to a bug in plotly.
-    if (include_critical) {
-        plot_dat$Species <- interaction(plot_dat$Species, plot_dat$Type)
-        p <- ggplot(plot_dat, aes(group = Species,
-                                  alpha = Type)) +
-            scale_discrete_manual("alpha", name = "Feeding Level",
-                                  values = c(actual = 1, critical = 0.5))
-    } else {
-        p <- ggplot(plot_dat, aes(group = Species))
+        p <- plot(feed, species = species, all.sizes = all.sizes,
+                  highlight = highlight, wlim = wlim, llim = llim,
+                  size_axis = size_axis, log_x = log_x, log_y = log_y, ...)
+        if (!log_y) {
+            p <- p + coord_cartesian(ylim = c(0, 1))
+        }
+        return(p)
+    }
+    crit <- getCriticalFeedingLevel(params)
+    df_actual <- plot(feed, species = species, all.sizes = all.sizes,
+                      wlim = wlim, llim = llim, size_axis = size_axis,
+                      log_x = log_x, log_y = log_y,
+                      return_data = TRUE, ...)
+    df_crit <- plot(crit, species = species, all.sizes = all.sizes,
+                    wlim = wlim, llim = llim, size_axis = size_axis,
+                    log_x = log_x, log_y = log_y,
+                    return_data = TRUE, ...)
+    names(df_crit)[2] <- names(df_actual)[2]
+    df_actual$Type <- "actual"
+    df_crit$Type <- "critical"
+    if (return_data) {
+        return(rbind(df_actual, df_crit))
+    }
+    p <- plot(feed, species = species, all.sizes = all.sizes,
+              highlight = highlight, wlim = wlim, llim = llim,
+              size_axis = size_axis, log_x = log_x, log_y = log_y, ...)
+    p <- addPlot(p, crit, species = species, all.sizes = all.sizes,
+                 wlim = wlim, llim = llim, size_axis = size_axis,
+                 alpha = 0.5)
+    if (!log_y) {
+        feeding_level_ylim <- array_ylim(feed, c(NA, NA), log_y,
+                                         c(df_actual[[2]], df_crit[[2]]))
+        p <- p + coord_cartesian(ylim = feeding_level_ylim)
     }
-    make_mizer_plot(
-        p + geom_line(aes(x = .data[[x_var]], y = .data[["Feeding level"]],
-                          colour = Legend, linetype = Legend, linewidth = Legend)) +
-            scale_x_continuous(name = plot_size_xlab(size_axis),
-                               trans = if (log_x) "log10" else "identity",
-                               limits = plot_size_xlim(wlim, size_axis, llim)) +
-            scale_y_continuous(name = "Feeding Level",
-                               trans = if (log_y) "log10" else "identity") +
-            coord_cartesian(ylim = feeding_level_ylim) +
-            scale_colour_manual(values = params@linecolour[legend_levels]) +
-            scale_linetype_manual(values = params@linetype[legend_levels]) +
-            scale_discrete_manual("linewidth", values = linesize),
-        plot_size_tooltip(size_axis, before = "Species", after = "Feeding level")
-    )
+    p
 }
 
 #' @rdname plotFeedingLevel
 #' @usage NULL
 #' @export
-plotlyFeedingLevel <- function(object,
-                             species = NULL,
-                             time_range,
-                             all.sizes = FALSE,
-                             highlight = NULL,
-                             include_critical = FALSE,
-                             wlim = c(NA, NA), llim = c(NA, NA),
-                             size_axis = c("w", "l"),
-                             log_x = TRUE, log_y = FALSE, log = NULL, ...) {
-    size_axis <- plot_size_axis(size_axis)
-    argg <- as.list(environment())
-    p <- plotHover(do.call("plotFeedingLevel", argg))
-
-    # When critical feeding level is included, ggplotly creates traces split by the
-    # interaction of Species and Type, which produces a very long combined legend.
-    # The code below reshapes the legend to mirror the ggplot output:
-    # - Species appear once under a "Legend" group
-    # - A separate "Feeding Level" group shows "actual" and "critical"
-    if (isTRUE(include_critical)) {
-        species_seen <- character(0)
-        for (i in seq_along(p$x$data)) {
-            tr <- p$x$data[[i]]
-            # Only adjust line traces with names like "(actual, Cod)"
-            if (!identical(tr$type, "scatter") ||
-                is.null(tr$mode) || !grepl("lines", tr$mode)) next
-            nm <- tr$name
-            if (!is.null(nm) && grepl("^\\(", nm)) {
-                nm_clean <- gsub("^\\(|\\)$", "", nm)
-                parts <- strsplit(nm_clean, ",\\s*")[[1]]
-                if (length(parts) >= 2) {
-                    typ <- parts[[1]]
-                    sp <- paste(parts[-1], collapse = ",")
-                    # Rename the trace to species name and group under "Legend"
-                    p$x$data[[i]]$name <- sp
-                    p$x$data[[i]]$legendgroup <- "Legend"
-                    # Add group title once
-                    if (length(species_seen) == 0) {
-                        p$x$data[[i]]$legendgrouptitle <- list(text = "Legend")
-                    }
-                    # Hide duplicate legend entries for "critical" traces
-                    if (identical(typ, "critical")) {
-                        p$x$data[[i]]$showlegend <- FALSE
-                    }
-                    species_seen <- unique(c(species_seen, sp))
-                }
-            }
-        }
-
-        # Add two legend-only traces to show the "Feeding Level" group
-        p <- plotly::add_trace(
-            p,
-            x = c(0, 1), y = c(0, 1),
-            type = "scatter", mode = "lines",
-            name = "actual",
-            legendgroup = "Feeding Level",
-            legendgrouptitle = list(text = "Feeding Level"),
-            line = list(color = "blue"),
-            opacity = 1,
-            hoverinfo = "skip",
-            visible = "legendonly",
-            showlegend = TRUE,
-            inherit = FALSE
-        )
-        p <- plotly::add_trace(
-            p,
-            x = c(0, 1), y = c(0, 1),
-            type = "scatter", mode = "lines",
-            name = "critical",
-            legendgroup = "Feeding Level",
-            line = list(color = "blue"),
-            opacity = 0.5,
-            hoverinfo = "skip",
-            visible = "legendonly",
-            showlegend = TRUE,
-            inherit = FALSE
-        )
-    }
-    p
+plotlyFeedingLevel <- function(object, ...) {
+    plotHover(plotFeedingLevel(object, ...))
 }
 
 #' Plot predation mortality rate of each species against size
diff --git a/inst/skills/analyse-and-plot/SKILL.md b/inst/skills/analyse-and-plot/SKILL.md
index 73bee0d23..c1e3c8755 100644
--- a/inst/skills/analyse-and-plot/SKILL.md
+++ b/inst/skills/analyse-and-plot/SKILL.md
@@ -249,10 +249,9 @@ Which arguments apply depends on the array's shape:
   `ArraySpeciesBySize` method, so it accepts everything that method does plus
   `time` (default: the last time step). It has no `tlim`: only one time is
   shown.
-- `plot(<ArrayResourceBySize>)` accepts `log_x`, `log_y`, `wlim`, `ylim` only.
-  The resource is a single spectrum, so there is nothing for `species`,
-  `highlight`, `total` or `background` to select, and no `size_axis`/`llim`
-  because the resource has no length-weight relationship.
+- `plot(<ArrayResourceBySize>)` accepts `log_x`, `log_y`, `wlim`, `llim`, `ylim`,
+  `size_axis`, `per_log_size`. The resource is a single spectrum, so there is
+  nothing for `species`, `highlight`, `total` or `background` to select.
 - `plot(<ArrayTimeByResourceBySize>)` accepts the same as
   `ArrayResourceBySize` plus `time`.
 
@@ -315,9 +314,10 @@ plot(initialN(params), size_axis = "l", per_log_size = TRUE) # per log length
 plot(initialNResource(params), per_log_size = TRUE)          # resource too
 ```
 
-`per_log_size` needs no weight-length relationship, only `size_axis` does, so
-unlike a length axis it is available for the resource classes. Asking for it on
-an array that does not hold a density is an error rather than being ignored.
+Both `size_axis` and `per_log_size` are available for the resource classes as
+well, using the resource's weight-length relationship from `resource_params()`.
+Asking for `per_log_size` on an array that does not hold a density is an error
+rather than being ignored.
 
 ### Plotting proportions
 
diff --git a/man/parseTimePlotLog.Rd b/man/parseTimePlotLog.Rd
deleted file mode 100644
index 24483a7c4..000000000
--- a/man/parseTimePlotLog.Rd
+++ /dev/null
@@ -1,24 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/plots.R
-\name{parseTimePlotLog}
-\alias{parseTimePlotLog}
-\title{Parse the \code{log} argument for time-series plots}
-\usage{
-parseTimePlotLog(log, log_x, log_y)
-}
-\arguments{
-\item{log}{\code{NULL}, a single logical value, or a character string containing
-only \code{"x"} and/or \code{"y"}.}
-
-\item{log_x, log_y}{Default logical flags used when \code{log} is \code{NULL}.}
-}
-\value{
-A list with logical elements \code{log_x} and \code{log_y}.
-}
-\description{
-Translates the \code{log} argument of the time-series plotting functions into a
-list of logical flags for the x and y axes, falling back to \code{log_x} and
-\code{log_y} when \code{log} is \code{NULL}. For backwards compatibility a single logical
-\code{log} value sets only \code{log_y}.
-}
-\keyword{internal}
diff --git a/man/plot_feeding_level.Rd b/man/plot_feeding_level.Rd
index 5b4e8f5b7..7991878a0 100644
--- a/man/plot_feeding_level.Rd
+++ b/man/plot_feeding_level.Rd
@@ -11,12 +11,14 @@ plot_feeding_level(
   highlight,
   all.sizes,
   include_critical,
-  log_x,
-  log_y,
   wlim,
   llim,
   size_axis,
-  return_data
+  return_data,
+  log_x = TRUE,
+  log_y = FALSE,
+  log = NULL,
+  ...
 )
 }
 \arguments{
@@ -33,8 +35,6 @@ are removed.}
 
 \item{include_critical}{Whether to also plot the critical feeding level.}
 
-\item{log_x, log_y}{Logical flags for log10 axes.}
-
 \item{wlim, llim}{Numeric vectors of length two giving the weight and length
 limits.}
 
@@ -42,6 +42,12 @@ limits.}
 
 \item{return_data}{If \code{TRUE}, return the plotting data frame instead of the
 plot.}
+
+\item{log_x, log_y}{Logical flags for log10 axes.}
+
+\item{log}{Optional base-R log argument string or boolean.}
+
+\item{...}{Additional arguments passed to \code{\link[=plot.ArraySpeciesBySize]{plot.ArraySpeciesBySize()}}.}
 }
 \value{
 A \code{mizer_plot} (ggplot2) object, or the plotting data frame if
diff --git a/tests/testthat/test-ArrayResourceBySize-class.R b/tests/testthat/test-ArrayResourceBySize-class.R
index 2dd215e3f..5ca6cda55 100644
--- a/tests/testthat/test-ArrayResourceBySize-class.R
+++ b/tests/testthat/test-ArrayResourceBySize-class.R
@@ -292,3 +292,45 @@ test_that("the resource joins the species on a length-axis spectrum", {
     rp <- resource_params(params)
     expect_equal(res_l$l, (res_w$w / rp$a)^(1 / rp$b))
 })
+
+test_that("plot2, plotRelative and addPlot support length axis and per_log_size for resource arrays", {
+    params <- setResource(NS_params_small)
+    r1 <- initialNResource(params)
+    r2 <- r1
+    r2[] <- r1 * 1.5
+
+    # ArrayResourceBySize methods
+    p2_l <- plot2(r1, r2, size_axis = "l")
+    expect_s3_class(p2_l, "ggplot")
+    expect_identical(p2_l$scales$get_scales("x")$name, "Length [cm]")
+    expect_identical(p2_l$scales$get_scales("y")$name, "Number density [1/cm]")
+
+    p2_log <- plot2(r1, r2, per_log_size = TRUE)
+    expect_s3_class(p2_log, "ggplot")
+    expect_identical(p2_log$scales$get_scales("y")$name, "Number density in log weight")
+
+    pr_l <- plotRelative(r1, r2, size_axis = "l")
+    expect_s3_class(pr_l, "ggplot")
+    expect_identical(pr_l$scales$get_scales("x")$name, "Length [cm]")
+
+    pa_l <- addPlot(plot(r1, size_axis = "l"), r2, size_axis = "l")
+    expect_s3_class(pa_l, "ggplot")
+    expect_equal(length(pa_l$layers), 2)
+
+    # ArrayTimeByResourceBySize methods
+    nr1 <- NResource(NS_sim_small)
+    nr2 <- nr1
+    nr2[] <- nr1 * 1.2
+    t2_l <- plot2(nr1, nr2, size_axis = "l")
+    expect_s3_class(t2_l, "ggplot")
+    expect_identical(t2_l$scales$get_scales("x")$name, "Length [cm]")
+
+    tr_l <- plotRelative(nr1, nr2, size_axis = "l")
+    expect_s3_class(tr_l, "ggplot")
+    expect_identical(tr_l$scales$get_scales("x")$name, "Length [cm]")
+
+    ta_l <- addPlot(plot(nr1, size_axis = "l"), nr2, size_axis = "l")
+    expect_s3_class(ta_l, "ggplot")
+    expect_equal(length(ta_l$layers), 2)
+})
+
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index 803c3720b..0993d6e62 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -998,8 +998,8 @@ test_that("yield plotting helpers validate comparison and gear selection", {
 })
 
 test_that("yield plotly wrappers return plotly objects", {
-    expect_s3_class(plotlyYield(sim, species = species, ylim = c(1e-5, 1)), "plotly")
-    expect_s3_class(plotlyYieldGear(sim, species = species, ylim = c(1e-5, 1)), "plotly")
+    expect_s3_class(plotlyYield(sim, species = species, ylim = c(1e5, 1e12)), "plotly")
+    expect_s3_class(plotlyYieldGear(sim, species = species, ylim = c(1e5, 1e12)), "plotly")
 })
 
 test_that("yield plotting helpers accept ylim", {
@@ -1041,7 +1041,7 @@ test_that("return_data is identical",{
                                  return_data = TRUE)), c(22, 4))
 
     expect_equal(dim(plotFeedingLevel(sim, species = species,
-                                      return_data = TRUE)), c(40, 3))
+                                      return_data = TRUE)), c(40, 4))
 
     expect_equal(dim(plotPredMort(sim, species = species,
                                   return_data = TRUE)), c(40, 4))
@@ -1515,3 +1515,60 @@ test_that("plotBiomass validates time range and can include total", {
     d <- plotBiomass(sim, species = "Cod", total = TRUE, return_data = TRUE)
     expect_true("Total" %in% d$Species)
 })
+
+test_that("plotCDF and comparisons work with length axis and total", {
+    params <- NS_params_small
+    sim <- NS_sim_small
+
+    # plotCDF on length axis with total
+    cdf_l <- plotCDF(params, size_axis = "l", total = TRUE, return_data = TRUE)
+    expect_identical(names(cdf_l)[[1]], "l")
+    expect_true("Total" %in% cdf_l$Species)
+    expect_equal(max(cdf_l[["Cumulative proportion of biomass"]][cdf_l$Species == "Total"]), 1)
+
+    # plotCDF2 with size_axis = "l" and total = TRUE
+    p2 <- plotCDF2(params, params, size_axis = "l", total = TRUE)
+    expect_s3_class(p2, "ggplot")
+    expect_identical(p2$scales$get_scales("x")$name, "Length [cm]")
+    expect_true("Total" %in% levels(p2$data$Legend))
+
+    # plotSpectra2 with size_axis = "l" and total = TRUE
+    s2 <- plotSpectra2(params, params, size_axis = "l", total = TRUE)
+    expect_s3_class(s2, "ggplot")
+    expect_identical(s2$scales$get_scales("x")$name, "Length [cm]")
+    expect_true("Total" %in% levels(s2$data$Legend))
+
+    # plotSpectraRelative with size_axis = "l" and total = TRUE
+    sr <- plotSpectraRelative(params, params, size_axis = "l", total = TRUE)
+    expect_s3_class(sr, "ggplot")
+    expect_identical(sr$scales$get_scales("x")$name, "Length [cm]")
+    expect_true("Total" %in% levels(sr$data$Legend))
+})
+
+test_that("plotFeedingLevel and plotlyFeedingLevel with include_critical", {
+    params <- NS_params_small
+
+    p <- plotFeedingLevel(params, include_critical = TRUE)
+    expect_s3_class(p, "ggplot")
+    expect_equal(length(p$layers), 2)
+
+    df <- plotFeedingLevel(params, include_critical = TRUE, return_data = TRUE)
+    expect_true(all(c("actual", "critical") %in% df$Type))
+
+    pl <- plotlyFeedingLevel(params, include_critical = TRUE)
+    expect_s3_class(pl, "plotly")
+})
+
+test_that("plotYield delegates to plot(getYield()) and warns on sim2", {
+    sim <- NS_sim_small
+
+    # Single sim plot
+    p <- plotYield(sim)
+    expect_s3_class(p, "ggplot")
+    expect_identical(p$scales$get_scales("x")$name, "Year")
+    expect_identical(p$scales$get_scales("y")$name, "Yield [g/year]")
+
+    # sim2 deprecation warning
+    lifecycle::expect_deprecated(plotYield(sim, sim))
+})
+
diff --git a/vignettes/cheatsheet-analysis-and-plotting.Rmd b/vignettes/cheatsheet-analysis-and-plotting.Rmd
index 2d3ac8744..5a5e81c79 100644
--- a/vignettes/cheatsheet-analysis-and-plotting.Rmd
+++ b/vignettes/cheatsheet-analysis-and-plotting.Rmd
@@ -262,10 +262,9 @@ Which arguments apply depends on the array's shape:
   `ArraySpeciesBySize` method, so it accepts everything that method does plus
   `time` (default: the last time step). It has no `tlim`: only one time is
   shown.
-- `plot(<ArrayResourceBySize>)` accepts `log_x`, `log_y`, `wlim`, `ylim` only.
-  The resource is a single spectrum, so there is nothing for `species`,
-  `highlight`, `total` or `background` to select, and no `size_axis`/`llim`
-  because the resource has no length-weight relationship.
+- `plot(<ArrayResourceBySize>)` accepts `log_x`, `log_y`, `wlim`, `llim`, `ylim`,
+  `size_axis`, `per_log_size`. The resource is a single spectrum, so there is
+  nothing for `species`, `highlight`, `total` or `background` to select.
 - `plot(<ArrayTimeByResourceBySize>)` accepts the same as
   `ArrayResourceBySize` plus `time`.
 
@@ -328,9 +327,10 @@ plot(initialN(params), size_axis = "l", per_log_size = TRUE) # per log length
 plot(initialNResource(params), per_log_size = TRUE)          # resource too
 ```
 
-`per_log_size` needs no weight-length relationship, only `size_axis` does, so
-unlike a length axis it is available for the resource classes. Asking for it on
-an array that does not hold a density is an error rather than being ignored.
+Both `size_axis` and `per_log_size` are available for the resource classes as
+well, using the resource's weight-length relationship from [`resource_params()`](../reference/resource_params.html).
+Asking for `per_log_size` on an array that does not hold a density is an error
+rather than being ignored.
 
 ### Plotting proportions
 
@@ -372,7 +372,7 @@ plotSpectra(sim, biomass = TRUE, size_axis = "l")     # biomass density against
 **The resource has its own length convention.** It is a composite of many taxa,
 so instead of a taxonomic weight-length relationship it uses the equivalent
 spherical diameter of an organism with the density of water (`a = pi/6`,
-`b = 3`, in [`resource_params()`](../reference/resource_params.html)). It therefore appears on a length axis, but
+`b = 3`, in `resource_params()`). It therefore appears on a length axis, but
 measured differently from the fish: a fish of a given weight is about 3.7 times
 longer than a sphere of that weight. That gap at the resource-consumer boundary
 is real biology, not an artefact.