diff --git a/NAMESPACE b/NAMESPACE
index fcc5e98..6e61043 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -16,6 +16,8 @@ export(silv_biomass)
 export(silv_density_hart)
 export(silv_density_ntrees_ha)
 export(silv_density_sdi)
+export(silv_density_sdi_auto)
+export(silv_density_sdi_class)
 export(silv_diametric_class)
 export(silv_dominant_height)
 export(silv_lorey_height)
diff --git a/R/data.R b/R/data.R
index b739df3..f8098c8 100644
--- a/R/data.R
+++ b/R/data.R
@@ -107,3 +107,21 @@
 #' MITECO. 4th Spanish National Forest Inventory - SIG database codes.
 #' \url{https://www.miteco.gob.es/content/dam/miteco/es/biodiversidad/temas/inventarios-nacionales/documentador_sig_tcm30-536622.pdf}
 "snfi4_volume_coefficients"
+
+
+#' SDI beta coefficients
+#'
+#' Specific beta coefficients for Reineke's Stand Density Index (SDI) per
+#' species, country, and region.
+#'
+#' @format A `tibble`
+#' \describe{
+#'   \item{article_id}{Character. Short identifier of the source article.}
+#'   \item{title}{Character. Full title of the source article.}
+#'   \item{doi_url}{Character. DOI URL of the source article.}
+#'   \item{country}{Character. Country where the study was conducted.}
+#'   \item{region}{Character. Region within the country.}
+#'   \item{species}{Character. Scientific name of the tree species.}
+#'   \item{beta}{Numeric. Beta coefficient for SDI calculation.}
+#' }
+"sdi_coefficients"
diff --git a/R/metrics-stand-density.R b/R/metrics-stand-density.R
index 6af2fb1..28e3856 100644
--- a/R/metrics-stand-density.R
+++ b/R/metrics-stand-density.R
@@ -1,5 +1,3 @@
-
-
 #' Calculates number of trees per hectare
 #'
 #' Calculates number of trees per hectare for a given plot size and shape
@@ -28,12 +26,11 @@
 #'       n,
 #'       plot_size = c(10, 15),
 #'       plot_shape = "rectangular"
-#'      )
+#'     )
 #'   )
 silv_density_ntrees_ha <- function(ntrees,
-                           plot_size,
-                           plot_shape = "circular") {
-
+                                   plot_size,
+                                   plot_shape = "circular") {
   # 0. Handle errors
   stopifnot(plot_shape %in% c("circular", "rectangular"))
   if (length(plot_size) == 1 && plot_size <= 0) cli::cli_abort("`plot_size` has to be greater than 0")
@@ -44,8 +41,6 @@ silv_density_ntrees_ha <- function(ntrees,
   } else {
     ntrees * 10000 / prod(plot_size)
   }
-
-
 }
 
 
@@ -54,72 +49,417 @@ silv_density_ntrees_ha <- function(ntrees,
 
 #' Calculates the Stand Density Index
 #'
-#' The Stand Density Index (SDI) is relationship between the average tree size and 
+#' The Stand Density Index (SDI) is the relationship between the average tree size and
 #' density of trees per hectare.
 #'
 #' @template ntrees
 #' @template dg
-#' @param classify whether to classify the values using USDA thresholds
-#' @param max_sdi used when \code{classify = TRUE}. The maximum SDi, which depends
-#' on the species, stand type, and site
-#' 
-#' @return A numeric vector
+#' @param beta The Stand Density Index exponent (default is \code{1.605}).
+#'
+#' @return A numeric vector representing the absolute SDI.
 #' @export
-#' 
+#'
 #' @details
-#' The SDI has different interpretation depending on the species, location, and also
+#' The SDI has different interpretations depending on the species, location, and also
 #' the management type (even-aged, uneven-aged...). The value of maximum SDI must
-#' be determined from the literature and used carefully. The option \code{classify = TRUE}
-#' will use this value to classify the SDI in low density (<24%), moderate density (24-35%),
-#' high density (34-55%), and extremely high density (>55%).
+#' be determined from the literature and used carefully. The \code{beta} exponent allows
+#' adjustments for different species or mixed stands.
+#'
+#' @references Reineke, L. H. (1933). Perfecting a stand-density index for even-aged forests.
+#'   Journal of Agricultural Research, 46(7), 627-638. URL: https://research.fs.usda.gov/download/treesearch/60134.pdf
 #'
 #' @examples
-#' ## calculate SDI for a Pinus sulvestris stand (max 990)
-#' silv_density_sdi(ntrees = 800, dg = 23.4, max_sdi = 990)
-#' 
-#' ## check base classification (other can be used)
-#' silv_density_sdi(ntrees = 800, dg = 23.4, classify = TRUE, max_sdi = 990)
+#' ## calculate SDI for a Pinus sylvestris stand (beta = 1.605)
+#' silv_density_sdi(ntrees = 800, dg = 23.4)
+#'
+#' ## calculate SDI with custom beta
+#' silv_density_sdi(ntrees = 800, dg = 23.4, beta = 1.7)
 silv_density_sdi <- function(
-  ntrees, 
-  dg, 
-  classify = FALSE, 
-  max_sdi = NULL
+  ntrees,
+  dg,
+  beta = 1.605
 ) {
-
-  # 0. Validate inputs
+  # 0. validate inputs
   assert_positive_numeric(ntrees, "ntrees")
   assert_positive_numeric(dg, "dg")
-  assert_logical(classify, "classify")
+  if (!is.numeric(beta)) cli::cli_abort("{.arg beta} has to be a numeric vector.")
   assert_same_length(ntrees, dg, names = c("ntrees", "dg"))
 
+  # 1. calculate sdi
+  sdi <- ntrees * ((25.4 / dg)**-abs(beta)) # note: abs() avoids errors with signs
+  return(sdi)
+}
 
-  # 1. Calculate SDI
-  sdi <- ntrees * ((dg / 25.4)) ** 1.605
+# --- Internal Auto Selector Helper ---
 
-  # 2. Classify?
-  if (classify) {
+#' @noRd
+.auto_select_sdi_beta <- function(species, country = NULL, region = NULL) {
+  
+  sdi_coefficients <- silviculture::sdi_coefficients
+  
+  # 1. Try exact species + country + region match
+  if (!is.null(country) && !is.null(region)) {
+    sel <- sdi_coefficients[sdi_coefficients$species == species & 
+                            sdi_coefficients$country == country & 
+                            sdi_coefficients$region == region, ]
+    if (nrow(sel) > 0) {
+      return(list(
+        model = sel$article_id[1], 
+        beta = sel$beta[1], 
+        matched_species = species, 
+        matched_country = country,
+        matched_region = region,
+        is_fallback = FALSE,
+        fallback_type = NA_character_,
+        model_desc = paste0(sel$article_id[1], " (", country, ", ", region, ")")
+      ))
+    }
+  }
+  
+  # 2. Try species + country + "all" regions match
+  if (!is.null(country)) {
+    sel <- sdi_coefficients[sdi_coefficients$species == species & 
+                            sdi_coefficients$country == country & 
+                            sdi_coefficients$region == "all", ]
+    if (nrow(sel) > 0) {
+      return(list(
+        model = sel$article_id[1], 
+        beta = sel$beta[1], 
+        matched_species = species, 
+        matched_country = country,
+        matched_region = "all",
+        is_fallback = !is.null(region),
+        fallback_type = "region",
+        model_desc = paste0(sel$article_id[1], " (", country, ", all regions)")
+      ))
+    }
+  }
 
-    ## assert inputs
-    if (is.null(max_sdi)) cli::cli_abort("You must specify <max_sdi> when <classify = TRUE>")
-    assert_positive_numeric(max_sdi, "max_sdi")
-
-    ## calculate
-    sdi <- (sdi / max_sdi) * 100
-    sdi <- dplyr::case_when(
-      sdi <= 24            ~ "Low density",
-      sdi > 24 & sdi <= 34 ~ "Moderate density",
-      sdi > 34 & sdi <= 55 ~ "High density",
-      sdi > 55             ~ "Extremely high density" 
-    )
-  } else if (!is.null(max_sdi)) {
-    sdi <- (sdi / max_sdi) * 100
+  # 3. Fallback when country is not found/specified, but we find the species in some other country
+  sel <- sdi_coefficients[sdi_coefficients$species == species, ]
+  if (nrow(sel) > 0) {
+    if (!is.null(region)) {
+      sel_reg <- sel[sel$region == region, ]
+      if (nrow(sel_reg) > 0) {
+        return(list(
+          model = sel_reg$article_id[1], 
+          beta = sel_reg$beta[1], 
+          matched_species = species, 
+          matched_country = sel_reg$country[1],
+          matched_region = region,
+          is_fallback = TRUE,
+          fallback_type = "country",
+          model_desc = paste0(sel_reg$article_id[1], " (", sel_reg$country[1], ", ", region, ")")
+        ))
+      }
+    }
+    
+    sel_all <- sel[sel$region == "all", ]
+    if (nrow(sel_all) > 0) {
+      return(list(
+        model = sel_all$article_id[1], 
+        beta = sel_all$beta[1], 
+        matched_species = species, 
+        matched_country = sel_all$country[1],
+        matched_region = "all",
+        is_fallback = TRUE,
+        fallback_type = "region",
+        model_desc = paste0(sel_all$article_id[1], " (", sel_all$country[1], ", all regions)")
+      ))
+    }
+    
+    return(list(
+      model = sel$article_id[1], 
+      beta = sel$beta[1], 
+      matched_species = species, 
+      matched_country = sel$country[1],
+      matched_region = sel$region[1],
+      is_fallback = TRUE,
+      fallback_type = "region",
+      model_desc = paste0(sel$article_id[1], " (", sel$country[1], ", ", sel$region[1], ")")
+    ))
   }
-  return(sdi)
+  
+  # 4. Try genus fallback (genus spp.)
+  genus <- strsplit(species, " ")[[1]][1]
+  genus_spp <- paste0(genus, " spp.")
+  
+  if (!is.null(country)) {
+    if (!is.null(region)) {
+      sel <- sdi_coefficients[sdi_coefficients$species == genus_spp & 
+                              sdi_coefficients$country == country & 
+                              sdi_coefficients$region == region, ]
+      if (nrow(sel) > 0) {
+        return(list(
+          model = sel$article_id[1], 
+          beta = sel$beta[1], 
+          matched_species = genus_spp, 
+          matched_country = country,
+          matched_region = region,
+          is_fallback = TRUE,
+          fallback_type = "genus",
+          model_desc = paste0(sel$article_id[1], " (genus fallback: ", country, ", ", region, ")")
+        ))
+      }
+    }
+    
+    sel <- sdi_coefficients[sdi_coefficients$species == genus_spp & 
+                            sdi_coefficients$country == country & 
+                            sdi_coefficients$region == "all", ]
+    if (nrow(sel) > 0) {
+      return(list(
+        model = sel$article_id[1], 
+        beta = sel$beta[1], 
+        matched_species = genus_spp, 
+        matched_country = country,
+        matched_region = "all",
+        is_fallback = TRUE,
+        fallback_type = "genus",
+        model_desc = paste0(sel$article_id[1], " (genus fallback: ", country, ", all regions)")
+      ))
+    }
+  }
+  
+  sel <- sdi_coefficients[sdi_coefficients$species == genus_spp, ]
+  if (nrow(sel) > 0) {
+    return(list(
+      model = sel$article_id[1], 
+      beta = sel$beta[1], 
+      matched_species = genus_spp, 
+      matched_country = sel$country[1],
+      matched_region = sel$region[1],
+      is_fallback = TRUE,
+      fallback_type = "genus",
+      model_desc = paste0(sel$article_id[1], " (genus fallback: ", sel$country[1], ", ", sel$region[1], ")")
+    ))
+  }
+  
+  # 5. Total fallback (default)
+  sel <- sdi_coefficients[sdi_coefficients$species == "default", ]
+  if (nrow(sel) == 0) {
+    default_beta <- -1.605
+    default_desc <- "reineke-1933 (-1.605)"
+  } else {
+    default_beta <- sel$beta[1]
+    default_desc <- paste0(sel$article_id[1], " (", default_beta, ")")
+  }
+  
+  return(list(
+    model = sel$article_id[1], 
+    beta = default_beta, 
+    matched_species = "default", 
+    matched_country = "default",
+    matched_region = "default",
+    is_fallback = TRUE,
+    fallback_type = "default",
+    model_desc = default_desc
+  ))
+}
 
+#' Predict Stand Density Index automatically
+#'
+#' @description
+#' `silv_density_sdi_auto()` is a vectorized function that automatically selects
+#' the best available Stand Density Index exponent (\code{beta}) for each row
+#' based on a provided species, country, and region from the internal \code{sdi_coefficients} database.
+#' 
+#' If an exact species, country, and region match is not found, the function falls back to a 
+#' country-wide species model (\code{region = "all"}), then searches other countries, then falls
+#' back to a genus-level fallback (e.g., "Pinus spp."), and finally to the default SDI exponent 
+#' (\code{beta = -1.605} from Reineke 1933).
+#'
+#' @template ntrees
+#' @template dg
+#' @param species A character string or vector of tree species (e.g., `"Pinus sylvestris"`).
+#' @param country A character string or vector of the country (e.g., `"Spain"`).
+#'   Defaults to `NULL` (no country specified).
+#' @param region A character string or vector of the region (e.g., `"Castilla y León"`).
+#'   Defaults to `NULL` (no region specified).
+#' @param quiet Logical. If `FALSE`, informs the user about fallbacks to genus or default models.
+#'
+#' @return A `data.frame` with three columns:
+#'   - `sdi`: The computed absolute Stand Density Index.
+#'   - `beta`: The beta exponent used for the calculation.
+#'   - `sdi_model`: The model used (e.g., `"del-rio-2006 (Spain, Castilla y León)"`, 
+#'     `"reineke-1933 (-1.605)"`, etc.).
+#'
+#' @name silv_density_sdi_auto
+#'
+#' @examples
+#' # Calculate SDI with automatic selection
+#' silv_density_sdi_auto(
+#'   ntrees = 800,
+#'   dg = 23.4,
+#'   species = "Pinus sylvestris",
+#'   region = "Castilla y León"
+#' )
+#'
+#' # Fallback to default
+#' silv_density_sdi_auto(
+#'   ntrees = 800,
+#'   dg = 23.4,
+#'   species = "Unknown species"
+#' )
+#'
+#' @export
+silv_density_sdi_auto <- function(
+  ntrees,
+  dg,
+  species,
+  country = NULL,
+  region = NULL,
+  quiet = FALSE
+) {
+  # Validations
+  n_trees <- length(ntrees)
+  assert_positive_numeric(ntrees, "ntrees")
+  assert_positive_numeric(dg, "dg")
+  assert_same_length(ntrees, dg, names = c("ntrees", "dg"))
+
+  if (length(species) == 1) {
+    species <- rep(species, n_trees)
+  } else if (length(species) != n_trees) {
+    cli::cli_abort("{.arg species} must be of length 1 or the same length as {.arg ntrees}.")
+  }
+
+  if (is.null(country)) {
+    country <- rep(NA_character_, n_trees)
+  } else if (length(country) == 1) {
+    country <- rep(country, n_trees)
+  } else if (length(country) != n_trees) {
+    cli::cli_abort("{.arg country} must be of length 1 or the same length as {.arg ntrees}.")
+  }
+
+  if (is.null(region)) {
+    region <- rep(NA_character_, n_trees)
+  } else if (length(region) == 1) {
+    region <- rep(region, n_trees)
+  } else if (length(region) != n_trees) {
+    cli::cli_abort("{.arg region} must be of length 1 or the same length as {.arg ntrees}.")
+  }
+
+  sdi_values <- rep(NA_real_, n_trees)
+  beta_values <- rep(NA_real_, n_trees)
+  sdi_models_used <- rep(NA_character_, n_trees)
+  
+  unique_combos <- unique(data.frame(
+    species = species, 
+    country = country,
+    region = region, 
+    stringsAsFactors = FALSE
+  ))
+  
+  for (i in seq_len(nrow(unique_combos))) {
+    sp <- unique_combos$species[i]
+    cnt <- unique_combos$country[i]
+    reg <- unique_combos$region[i]
+    
+    idx_sp <- species == sp
+    
+    if (is.na(cnt)) {
+      idx_cnt <- is.na(country)
+      cnt_arg <- NULL
+    } else {
+      idx_cnt <- country == cnt
+      cnt_arg <- cnt
+    }
+    
+    if (is.na(reg)) {
+      idx_reg <- is.na(region)
+      reg_arg <- NULL
+    } else {
+      idx_reg <- region == reg
+      reg_arg <- reg
+    }
+    
+    idx <- which(idx_sp & idx_cnt & idx_reg)
+    
+    best_model_info <- .auto_select_sdi_beta(sp, cnt_arg, reg_arg)
+    
+    if (best_model_info$is_fallback && !quiet) {
+      if (best_model_info$fallback_type == "default") {
+         cli::cli_alert_info("Exact model for {.val {sp}} not found. Using default beta {.val {best_model_info$beta}}.")
+      } else if (best_model_info$fallback_type == "genus") {
+         cli::cli_alert_info("Exact species {.val {sp}} not found. Using genus fallback {.val {best_model_info$matched_species}} from {.val {best_model_info$model}}.")
+      } else if (best_model_info$fallback_type == "region") {
+         cli::cli_alert_info("Exact region {.val {reg_arg}} not found for {.val {sp}}. Using fallback region {.val {best_model_info$matched_region}} from {.val {best_model_info$model}}.")
+      } else if (best_model_info$fallback_type == "country") {
+         cli::cli_alert_info("Exact country {.val {cnt_arg}} not found for {.val {sp}}. Using fallback country {.val {best_model_info$matched_country}} from {.val {best_model_info$model}}.")
+      }
+    }
+    
+    sdi_values[idx] <- silv_density_sdi(ntrees[idx], dg[idx], beta = best_model_info$beta)
+    beta_values[idx] <- best_model_info$beta
+    sdi_models_used[idx] <- best_model_info$model_desc
+  }
+  
+  return(data.frame(
+    sdi = sdi_values,
+    beta = beta_values,
+    sdi_model = sdi_models_used,
+    stringsAsFactors = FALSE
+  ))
 }
 
 
+#' Classifies the Stand Density Index
+#'
+#' Classifies the Stand Density Index (SDI) into density classes or calculates the relative SDI
+#' percentage based on USDA thresholds.
+#'
+#' @param sdi A numeric vector representing the Stand Density Index.
+#' @param max_sdi A numeric vector representing the maximum SDI for the species/site.
+#' @param classify A logical value indicating whether to classify the values into density classes
+#'   (default is \code{TRUE}). If \code{FALSE}, it returns the relative SDI as a percentage.
+#'
+#' @return A character vector with the density classes if \code{classify = TRUE}, or a numeric vector
+#'   with the relative SDI percentage if \code{classify = FALSE}.
+#' @export
+#'
+#' @details
+#' The option \code{classify = TRUE} will use the \code{max_sdi} value to classify the SDI into
+#' four competitive and growth conditions: low density (<24%), moderate density (24-35%),
+#' high density (34-55%), and extremely high density (>55%).
+#'
+#' @references USDA Forest Service. (n.d.). Stand Density Index.
+#'   https://www.fs.usda.gov/Internet/FSE_DOCUMENTS/stelprdb5270993.pdf
+#'
+#' @examples
+#' ## calculate SDI for a Pinus sylvestris stand (max 990)
+#' sdi_val <- silv_density_sdi(ntrees = 800, dg = 23.4)
+#'
+#' ## check base classification
+#' silv_density_sdi_class(sdi = sdi_val, max_sdi = 990)
+#'
+#' ## get relative SDI percentage
+#' silv_density_sdi_class(sdi = sdi_val, max_sdi = 990, classify = FALSE)
+silv_density_sdi_class <- function(
+  sdi,
+  max_sdi,
+  classify = TRUE
+) {
+  # 0. validate inputs
+  assert_positive_numeric(sdi, "sdi")
+  assert_positive_numeric(max_sdi, "max_sdi")
+  assert_logical(classify, "classify")
+  assert_same_length(sdi, max_sdi, names = c("sdi", "max_sdi"))
+
+  # 1. calculate relative sdi
+  rel_sdi <- (sdi / max_sdi) * 100
 
+  # 2. classify or return percentage
+  if (classify) {
+    res <- dplyr::case_when(
+      rel_sdi <= 24 ~ "Low density",
+      rel_sdi > 24 & rel_sdi <= 34 ~ "Moderate density",
+      rel_sdi > 34 & rel_sdi <= 55 ~ "High density",
+      rel_sdi > 55 ~ "Extremely high density"
+    )
+  } else {
+    res <- rel_sdi
+  }
+  return(res)
+}
 
 
 #' Hart or Hart-Becking spacing index
@@ -163,7 +503,6 @@ silv_density_hart <- function(
   ntrees,
   which = c("hart", "hart-becking")
 ) {
-
   # 0. Validate inputs
   assert_positive_numeric(h0, "h0")
   assert_positive_numeric(ntrees, "ntrees")
@@ -176,9 +515,4 @@ silv_density_hart <- function(
     "hart-becking" = sqrt(20000 / (ntrees * sqrt(3))) / h0 * 100,
     cli::cli_abort("`which` must be either <hart> or <hart-becking>")
   )
-
 }
-
-
-
-
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 2b7b908..635f8ca 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -77,6 +77,8 @@ reference:
     - silv_density_ntrees_ha
     - silv_density_hart
     - silv_density_sdi
+    - silv_density_sdi_auto
+    - silv_density_sdi_class
     - silv_spacing_index
     - silv_ntrees_ha
 
@@ -103,6 +105,7 @@ reference:
     contents:
     - biomass_models
     - carbon_models
+    - sdi_coefficients
     - snfi3_volume_coefficients
     - snfi4_volume_coefficients
     - inventory_samples
diff --git a/data/sdi_coefficients.rda b/data/sdi_coefficients.rda
new file mode 100644
index 0000000..a8b8fb2
Binary files /dev/null and b/data/sdi_coefficients.rda differ
diff --git a/man/sdi_coefficients.Rd b/man/sdi_coefficients.Rd
new file mode 100644
index 0000000..1f5b7cb
--- /dev/null
+++ b/man/sdi_coefficients.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{sdi_coefficients}
+\alias{sdi_coefficients}
+\title{SDI beta coefficients}
+\format{
+A \code{tibble}
+\describe{
+\item{article_id}{Character. Short identifier of the source article.}
+\item{title}{Character. Full title of the source article.}
+\item{doi_url}{Character. DOI URL of the source article.}
+\item{country}{Character. Country where the study was conducted.}
+\item{region}{Character. Region within the country.}
+\item{species}{Character. Scientific name of the tree species.}
+\item{beta}{Numeric. Beta coefficient for SDI calculation.}
+}
+}
+\usage{
+sdi_coefficients
+}
+\description{
+Specific beta coefficients for Reineke's Stand Density Index (SDI) per
+species, country, and region.
+}
+\keyword{datasets}
diff --git a/man/silv_density_ntrees_ha.Rd b/man/silv_density_ntrees_ha.Rd
index a9361c9..35dfd4b 100644
--- a/man/silv_density_ntrees_ha.Rd
+++ b/man/silv_density_ntrees_ha.Rd
@@ -37,6 +37,6 @@ inventory_samples |>
       n,
       plot_size = c(10, 15),
       plot_shape = "rectangular"
-     )
+    )
   )
 }
diff --git a/man/silv_density_sdi.Rd b/man/silv_density_sdi.Rd
index 422bc24..6800fb2 100644
--- a/man/silv_density_sdi.Rd
+++ b/man/silv_density_sdi.Rd
@@ -4,7 +4,7 @@
 \alias{silv_density_sdi}
 \title{Calculates the Stand Density Index}
 \usage{
-silv_density_sdi(ntrees, dg, classify = FALSE, max_sdi = NULL)
+silv_density_sdi(ntrees, dg, beta = 1.605)
 }
 \arguments{
 \item{ntrees}{Numeric vector with number of trees of the diameter class per
@@ -13,29 +13,29 @@ corresponds to only one tree}
 
 \item{dg}{Numeric vector of quadratic mean diameters}
 
-\item{classify}{whether to classify the values using USDA thresholds}
-
-\item{max_sdi}{used when \code{classify = TRUE}. The maximum SDi, which depends
-on the species, stand type, and site}
+\item{beta}{The Stand Density Index exponent (default is \code{1.605}).}
 }
 \value{
-A numeric vector
+A numeric vector representing the absolute SDI.
 }
 \description{
-The Stand Density Index (SDI) is relationship between the average tree size and
+The Stand Density Index (SDI) is the relationship between the average tree size and
 density of trees per hectare.
 }
 \details{
-The SDI has different interpretation depending on the species, location, and also
+The SDI has different interpretations depending on the species, location, and also
 the management type (even-aged, uneven-aged...). The value of maximum SDI must
-be determined from the literature and used carefully. The option \code{classify = TRUE}
-will use this value to classify the SDI in low density (<24\%), moderate density (24-35\%),
-high density (34-55\%), and extremely high density (>55\%).
+be determined from the literature and used carefully. The \code{beta} exponent allows
+adjustments for different species or mixed stands.
 }
 \examples{
-## calculate SDI for a Pinus sulvestris stand (max 990)
-silv_density_sdi(ntrees = 800, dg = 23.4, max_sdi = 990)
+## calculate SDI for a Pinus sylvestris stand (beta = 1.605)
+silv_density_sdi(ntrees = 800, dg = 23.4)
 
-## check base classification (other can be used)
-silv_density_sdi(ntrees = 800, dg = 23.4, classify = TRUE, max_sdi = 990)
+## calculate SDI with custom beta
+silv_density_sdi(ntrees = 800, dg = 23.4, beta = 1.7)
+}
+\references{
+Reineke, L. H. (1933). Perfecting a stand-density index for even-aged forests.
+Journal of Agricultural Research, 46(7), 627-638. URL: https://research.fs.usda.gov/download/treesearch/60134.pdf
 }
diff --git a/man/silv_density_sdi_auto.Rd b/man/silv_density_sdi_auto.Rd
new file mode 100644
index 0000000..1c08cb7
--- /dev/null
+++ b/man/silv_density_sdi_auto.Rd
@@ -0,0 +1,68 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/metrics-stand-density.R
+\name{silv_density_sdi_auto}
+\alias{silv_density_sdi_auto}
+\title{Predict Stand Density Index automatically}
+\usage{
+silv_density_sdi_auto(
+  ntrees,
+  dg,
+  species,
+  country = NULL,
+  region = NULL,
+  quiet = FALSE
+)
+}
+\arguments{
+\item{ntrees}{Numeric vector with number of trees of the diameter class per
+hectare. If \code{ntrees = NULL}, the function will assume that each diameter
+corresponds to only one tree}
+
+\item{dg}{Numeric vector of quadratic mean diameters}
+
+\item{species}{A character string or vector of tree species (e.g., \code{"Pinus sylvestris"}).}
+
+\item{country}{A character string or vector of the country (e.g., \code{"Spain"}).
+Defaults to \code{NULL} (no country specified).}
+
+\item{region}{A character string or vector of the region (e.g., \code{"Castilla y León"}).
+Defaults to \code{NULL} (no region specified).}
+
+\item{quiet}{Logical. If \code{FALSE}, informs the user about fallbacks to genus or default models.}
+}
+\value{
+A \code{data.frame} with three columns:
+\itemize{
+\item \code{sdi}: The computed absolute Stand Density Index.
+\item \code{beta}: The beta exponent used for the calculation.
+\item \code{sdi_model}: The model used (e.g., \code{"del-rio-2006 (Spain, Castilla y León)"},
+\code{"reineke-1933 (-1.605)"}, etc.).
+}
+}
+\description{
+\code{silv_density_sdi_auto()} is a vectorized function that automatically selects
+the best available Stand Density Index exponent (\code{beta}) for each row
+based on a provided species, country, and region from the internal \code{sdi_coefficients} database.
+
+If an exact species, country, and region match is not found, the function falls back to a
+country-wide species model (\code{region = "all"}), then searches other countries, then falls
+back to a genus-level fallback (e.g., "Pinus spp."), and finally to the default SDI exponent
+(\code{beta = -1.605} from Reineke 1933).
+}
+\examples{
+# Calculate SDI with automatic selection
+silv_density_sdi_auto(
+  ntrees = 800,
+  dg = 23.4,
+  species = "Pinus sylvestris",
+  region = "Castilla y León"
+)
+
+# Fallback to default
+silv_density_sdi_auto(
+  ntrees = 800,
+  dg = 23.4,
+  species = "Unknown species"
+)
+
+}
diff --git a/man/silv_density_sdi_class.Rd b/man/silv_density_sdi_class.Rd
new file mode 100644
index 0000000..f92a509
--- /dev/null
+++ b/man/silv_density_sdi_class.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/metrics-stand-density.R
+\name{silv_density_sdi_class}
+\alias{silv_density_sdi_class}
+\title{Classifies the Stand Density Index}
+\usage{
+silv_density_sdi_class(sdi, max_sdi, classify = TRUE)
+}
+\arguments{
+\item{sdi}{A numeric vector representing the Stand Density Index.}
+
+\item{max_sdi}{A numeric vector representing the maximum SDI for the species/site.}
+
+\item{classify}{A logical value indicating whether to classify the values into density classes
+(default is \code{TRUE}). If \code{FALSE}, it returns the relative SDI as a percentage.}
+}
+\value{
+A character vector with the density classes if \code{classify = TRUE}, or a numeric vector
+with the relative SDI percentage if \code{classify = FALSE}.
+}
+\description{
+Classifies the Stand Density Index (SDI) into density classes or calculates the relative SDI
+percentage based on USDA thresholds.
+}
+\details{
+The option \code{classify = TRUE} will use the \code{max_sdi} value to classify the SDI into
+four competitive and growth conditions: low density (<24\%), moderate density (24-35\%),
+high density (34-55\%), and extremely high density (>55\%).
+}
+\examples{
+## calculate SDI for a Pinus sylvestris stand (max 990)
+sdi_val <- silv_density_sdi(ntrees = 800, dg = 23.4)
+
+## check base classification
+silv_density_sdi_class(sdi = sdi_val, max_sdi = 990)
+
+## get relative SDI percentage
+silv_density_sdi_class(sdi = sdi_val, max_sdi = 990, classify = FALSE)
+}
+\references{
+USDA Forest Service. (n.d.). Stand Density Index.
+https://www.fs.usda.gov/Internet/FSE_DOCUMENTS/stelprdb5270993.pdf
+}
diff --git a/tests/testthat/test-stand-density.R b/tests/testthat/test-stand-density.R
index c54fd2c..ba02281 100644
--- a/tests/testthat/test-stand-density.R
+++ b/tests/testthat/test-stand-density.R
@@ -77,5 +77,86 @@ test_that("Errors work", {
   )
   expect_error(silv_density_hart("17.8", 400))
   expect_error(silv_density_hart(17.8, "400"))
-  expect_error(silv_density_hart(c(17.8, 20.5), 400))
+  expect_error(
+    silv_density_hart(c(17.8, 20.5), 400)
+  )
+})
+
+
+# 3. silv_density_sdi ----------------------------------------------------
+
+## Tests
+test_that("Stand Density Index is well calculated", {
+  # default beta (1.605)
+  expect_equal(
+    silv_density_sdi(ntrees = 800, dg = 23.4),
+    702.40,
+    tolerance = 0.01
+  )
+
+  # custom beta
+  expect_equal(
+    silv_density_sdi(ntrees = 800, dg = 23.4, beta = 1.7),
+    692.68,
+    tolerance = 0.01
+  )
+
+  # negative beta
+  expect_equal(
+    silv_density_sdi(ntrees = 800, dg = 23.4, beta = -1.605),
+    702.40,
+    tolerance = 0.01
+  )
+
+  # with max_sdi (returns percentage) using silv_density_sdi_class
+  sdi_val <- silv_density_sdi(ntrees = 800, dg = 23.4)
+  expect_equal(
+    silv_density_sdi_class(sdi = sdi_val, max_sdi = 990, classify = FALSE),
+    70.95,
+    tolerance = 0.01
+  )
+
+  # with classification using silv_density_sdi_class
+  expect_equal(
+    silv_density_sdi_class(sdi = sdi_val, max_sdi = 990),
+    "Extremely high density"
+  )
+})
+
+## Test errors
+test_that("Errors work in silv_density_sdi", {
+  expect_error(silv_density_sdi(800, 23.4, beta = "1.605"))
+})
+
+test_that("Errors work in silv_density_sdi_class", {
+  expect_error(silv_density_sdi_class(700, "990"))
+  expect_error(silv_density_sdi_class(700, 990, classify = "TRUE"))
+})
+
+
+# 4. silv_density_sdi_auto -----------------------------------------------
+
+test_that("silv_density_sdi_auto calculates and falls back correctly", {
+  
+  # Exact match
+  res1 <- silv_density_sdi_auto(800, 23.4, "Pinus sylvestris", country = "Spain", region = "Castilla y León", quiet = TRUE)
+  expect_equal(res1$sdi_model, "del-rio-2006 (Spain, Castilla y León)")
+  expect_equal(res1$beta, -1.75)
+  expect_true(is.numeric(res1$sdi))
+  
+  # Region fallback ("all")
+  res2 <- silv_density_sdi_auto(800, 23.4, "Pinus pinaster", country = "Spain", region = "Unknown region", quiet = TRUE)
+  expect_equal(res2$sdi_model, "aguirre-2017 (Spain, all regions)")
+  expect_equal(res2$beta, -1.9477)
+  
+  # Default fallback
+  res3 <- silv_density_sdi_auto(800, 23.4, "Unknown species", quiet = TRUE)
+  expect_equal(res3$sdi_model, "reineke-1933 (-1.605)")
+  expect_equal(res3$beta, -1.605)
+  
+})
+
+test_that("Errors work in silv_density_sdi_auto", {
+  expect_error(silv_density_sdi_auto(800, 23.4, species = 123))
+  expect_error(silv_density_sdi_auto(c(800, 700), c(23.4, 25.1), species = c("Pinus sylvestris", "Pinus pinaster", "Quercus robur")))
 })