diff --git a/DESCRIPTION b/DESCRIPTION
index 7fdadd0c4d72f5780370a15e20025bc2284ff4f2..81120399d7984b88c00c3a81c8ffb33624612b25 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -12,13 +12,14 @@ Authors@R: c(
   person("Marcos", "Raül", , "raul.marcos@bsc.es", role = "ctb"),
   person("Palma", "Lluis", , "lluis.palma@bsc.es", role = "ctb"),
   person("An-Chi", "Ho", , "an.ho@bsc.es", role = c("ctb")),
+  person("Javier", "Corvillo", , "javier.corvillo@bsc.es", role = c("ctb")),
   person("BSC-CNS", role = "cph"))
 Description: Set of generalised tools for the flexible computation of climate 
   related indicators defined by the user. Each method represents a specific 
   mathematical approach which is combined with the possibility to select an 
   arbitrary time period to define the indicator. This enables a wide range of 
   possibilities to tailor the most suitable indicator for each particular climate 
-  service application (agriculture, food security, energy, water management, ...). 
+  service application (agriculture, food security, energy, water management, health...). 
   This package is intended for sub-seasonal, seasonal and decadal climate 
   predictions, but its methods are also applicable to other time-scales, 
   provided the dimensional structure of the input is maintained. Additionally, 
diff --git a/R/RNoughtIndices.R b/R/RNoughtIndices.R
new file mode 100644
index 0000000000000000000000000000000000000000..204408de17c7c3647f6a78f5fd2ceea753a76a07
--- /dev/null
+++ b/R/RNoughtIndices.R
@@ -0,0 +1,528 @@
+#'R0 Index computation on s2dv_cube objects
+#'
+#'@author Javier Corvillo, \email{javier.corvillo@bsc.es}
+#'@description This function computes the environmental contribution to Aedes-borne
+#' disease transmissibility. Values higher than 1 imply that the environmental conditions allow for 
+#' the disease to be transmitted to more than 1 person (the infection can spread) while R0 values lower
+#'  than 1 imply that the environmental conditions are not suitable for the disease to spread.
+#' 
+#' This function utilizes four possible ento-epidemiological models, each of them stated
+#' in Caminade el al., 2015; Liu-Helmerssohn et al., 2014; and Mordecai et
+#' al., 2017 and Wesolowski et al., 2015 respectively. Additionally, an adjustment to
+#' real life data, using recorded DENV data recorded between 2014-2017 in the Caribbean, can also be performed.
+#' 
+#'@references Caminade, Cyril et al. (Jan. 2017). ‘Global risk model for vector-borne 
+#' transmission of Zika virus reveals the role of El Niño 2015’.
+#' In: Proceedings of the National Academy of Sciences 114.1.
+#' Publisher: Proceedings of the National Academy of 
+#' Sciences, pp. 119–124. doi: 10.1073/pnas.1614303114. 
+#' url: https://www.pnas.org/doi/full/10.1073/pnas.1614303114.
+#' 
+#'@references Liu-Helmersson, Jing et al. (Mar. 2014). ‘Vectorial Capacity of Aedes
+#' aegypti: Effects of Temperature and Implications for Global Dengue Epidemic Potential’. en.
+#' In: PLOS ONE 9.3.
+#' Publisher: Public Library of Science, e89783.
+#' issn: 1932-6203. doi: 10.1371/journal.pone.0089783.
+#' url: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0089783.
+#' 
+#' @references Mordecai, Erin A. et al. (Apr. 2017). ‘Detecting the impact of temperature on
+#' transmission of Zika, dengue, and chikungunya using mechanistic models’. en. 
+#' In: PLOS Neglected Tropical Diseases 11.4.
+#' Publisher: Public Library of Science, e0005568.
+#' issn: 1935-2735. doi: 10.1371/journal.pntd.0005568.
+#' url: https://journals.plos.org/plosntds/article?id=10.1371/journal.pntd.0005568.
+#' 
+#' @references Wesolowski, Amy et al. (Sept. 2015). ‘Impact of human mobility on the emergence
+#' of dengue epidemics in Pakistan’. 
+#' In: Proceedings of the National Academy of Sciences 112.38.
+#' Publisher: Proceedings of the National Academy of Sciences, pp. 11887–11892.
+#' doi: 10.1073/pnas.1504964112.
+#' url: https://www.pnas.org/doi/full/10.1073/pnas.1504964112.
+#'
+#'@param temp An s2dv_cube object with temperature values, expressed in degrees Celsius. If
+#' "caminade" is selected in "method", named spatial dimensions are required in the s2dv_cube object, 
+#' ordered as ("latitude", "longitude"). Valid dimension names are "lon", "longitude", "lat", "latitude".
+#' 
+#' @param method a string indicating the R0 ento-epidemiological
+#' model used to obtain the R0 outputs. Methods include "caminade",
+#' "wesolowski", "liuhelmersson" and "mordecai". Alternatively, an 
+#' extrapolation to real life data can be made by setting "method" to "empirical".
+#' 
+#' @param mm A Kramer probabilty matrix of dimensions (2, longitude, latitude) to be passed
+#' if "caminade" is selected as the computation method. The spatial dimensions must
+#' match those of temp and be equal in length. Valid dimension names are "lon", "longitude",
+#' "lat", "latitude". 
+#' 
+#' @param lon_dim a character vector indicating the longitude dimension name in
+#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
+#'
+#' @param lat_dim a character vector indicating the latitude dimension name in
+#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
+#' 
+#'@param ncores An integer indicating the number of cores to use in parallel
+#'  computation for temporal subsetting.
+#'@return An s2dv_cube object containing the R0 Indices (unitless).
+#'
+#'@examples
+
+# dims <- c(time = 100, lat = 5, lon = 4)
+# temp_cube <- list(
+#   data = array(
+#     runif(prod(dims), min = 15, max = 35),
+#     dim = dims
+#   ),
+#   attrs = list(
+#     Variable = list(varName = '2m Temperature')
+#   )
+# )
+  
+# class(temp_cube) <- 's2dv_cube'
+
+# R0 <- CST_RNoughtIndices(temp_cube, method = "mordecai", mm = NULL, lon_dim = NULL, lat_dim = NULL, ncores = NULL)
+
+# Caminade method (requires Kramer matrix)
+# mm <- array(runif(2 * 5 * 4, 0.1, 1.0), 
+#             dim = c(2, 4, 5))
+# names(dim(mm)) <- c("", "lon", "lat")
+# R0 <- CST_RNoughtIndices(temp_cube, method = "caminade", mm = mm, lon_dim = "lon", lat_dim = "lat", ncores = NULL)
+
+#'@import multiApply
+#'@export
+
+CST_RNoughtIndices <- function(
+  temp,
+  method,
+  mm,
+  lon_dim = lon_dim,
+  lat_dim = lat_dim,
+  ncores = NULL) {
+    
+  # Check "s2dv_cube"
+  if (!inherits(temp, "s2dv_cube")) {
+    stop("Parameter 'temp' must be of the class 's2dv_cube'.")
+  } 
+  
+  RNought <- RNoughtIndices(
+    temp = temp$data,
+    method = method,
+    mm = mm,
+    lon_dim = lon_dim,
+    lat_dim = lat_dim,
+    ncores = ncores
+  )
+  
+  temp$data <- RNought
+
+  if ("Variable" %in% names(temp$attrs)) {
+    if ("varName" %in% names(temp$attrs$Variable)) {
+      temp$attrs$Variable$varName <- "RNoughtIndex"
+    }
+  }
+
+  return(temp)
+}
+
+#'R0 Index computation on s2dv_cube objects
+#'
+#'@author Javier Corvillo, \email{javier.corvillo@bsc.es}
+#'@description This function computes the environmental contribution to Aedes-borne
+#' disease transmissibility. Values higher than 1 imply that the environmental conditions allow for 
+#' the disease to be transmitted to more than 1 person (the infection can spread) while R0 values lower
+#'  than 1 imply that the environmental conditions are not suitable for the disease to spread.
+#' 
+#' This function utilizes four possible ento-epidemiological models, each of them stated
+#' in Caminade el al., 2015; Liu-Helmerssohn et al., 2014; and Mordecai et
+#' al., 2017 and Wesolowski et al., 2015 respectively. Additionally, an adjustment to
+#' real life data, using recorded DENV data recorded between 2014-2017 in the Caribbean, can also be performed.
+#' 
+#'@references Caminade, Cyril et al. (Jan. 2017). ‘Global risk model for vector-borne 
+#' transmission of Zika virus reveals the role of El Niño 2015’.
+#' In: Proceedings of the National Academy of Sciences 114.1.
+#' Publisher: Proceedings of the National Academy of 
+#' Sciences, pp. 119–124. doi: 10.1073/pnas.1614303114. 
+#' url: https://www.pnas.org/doi/full/10.1073/pnas.1614303114.
+#' 
+#'@references Liu-Helmersson, Jing et al. (Mar. 2014). ‘Vectorial Capacity of Aedes
+#' aegypti: Effects of Temperature and Implications for Global Dengue Epidemic Potential’. en.
+#' In: PLOS ONE 9.3.
+#' Publisher: Public Library of Science, e89783.
+#' issn: 1932-6203. doi: 10.1371/journal.pone.0089783.
+#' url: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0089783.
+#' 
+#' @references Mordecai, Erin A. et al. (Apr. 2017). ‘Detecting the impact of temperature on
+#' transmission of Zika, dengue, and chikungunya using mechanistic models’. en. 
+#' In: PLOS Neglected Tropical Diseases 11.4.
+#' Publisher: Public Library of Science, e0005568.
+#' issn: 1935-2735. doi: 10.1371/journal.pntd.0005568.
+#' url: https://journals.plos.org/plosntds/article?id=10.1371/journal.pntd.0005568.
+#' 
+#' @references Wesolowski, Amy et al. (Sept. 2015). ‘Impact of human mobility on the emergence
+#' of dengue epidemics in Pakistan’. 
+#' In: Proceedings of the National Academy of Sciences 112.38.
+#' Publisher: Proceedings of the National Academy of Sciences, pp. 11887–11892.
+#' doi: 10.1073/pnas.1504964112.
+#' url: https://www.pnas.org/doi/full/10.1073/pnas.1504964112.
+#'
+#'@param temp An array with temperature values, expressed in degrees Celsius. If
+#' "caminade" is selected in "method", named spatial dimensions are required in the s2dv_cube object, 
+#' ordered as ("latitude", "longitude"). Valid dimension names are "lon", "longitude", "lat", "latitude".
+#' 
+#' @param method a string indicating the R0 ento-epidemiological
+#' model used to obtain the R0 outputs. Methods include "caminade",
+#' "wesolowski", "liuhelmersson" and "mordecai". Alternatively, an 
+#' extrapolation to real life data can be made by setting "method" to "empirical".
+#' 
+#' @param mm A Kramer probabilty matrix of dimensions (2, longitude, latitude) to be passed
+#' if "caminade" is selected as the computation method. The spatial dimensions must
+#' match those of temp and be equal in length. Valid dimension names are "lon", "longitude",
+#' "lat", "latitude". 
+#' 
+#' @param lon_dim a character vector indicating the longitude dimension name in
+#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
+#'
+#' @param lat_dim a character vector indicating the latitude dimension name in
+#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
+#' 
+#'@param ncores An integer indicating the number of cores to use in parallel
+#'  computation for temporal subsetting.
+#'@return An s2dv_cube object containing the R0 Indices (unitless).
+#'
+#'@examples
+
+# dims <- c(time = 100, lat = 5, lon = 4)
+# temp <- array(
+#     runif(prod(dims), min = 15, max = 35),
+#     dim = dims
+#   )
+  
+# R0 <- RNoughtIndices(temp, method = "mordecai", mm = NULL, lon_dim = NULL, lat_dim = NULL, ncores = NULL)
+
+# Caminade method (requires Kramer matrix)
+# mm <- array(runif(2 * 5 * 4, 0.1, 1.0), 
+#             dim = c(2, 4, 5))
+# names(dim(mm)) <- c("", "lon", "lat")
+# R0 <- RNoughtIndices(temp, method = "caminade", mm = mm, lon_dim = "lon", lat_dim = "lat", ncores = NULL)
+
+#'@import multiApply
+#'@export
+
+RNoughtIndices <-  function(
+  temp,
+  method,
+  mm,
+  lon_dim, 
+  lat_dim,
+  ncores = NULL) {
+
+  if (!(method %in% c("caminade", "wesolowski", "liuhelmerssohn", "mordecai",
+  "empirical"))) {
+    stop("Invalid value for method. Accepted methods include 'caminade',
+    'wesolowski', 'liuhelmerssohn' and 'mordecai'. Alternatively, an empirical
+    calculation can be performed with 'empirical'")
+  
+  } else if (method == "caminade") {
+
+    if (is.null(mm)) {
+      stop("The method 'caminade' has been selected but no Kramer probability matrix", 
+      " in 'mm' has been provided.")
+    }
+
+    # Ensure the input probability matrix has dimensions (2, lon_dim, lat_dim)
+    if (length(dim(mm)) != 3 || dim(mm)[1] != 2) {
+      stop("The probability matrix 'mm' must be a 3D array with dimensions (2, lon_dim, lat_dim).")
+    }
+    
+    # Ensure the input temperature array is at least 2D with lon_dim and lat_dim
+    if (!is.array(temp) || length(dim(temp)) < 2) {
+      stop("The temperature array must be at least 2D with spatial dimensions (lat_dim, lon_dim).")
+    }
+    
+    if (is.null(names(dim(mm))) || is.null(names(dim(temp))) || 
+        !all(c(lon_dim, lat_dim) %in% names(dim(mm))) ||
+        !all(c(lon_dim, lat_dim) %in% names(dim(temp)))) {
+      stop("Both 'mm' and 'temperature' must have the spatial dimensions set in 'lon_dim'",
+      " and 'lat_dim'.")
+    }
+
+    if (dim(mm)[which(names(dim(mm)) == lat_dim)] != dim(temp)[which(names(dim(temp)) == lat_dim)] ||
+      dim(mm)[which(names(dim(mm)) == lon_dim)] != dim(temp)[which(names(dim(temp)) == lon_dim)]) {
+      stop("The arrays 'mm' and 'temperature' don't have matching spatial dimension lengths")
+    }
+
+    env_index <- Apply(
+      data = array(temp, c(index = 1, dim(temp))),
+      target_dims = c(lat_dim, lon_dim), 
+      fun = .caminade, 
+      mm = mm, 
+      nlat = dim(mm)[which(names(dim(mm)) == lat_dim)], 
+      nlon = dim(mm)[which(names(dim(mm)) == lon_dim)], 
+      ncores = ncores
+    )$output1
+
+    env_index <- drop(aperm(env_index, c(3:length(dim(env_index)), 1, 2)))
+  
+  } else if (method == "liuhelmerssohn") {
+      
+    env_index <- Apply(
+      data = array(temp, c(index = 1, dim(temp))),
+      target_dims = c("index"),
+      fun = .liuhelmerssohn,
+      ncores = ncores
+    )$output1
+
+    env_index <- drop(array(env_index, dim = dim(env_index)[-1]))
+
+  } else if (method == "mordecai") {
+      
+    env_index <- Apply(
+      data = array(temp, c(index = 1, dim(temp))),
+      target_dims = c("index"),
+      fun = .mordecai,
+      ncores = ncores
+    )$output1
+
+    env_index <- drop(array(env_index, dim = dim(env_index)[-1]))
+
+  } else if (method == "wesolowski") {
+      
+    env_index <- Apply(
+      data = array(temp, c(index = 1, dim(temp))),
+      target_dims = c("index"),
+      fun = .wesolowski,
+      ncores = ncores
+    )$output1
+
+    env_index <- drop(array(env_index, dim = dim(env_index)[-1]))
+
+  } else if (method == "empirical") {
+
+    empirical <- read.table(
+      system.file(
+        "index_empirical", "empirical.txt", package = "CSIndicators", mustWork = TRUE
+      )
+    )
+
+    env_index <- Apply(
+      data = array(temp, c(index = 1, dim(temp))),
+      target_dims = c("index"),
+      fun = .empirical,
+      empirical = empirical,
+      ncores = ncores
+    )$output1
+      
+  }
+  return(env_index)
+}
+
+.caminade <- function(temp, mm, nlat, nlon) {
+
+  # Setting parameter values for Caminade et al., 2015:
+  a <- array(NA, dim = c(2, nlat, nlon))
+  a[1, , ] <- 0.0043 * temp + 0.0943  #biting rates per day
+  a[2, , ] <- 0.5 * (0.0043 * temp + 0.0943)
+
+  # Vector preference (0-1); typical: 0.88-1.0
+  phi <- NULL
+  phi[1] <- 1
+  #typical: 0.24-1.0
+  phi[2] <- 0.5
+
+  # Transmission probability (vector to host); typical: 0.1-0.75
+  b <- NULL
+  b[1] <- 0.5
+  b[2] <- 0.5
+
+  # Transmission probability (host to vector); typical: 0.1-0.75
+  beta <- NULL
+  beta[1] <- 0.1
+  beta[2] <- 0.033
+  mu <- array(NA, dim = c(2, nlat, nlon))
+
+
+  for (ix in 1:nlon) {
+    for (iy in 1:nlat) {
+      # Including option for NA
+      if (is.na(temp[iy, ix]) == TRUE) {
+          mu[1, iy, ix] <- NA
+          mu[2, iy, ix] <- NA
+      } else {
+        if (temp[iy, ix] < 22) {
+          mu[1, iy, ix] <- 1 / (1.22 + exp(-3.05 + 0.72 * temp[iy, ix])) + 0.196
+        } else {
+          mu[1, iy, ix] <- 1 / (1.24 + exp(35.14 - 0.905 * temp[iy, ix])) + 0.195
+        }
+        if (temp[iy, ix] < 15) {
+          mu[2, iy, ix] <- 1 / (1.1 + exp(-4.04 + 0.576 * temp[iy, ix])) + 0.12
+        } else if (15 <= temp[iy, ix] && temp[iy, ix] < 26.3) {
+          mu[2, iy, ix] <- 0.000339 * temp[iy, ix]^2 - 0.0189 * temp[iy, ix] + 0.336
+        } else if (temp[iy, ix] >= 26.3) {
+          mu[2, iy, ix] <- 1 / (1.065 + exp(32.2 - 0.92 * temp[iy, ix])) + 0.0747
+        }
+      }
+    }
+  }
+          
+  # Inverse of the EIP
+  eip <- array(NA, dim = c(2, nlat, nlon))
+  eip[1, , ] <- 4 + exp(5.15 - 0.123 * temp)
+  eip[2, , ] <- 1.03 * (4 + exp(5.15 - 0.123 * temp))
+
+  #Recovery rate
+  r <- 1 / 7
+
+  # Calculating the R0 parameters:
+  env_index <- array(NA, dim = c(2, nlat, nlon))
+  for (i in 1:2) {
+    for (ix in 1:nlon) {
+      for (iy in 1:nlat) {
+        if (is.na(temp[iy, ix]) == TRUE) {
+          env_index[i, iy, ix] <- NA
+        } else {
+          env_index[i, iy, ix] <- (b[i] * beta[i] * a[i, iy, ix]^2 / mu[i, iy, ix]) *
+            (1 / eip[i, iy, ix] / (1 / eip[i, iy, ix] + mu[i, iy, ix])) * (phi[i]^2 *
+            mm[i, ix, iy] / r)
+        }
+      }
+    }
+  }
+
+  env_index <- sqrt(env_index[1, , ] + env_index[2, , ])
+  
+  dim(env_index) <- dim(temp)
+  return(env_index)
+}
+
+
+.liuhelmerssohn <- function(temp) {
+  if (is.na(temp) == TRUE) {
+    env_index <- NA 
+  } else {        
+    if (temp < 12.4 || temp > 32) {
+      a <- 0
+    } else {
+      a <- 0.0043 * temp + 0.0943
+    }
+    if (temp < 12.4 || temp > 32.5) {
+      b <- 0
+    } else if (temp >= 12.4 && temp <= 26.1) {
+      b <- 0.0729 * temp - 0.9037
+    } else {
+      b <- 1
+    }
+    if (temp < 12.3 || temp > 32.5) {
+      c <- 0
+    } else {
+      c <- 0.001044 * temp * (temp - 12.286) * (32.461 - temp)^(1 / 2)
+      if (is.nan(c)) {
+        c <- 0
+      }
+    }
+    bc <- b * c
+    if (temp < 14 || temp > 32) {
+      mu <- 0.04168548
+    } else {
+      mu <- 0.8692 - 0.1590 * temp + 0.01116 * temp^2 - 3.408e-04 * temp^3 + 3.809e-06 * temp^4
+    }
+    if (temp < 12 || temp > 36) {
+      eip <- 0
+    } else {
+      eip <- 4 + exp(-0.123 * temp + 5.15)
+    }
+    env_index <- (a^2 * bc * exp(-mu * eip) / mu)
+  }
+  dim(env_index) <- dim(temp)
+  return(env_index)
+}
+
+.mordecai <- function(temp) {
+  ec <- 0.00001
+  rr <- 8.3144598 
+
+  if (is.na(temp) == TRUE) {
+    env_index <- NA
+  } else {
+    a <- -5.4 + 1.8 * temp - 0.2124 * temp^2 + 0.01015 * temp^3 - 1.515e-04 * temp^4
+    if (a < 0) {
+        a <- 0
+    }
+    b <- 0.0729 * temp - 0.97
+    if (b < 0) {
+        b <- 0
+    }
+    c <- 1.044e-03 * temp * (temp - 12.286) * (32.461 - temp)^(1 / 2)
+    if (is.nan(c)) {
+        c <- 0
+    }
+    bc <- b * c
+    x1 <- (0.8692 - 0.1599 * temp + 0.01116 * temp^2 - 3.408e-04 * temp^3 + 3.809e-06 * temp^4)
+    if (x1 < 0.01) {
+        x1 <- 0.01
+    }
+    lf <- 1 / x1
+    mu <- 1 / (lf + ec)
+    e2a <- exp(-(2.13 - 0.3787 * temp + 0.02457 * temp^2 - 6.778e-04 * temp^3 + 6.794e-06 * temp^4))
+    mdr <- 0.131 - 0.05723 * temp + 0.01164 * temp^2 - 0.001341 * temp^3 +
+      0.00008723 * temp^4 - 3.017e-06 * temp^5 + 5.153e-08 * temp^6 - 3.42e-10 * temp^7
+    pdr <- (1 + exp((6.203e21) / rr * (1 / (-2.176e30) - 1 / (temp + 273.15)))) /
+      ((3.3589e-03 * (temp + 273.15)) / 298 * exp((1.500 / rr) * (1 / 298 - 1 / (temp + 273.15))))
+    if (temp < 8.02 || temp > 35.65) {
+        efd <- 0
+    } else {
+        efd <- 4.88E-02 * (temp - 8.02) * (35.65 - temp)^0.5
+    }
+    env_index <- sqrt((a^2 * bc * (efd * e2a * mdr / (mu)^2) *
+      exp((-mu / (pdr + ec)))) / (mu)) / 1000
+  }
+  dim(env_index) <- dim(temp)           
+  return(env_index)
+}
+
+.wesolowski <- function(temp) {
+  # Setting parameter values for Wesolowski et al., 2015:
+  ec <- 0.00001
+  rr <- 8.3144598
+
+  if (is.na(temp) == TRUE) {
+    env_index <- NA
+  } else {
+    a <- -5.4 + 1.8 * temp - 0.2124 * temp^2 + 0.01015 * temp^3 - 1.515e-04 * temp^4
+    if (a < 0) {
+      a <- 0
+    }
+    b <- 0.0729 * temp - 0.97
+    if (b < 0) {
+      b <- 0
+    }
+    c <- 1.044e-03 * temp * (temp - 12.286) * (32.461 - temp)^(1 / 2)
+    if (is.nan(c)) {
+      c <- 0
+    }
+    bc <- b * c
+    e2a <- exp(-(2.13 - 0.3787 * temp + 0.02457 * temp^2 - 6.778e-04 * temp^3 + 6.794e-06 * temp^4))
+    x1 <- (0.8692 - 0.1599 * temp + 0.01116 * temp^2 - 3.408e-04 * temp^3 + 3.809e-06 * temp^4)
+    if (x1 < 0.01) {
+      x1 <- 0.01
+    }
+    lf <- 1 / x1
+    mdr <- 0.131 - 0.05723 * temp + 0.01164 * temp^2 - 0.001341 * temp^3 +
+      0.00008723 * temp^4 - 3.017e-06 * temp^5 + 5.153e-08 * temp^6 - 3.42e-10 * temp^7
+    pdr <- (1 + exp((6.203e21) / rr * (1 / (-2.176e30) - 1 / (temp + 273.15)))) / 
+      ((3.3589e-03 * (temp + 273.15)) / 298 * exp((1.500 / rr) * (1 / 298 - 1 / (temp + 273.15))))
+    mu <- 1 / (lf + ec)
+    env_index <- sqrt((a^2 * bc * (e2a * mdr / (mu)^2) * exp((-mu / (pdr + ec)))) / (mu)) / 1000
+  }
+  
+  dim(env_index) <- dim(temp)
+  return(env_index)
+}
+
+.empirical <- function(temp, empirical) {
+  env_index <- ifelse(
+    is.na(temp) == TRUE || temp < min(empirical$temp) || temp > max(empirical$temp),
+    NA,
+    sqrt(empirical$aegy[which.min(abs(empirical$temp - temp)) / 10] + 
+      empirical$albo[which.min(abs(empirical$temp - temp))]
+    )
+  )
+  return(env_index)
+}
\ No newline at end of file
diff --git a/inst/empirical.txt b/inst/empirical.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ad70549d6c5eadb6532b2f3f7634812de56ad1b1
--- /dev/null
+++ b/inst/empirical.txt
@@ -0,0 +1,402 @@
+temp	albo	aegy
+1	5	0	0
+2	5.1	0	0
+3	5.2	0	0
+4	5.3	0	0
+5	5.4	0	0
+6	5.5	0	0
+7	5.6	0	0
+8	5.7	0	0
+9	5.8	0	0
+10	5.9	0	0
+11	6	0	0
+12	6.1	0	0
+13	6.2	0	0
+14	6.3	0	0
+15	6.4	0	0
+16	6.5	0	0
+17	6.6	0	0
+18	6.7	0	0
+19	6.8	0	0
+20	6.9	0	0
+21	7	0	0
+22	7.1	0	0
+23	7.2	0	0
+24	7.3	0	0
+25	7.4	0	0
+26	7.5	0	0
+27	7.6	0	0
+28	7.7	0	0
+29	7.8	0	0
+30	7.9	0	0
+31	8	0	0
+32	8.1	0	0
+33	8.2	0	0
+34	8.3	0	0
+35	8.4	0	0
+36	8.5	0	0
+37	8.6	0	0
+38	8.7	0	0
+39	8.8	0	0
+40	8.9	0	0
+41	9	0	0
+42	9.1	0	0
+43	9.2	0	0
+44	9.3	0	0
+45	9.4	0	0
+46	9.5	0	0
+47	9.6	0	0
+48	9.7	0	0
+49	9.8	0	0
+50	9.9	0	0
+51	10	0	0
+52	10.1	0	7.26701110652832e-48
+53	10.2	0	1.99555189437007e-22
+54	10.3	0	6.49268415912656e-17
+55	10.4	0	1.95844479244999e-14
+56	10.5	0	5.39807607442572e-13
+57	10.6	0	4.90217958462951e-12
+58	10.7	0	2.4221021987658e-11
+59	10.8	0	8.25443535682249e-11
+60	10.9	0	2.20328093785689e-10
+61	11	0	5.04805611569711e-09
+62	11.1	0	1.47189779875918e-08
+63	11.2	0	3.3027529676603e-08
+64	11.3	0	6.46491099336551e-08
+65	11.4	0	1.15805258704875e-07
+66	11.5	0	1.94940525745036e-07
+67	11.6	0	3.14152945195542e-07
+68	11.7	0	4.9262890969407e-07
+69	11.8	0	7.62645770234994e-07
+70	11.9	0	1.16715912590015e-06
+71	12	0	1.75864001085107e-06
+72	12.1	0	2.60649603407274e-06
+73	12.2	0	3.94184080212757e-06
+74	12.3	3.14196404388935e-11	5.91210558178024e-06
+75	12.4	1.74869754004445e-10	8.79845154799995e-06
+76	12.5	2.24631012705507e-09	1.29591630826655e-05
+77	12.6	5.58235500971744e-09	2.08712773504148e-05
+78	12.7	1.10978416373972e-08	3.20826062480237e-05
+79	12.8	1.97331328557588e-08	5.10980761726375e-05
+80	12.9	3.26541362019039e-08	7.88913411592989e-05
+81	13	5.12532355034818e-08	0.000118201038481699
+82	13.1	7.72678881529682e-08	0.000174741278817335
+83	13.2	1.12860894254321e-07	0.000252640595988958
+84	13.3	1.68152404814637e-07	0.000362519239515314
+85	13.4	2.59159535417229e-07	0.00052268159339204
+86	13.5	4.34530126853778e-07	0.000737140258682322
+87	13.6	7.15988308077365e-07	0.00104280674254481
+88	13.7	1.10810098721519e-06	0.00146460594796231
+89	13.8	1.64030150318549e-06	0.00202682273943419
+90	13.9	2.56447091004182e-06	0.00274991529573702
+91	14	4.08669454479855e-06	0.00372238707325898
+92	14.1	6.42056166112727e-06	0.00498640748720935
+93	14.2	9.8970442092177e-06	0.00660177833279288
+94	14.3	1.44139370380086e-05	0.00874706177757016
+95	14.4	2.04868382403983e-05	0.0114536517778624
+96	14.5	2.87863306639688e-05	0.0148613798900532
+97	14.6	4.13493369784319e-05	0.0191832994538838
+98	14.7	5.80910531620301e-05	0.0246408789386216
+99	14.8	7.93800574474559e-05	0.0313570909381335
+100	14.9	0.000107828679535155	0.0397301246553432
+101	15	0.000143847478611779	0.0498428444244661
+102	15.1	0.000191117553287409	0.0621061468436647
+103	15.2	0.000254380558226787	0.0770132747433837
+104	15.3	0.00034202322685854	0.0952253602687086
+105	15.4	0.000451635632245419	0.116874203279141
+106	15.5	0.000592076173536784	0.142671816526039
+107	15.6	0.000769509605760271	0.173417177106153
+108	15.7	0.000984393737706117	0.209227433161616
+109	15.8	0.0012544157102959	0.251561865825694
+110	15.9	0.00159051318850117	0.301168051653919
+111	16	0.00201429434708981	0.359663734290125
+112	16.1	0.00252686222939845	0.427043260287913
+113	16.2	0.00317095547702531	0.505404019648344
+114	16.3	0.00396215825150203	0.596231144374168
+115	16.4	0.0048936906552879	0.701217652217518
+116	16.5	0.00604198293128486	0.820742880431283
+117	16.6	0.00739391904370352	0.956997709422477
+118	16.7	0.00897739530805604	1.11330467190701
+119	16.8	0.0109065790526767	1.28957967110784
+120	16.9	0.0131861392416037	1.48702508369891
+121	17	0.0158316862248661	1.71142201857145
+122	17.1	0.0189476905111509	1.9623756647912
+123	17.2	0.022639307696192	2.24050675673134
+124	17.3	0.0268943166114979	2.54998888644124
+125	17.4	0.0318364999193829	2.89213869180169
+126	17.5	0.0375214720828055	3.26748703411117
+127	17.6	0.0440315392492445	3.68414641881988
+128	17.7	0.051555332918417	4.14378322016399
+129	17.8	0.0602107062395169	4.64550462458758
+130	17.9	0.0699135999208867	5.18825849367906
+131	18	0.0810609234541401	5.77924231434693
+132	18.1	0.0936677276338084	6.41548493466512
+133	18.2	0.107838363512702	7.1063503376978
+134	18.3	0.123715421678103	7.8485876264103
+135	18.4	0.141356135955452	8.64426113944596
+136	18.5	0.161244645754316	9.4939047128634
+137	18.6	0.183309753716304	10.4040956510328
+138	18.7	0.207752947471398	11.3807669456181
+139	18.8	0.234670715499831	12.4207563412324
+140	18.9	0.26419692592572	13.5223346061691
+141	19	0.29699990724692	14.6891613895531
+142	19.1	0.332631451991198	15.923770453599
+143	19.2	0.37130383933373	17.2257024329454
+144	19.3	0.413766492094605	18.6000188958484
+145	19.4	0.459718181286087	20.04473975
+146	19.5	0.509203366743798	21.5661191321618
+147	19.6	0.562754149511879	23.1641220305504
+148	19.7	0.620492095414727	24.8377622195195
+149	19.8	0.682424407259472	26.5910937379638
+150	19.9	0.749150939428306	28.4288385027613
+151	20	0.82063142310609	30.3427447273267
+152	20.1	0.896975013345537	32.3340092913406
+153	20.2	0.978051229537782	34.4060799165644
+154	20.3	1.06374957602415	36.5588895302359
+155	20.4	1.15490698795063	38.7899217711649
+156	20.5	1.25111960115299	41.1046694926988
+157	20.6	1.35316926880471	43.4973257200633
+158	20.7	1.46098734000124	45.9696199506604
+159	20.8	1.57437924118154	48.5246117952479
+160	20.9	1.6935551096107	51.1582303539374
+161	21	1.81893495933483	53.8654831044457
+162	21.1	1.94990059351414	56.6532566867408
+163	21.2	2.0870333512605	59.5198331799156
+164	21.3	2.2299952500153	62.4602212465642
+165	21.4	2.37864753830364	65.4747158378878
+166	21.5	2.53335326546983	68.5638350079553
+167	21.6	2.69455009100578	71.7222549514207
+168	21.7	2.86179245035676	74.9546454224119
+169	21.8	3.03550730302727	78.2548143857356
+170	21.9	3.2151993291081	81.6193175802761
+171	22	3.4010035274587	85.040854275288
+172	22.1	3.59293553429209	88.5176675887392
+173	22.2	3.79126782028487	92.0552318953705
+174	22.3	3.99615926991393	95.6413974847115
+175	22.4	4.20765886298318	99.2699042340408
+176	22.5	4.42561189737981	102.941931576309
+177	22.6	4.64958321865548	106.65158864424
+178	22.7	4.87928762475173	110.397704591579
+179	22.8	5.11491798479008	114.170952399813
+180	22.9	5.35656693909105	117.967125697487
+181	23	5.60424233984691	121.783214145687
+182	23.1	5.85749394466384	125.608739367861
+183	23.2	6.11633496681088	129.437577367979
+184	23.3	6.38127969380391	133.264653172263
+185	23.4	6.65133086468265	137.084137733188
+186	23.5	6.92649769944134	140.889850255725
+187	23.6	7.20715429194241	144.674826310175
+188	23.7	7.49276142498194	148.43228330375
+189	23.8	7.78295498628382	152.155290146054
+190	23.9	8.07756157517245	155.836547381111
+191	24	8.37642121921268	159.471097491766
+192	24.1	8.67932310888436	163.048529054954
+193	24.2	8.98603357181722	166.563234144951
+194	24.3	9.29630131312483	170.006420206781
+195	24.4	9.60985934098338	173.368361949881
+196	24.5	9.92642591868661	176.640875935887
+197	24.6	10.2459259556872	179.816365690908
+198	24.7	10.5677723280057	182.886565391347
+199	24.8	10.8916563978521	185.843141407777
+200	24.9	11.2172677351661	188.677691082182
+201	25	11.5442634945313	191.382106168346
+202	25.1	11.8722880706697	193.948387755488
+203	25.2	12.2009702833037	196.36868193318
+204	25.3	12.5299266666884	198.634584993542
+205	25.4	12.8587622229721	200.738835006157
+206	25.5	13.1870710291405	202.673765392136
+207	25.6	13.514436761053	204.431072523135
+208	25.7	13.8404328706676	206.003416663432
+209	25.8	14.1646246480312	207.383063555628
+210	25.9	14.486568138682	208.563034819444
+211	26	14.8058112096168	209.536812138252
+212	26.1	15.1218947655244	210.298229912246
+213	26.2	15.4343527894637	210.841148332732
+214	26.3	15.7427130041459	211.15917058413
+215	26.4	16.0464920733374	211.247499381759
+216	26.5	16.3452101195749	211.101210401603
+217	26.6	16.6383854826497	210.715968526715
+218	26.7	16.925518528805	210.087730996625
+219	26.8	17.2061235687347	209.212766619546
+220	26.9	17.4797155533508	208.088131823885
+221	27	17.745802801227	206.711500672635
+222	27.1	18.0038910102781	205.081224277892
+223	27.2	18.2534850859162	203.196691590182
+224	27.3	18.4941076653118	201.057901900131
+225	27.4	18.725273273391	198.665035623928
+226	27.5	18.9465004904626	196.019405957757
+227	27.6	19.1573231937636	193.122935337843
+228	27.7	19.357276024177	189.979066062919
+229	27.8	19.5458952949288	186.590943550343
+230	27.9	19.7227284847839	182.962955087303
+231	28	19.8873481259785	179.100494482392
+232	28.1	20.0393225357807	175.011023865686
+233	28.2	20.178237206428	170.702404015189
+234	28.3	20.3036880173062	166.183121576334
+235	28.4	20.4152984184343	161.463520022466
+236	28.5	20.5126987352094	156.554353250729
+237	28.6	20.5955091781159	151.467765490333
+238	28.7	20.6634112982833	146.217030869804
+239	28.8	20.7160939928104	140.816223457413
+240	28.9	20.7532646585221	135.278428378365
+241	29	20.774641880637	129.619957604372
+242	29.1	20.7799638146496	123.856016776917
+243	29.2	20.7690295697177	118.003193567047
+244	29.3	20.7416475643345	112.077293738837
+245	29.4	20.6976366852697	106.096583090103
+246	29.5	20.6368446229951	100.079978575113
+247	29.6	20.5591501809307	94.0453698383962
+248	29.7	20.4644779699975	88.0155408797388
+249	29.8	20.3527679902312	82.0135143455154
+250	29.9	20.2239788853992	76.0662415191356
+251	30	20.0781063994079	70.1999178737101
+252	30.1	19.9151925801923	64.4433834812501
+253	30.2	19.7353198356195	58.8309232642251
+254	30.3	19.5385964916249	53.3983959309562
+255	30.4	19.3251625554253	48.1796194737268
+256	30.5	19.0951915630466	43.209247632402
+257	30.6	18.8488915647271	38.517378654949
+258	30.7	18.586504796528	34.1377388746358
+259	30.8	18.3083213865769	30.092813027491
+260	30.9	18.0146571226604	26.3975462375243
+261	31	17.7058557875912	23.0524962479585
+262	31.1	17.3823110989741	20.0560568082843
+263	31.2	17.0444441735365	17.3982637698392
+264	31.3	16.692717575385	15.053439626178
+265	31.4	16.3276230882856	12.993654190201
+266	31.5	15.9496850151796	11.1933264577378
+267	31.6	15.5594694425189	9.63258167302814
+268	31.7	15.1575751559031	8.29088515658698
+269	31.8	14.7446193527128	7.14133562536133
+270	31.9	14.3212714662343	6.15800896561084
+271	32	13.8881980732254	5.31939562492462
+272	32.1	13.4461166367272	4.60395079584917
+273	32.2	12.9957307751225	3.9922549126184
+274	32.3	12.5378165657429	3.46954928294953
+275	32.4	12.0731334650291	3.02233191588653
+276	32.5	11.6024831371139	2.63799900315302
+277	32.6	11.1266567352906	2.30747189250574
+278	32.7	10.6464346629302	2.02224150552596
+279	32.8	10.1626096844356	1.77374773541624
+280	32.9	9.67596223870091	1.55668473935375
+281	33	9.18721442726549	1.36626597490556
+282	33.1	8.69704258421091	1.19853358870185
+283	33.2	8.20603325615985	1.04925005682019
+284	33.3	7.71462359591067	0.91616624460685
+285	33.4	7.22306632238162	0.797486044480323
+286	33.5	6.73126370061339	0.691663754527126
+287	33.6	6.23856150854693	0.597444874799824
+288	33.7	5.74327051132003	0.513909887386837
+289	33.8	5.24150056025327	0.440022554554648
+290	33.9	4.7234608219282	0.37486440222695
+291	34	4.14441956459301	0.317555422973365
+292	34.1	3.3457934190837	0.267398700795832
+293	34.2	2.64029689947616	0.223895735004689
+294	34.3	2.05500762066432	0.186240986366627
+295	34.4	1.59549189487488	0.153691153583648
+296	34.5	1.21811975591742	0.125777288386402
+297	34.6	0.926211388416226	0.102096327349633
+298	34.7	0.703560252017681	0.0821732998307432
+299	34.8	0.525227625905333	0.0653912657250115
+300	34.9	0.398087337792727	0.0514715626511726
+301	35	0.301041486351743	0.0399938358044144
+302	35.1	0.221920551620596	0.0303490829527474
+303	35.2	0.158399308161674	0.0204050454865184
+304	35.3	0.11137240595035	0.0121926575691942
+305	35.4	0.0750022178641641	0.00871791527363357
+306	35.5	0.0497161680440828	0.00591367289041522
+307	35.6	0.0323922259513553	0.00393430854732474
+308	35.7	0.0196561006423747	0.00228726097768791
+309	35.8	0.0122379911295625	0.000521192531905449
+310	35.9	0.00753197877487249	0.000299158578473394
+311	36	0.00390458337631893	0.000138063326459654
+312	36.1	0.0023713454259117	2.11728627290423e-05
+313	36.2	0.00117379654065915	4.00672120858184e-06
+314	36.3	0.000512226871542284	1.06188790609357e-07
+315	36.4	0.000144970836057746	0
+316	36.5	5.75924730218119e-05	0
+317	36.6	2.37774473058921e-05	0
+318	36.7	1.78514268725238e-06	0
+319	36.8	0	0
+320	36.9	0	0
+321	37	0	0
+322	37.1	0	0
+323	37.2	0	0
+324	37.3	0	0
+325	37.4	0	0
+326	37.5	0	0
+327	37.6	0	0
+328	37.7	0	0
+329	37.8	0	0
+330	37.9	0	0
+331	38	0	0
+332	38.1	0	0
+333	38.2	0	0
+334	38.3	0	0
+335	38.4	0	0
+336	38.5	0	0
+337	38.6	0	0
+338	38.7	0	0
+339	38.8	0	0
+340	38.9	0	0
+341	39	0	0
+342	39.1	0	0
+343	39.2	0	0
+344	39.3	0	0
+345	39.4	0	0
+346	39.5	0	0
+347	39.6	0	0
+348	39.7	0	0
+349	39.8	0	0
+350	39.9	0	0
+351	40	0	0
+352	40.1	0	0
+353	40.2	0	0
+354	40.3	0	0
+355	40.4	0	0
+356	40.5	0	0
+357	40.6	0	0
+358	40.7	0	0
+359	40.8	0	0
+360	40.9	0	0
+361	41	0	0
+362	41.1	0	0
+363	41.2	0	0
+364	41.3	0	0
+365	41.4	0	0
+366	41.5	0	0
+367	41.6	0	0
+368	41.7	0	0
+369	41.8	0	0
+370	41.9	0	0
+371	42	0	0
+372	42.1	0	0
+373	42.2	0	0
+374	42.3	0	0
+375	42.4	0	0
+376	42.5	0	0
+377	42.6	0	0
+378	42.7	0	0
+379	42.8	0	0
+380	42.9	0	0
+381	43	0	0
+382	43.1	0	0
+383	43.2	0	0
+384	43.3	0	0
+385	43.4	0	0
+386	43.5	0	0
+387	43.6	0	0
+388	43.7	0	0
+389	43.8	0	0
+390	43.9	0	0
+391	44	0	0
+392	44.1	0	0
+393	44.2	0	0
+394	44.3	0	0
+395	44.4	0	0
+396	44.5	0	0
+397	44.6	0	0
+398	44.7	0	0
+399	44.8	0	0
+400	44.9	0	0
+401	45	0	0
diff --git a/tests/testthat/test-RNoughtIndices.R b/tests/testthat/test-RNoughtIndices.R
new file mode 100644
index 0000000000000000000000000000000000000000..5989e4d39d1a243db0b0b36650e71ffb2ed17268
--- /dev/null
+++ b/tests/testthat/test-RNoughtIndices.R
@@ -0,0 +1,486 @@
+# ============================================================================
+# COMPREHENSIVE UNIT TESTS FOR RNoughtIndices
+# ============================================================================
+# 
+# This test suite validates the integrity and functionality of:
+# - CST_RNoughtIndices (s2dv_cube wrapper)
+# - RNoughtIndices (main computation function) 
+#
+# Author: Javier Corvillo, \email{javier.corvillo@bsc.es}
+# ============================================================================
+
+# Load required libraries
+library(testthat)
+library(multiApply)
+library(CSIndicators)
+source("/esarchive/scratch/jcorvill/gitlab/bsc_functions/RNoughtIndices.R")
+
+# ============================================================================
+# INTEGRITY ANALYSIS FINDINGS
+# ============================================================================
+
+cat("=== FUNCTION INTEGRITY ANALYSIS ===\n")
+
+# ============================================================================
+# UTILITY FUNCTIONS FOR TESTING
+# ============================================================================
+
+# Create synthetic temperature data for testing
+create_test_temp_data <- function(dims = c(time = 100, lat = 5, lon = 4), 
+                                  temp_range = c(15, 35), 
+                                  add_nas = FALSE) {
+  n_total <- prod(dims)
+  temp_data <- array(
+    runif(n_total, min = temp_range[1], max = temp_range[2]),
+    dim = dims
+  )
+  
+  if (add_nas) {
+    # Add some NA values (5% of data)
+    na_indices <- sample(1:n_total, size = floor(0.05 * n_total))
+    temp_data[na_indices] <- NA
+  }
+  
+  return(temp_data)
+}
+
+# Create synthetic s2dv_cube object
+create_test_s2dv_cube <- function(dims = c(member = 2, sdate = 2, time = 100, lat = 5, lon = 4)) {
+  temp_cube <- list()
+  temp_cube$data <- create_test_temp_data(dims)
+  
+  # Create coordinates
+  temp_cube$coords <- list(
+    lat = seq(40, 44, length.out = dims["lat"]),
+    lon = seq(-5, -1, length.out = dims["lon"])
+  )
+  
+  # Create metadata
+  variable <- list(varName = 'temp', 
+                   metadata = list(temp = list(level = 'Surface')))
+  temp_cube$attrs <- list(
+    Variable = variable, 
+    Datasets = 'synthetic_test', 
+    when = Sys.time()
+  )
+  
+  # Create dates
+  if ("time" %in% names(dims)) {
+    n_time <- dims["time"]
+    dates <- seq(as.Date("2020-01-01"), by = "day", length.out = n_time)
+    if ("sdate" %in% names(dims)) {
+      dates <- array(rep(dates, dims["sdate"]), 
+                     dim = c(time = n_time, sdate = dims["sdate"]))
+    }
+    temp_cube$attrs$Dates <- dates
+  }
+  
+  class(temp_cube) <- 's2dv_cube'
+  return(temp_cube)
+}
+
+# Create test Kramer probability matrix for Caminade method
+create_test_kramer_matrix <- function(lat_dim = 5, lon_dim = 4) {
+  mm <- array(runif(2 * lat_dim * lon_dim, 0.1, 1.0), 
+              dim = c(2, lon_dim, lat_dim))
+  dimnames(mm) <- list(NULL, 
+                       paste0("lon_", 1:lon_dim), 
+                       paste0("lat_", 1:lat_dim))
+  names(dim(mm)) <- c("", "lon", "lat")
+  return(mm)
+}
+
+# ============================================================================
+# PARAMETER VALIDATION
+# ============================================================================
+
+test_parameter_validation <- function() {
+  cat("\n🧪 TESTING PARAMETER VALIDATION\n")
+  
+  # Test 1: Invalid method parameter
+  test_that("Invalid method throws error", {
+    temp_data <- create_test_temp_data()
+    expect_error(
+      RNoughtIndices(temp_data, method = "invalid_method", mm = NULL, 
+                     lon_dim = NULL, lat_dim = NULL),
+      "Invalid value for method"
+    )
+  })
+  
+  # Test 2: Caminade method without Kramer matrix
+  test_that("Caminade method requires Kramer matrix", {
+    temp_data <- create_test_temp_data()
+    expect_error(
+      RNoughtIndices(temp_data, method = "caminade", mm = NULL, 
+                     lon_dim = "lon", lat_dim = "lat"),
+      "no Kramer probability matrix"
+    )
+  })
+  
+  # Test 3: CST wrapper with non-s2dv_cube
+  test_that("CST wrapper requires s2dv_cube", {
+    temp_data <- create_test_temp_data()
+    expect_error(
+      CST_RNoughtIndices(temp_data, method = "mordecai", mm = NULL, 
+                         lon_dim = NULL, lat_dim = NULL),
+      "must be of the class 's2dv_cube'"
+    )
+  })
+  
+  cat("✅ Parameter validation tests completed\n")
+}
+
+# ============================================================================
+# MATHEMATICAL MODELS
+# ============================================================================
+
+test_mathematical_models <- function() {
+  cat("\n🧮 TESTING MATHEMATICAL MODELS\n")
+  
+  # Test temperatures for validation
+  test_temps <- c(10, 15, 20, 25, 30, 35, 40)
+  
+  test_that("Liu-Helmersson model produces reasonable outputs", {
+    for (temp in test_temps) {
+      result <- .liuhelmerssohn(temp)
+      
+      # Should return non-negative values for reasonable temperatures
+      if (temp >= 12.4 && temp <= 32.5) {
+        expect_true(result >= 0, 
+                    info = paste("Negative result for temp", temp))
+      }
+      
+      # Should return 0 for extreme temperatures
+      if (temp < 12.4 || temp > 32.5) {
+        expect_equal(result, 0, 
+                     info = paste("Non-zero result for extreme temp", temp))
+      }
+    }
+  })
+  
+  test_that("Mordecai model produces reasonable outputs", {
+    for (temp in test_temps) {
+      result <- .mordecai(temp)
+      
+      # Should be non-negative
+      expect_true(result >= 0, 
+                  info = paste("Negative result for temp", temp))
+      
+      # Should be finite for reasonable temperatures
+      if (temp >= 5 && temp <= 40) {
+        expect_true(is.finite(result), 
+                    info = paste("Non-finite result for temp", temp))
+      }
+    }
+  })
+  
+  test_that("Wesolowski model produces reasonable outputs", {
+    for (temp in test_temps) {
+      result <- .wesolowski(temp)
+      
+      # Should be non-negative
+      expect_true(result >= 0, 
+                  info = paste("Negative result for temp", temp))
+      
+      # Should be finite for reasonable temperatures
+      if (temp >= 5 && temp <= 40) {
+        expect_true(is.finite(result), 
+                    info = paste("Non-finite result for temp", temp))
+      }
+    }
+  })
+  
+  test_that("Caminade model produces reasonable outputs", {
+    # Create test temperature matrix
+    temp_matrix <- matrix(test_temps, nrow = length(test_temps), ncol = 1)
+    mm <- array(0.5, dim = c(2, ncol(temp_matrix), nrow(temp_matrix)))
+    
+    result <- .caminade(temp_matrix, mm, 
+                        nlon = ncol(temp_matrix),    
+                        nlat = nrow(temp_matrix)
+    )
+    # Should be non-negative
+    expect_true(all(result >= 0, na.rm = TRUE))
+    
+    # Should have same dimensions as input
+    expect_equal(dim(result), dim(temp_matrix))
+  })
+  
+  cat("✅ Mathematical model tests completed\n")
+}
+
+# ============================================================================
+# EDGE CASES
+# ============================================================================
+
+test_edge_cases <- function() {
+  cat("\n🔬 TESTING EDGE CASES\n")
+  
+  test_that("NA temperature handling", {
+    # Test each model with NA input
+    temp_na <- NA_real_
+    
+    expect_true(is.na(.liuhelmerssohn(temp_na)))
+    expect_true(is.na(.mordecai(temp_na)))
+    expect_true(is.na(.wesolowski(temp_na)))
+    
+    # Test Caminade with NA
+    temp_matrix <- matrix(c(25, NA, 30), nrow = 3, ncol = 1)
+    mm <- array(0.5, dim = c(2, 1, 3))
+    result <- .caminade(temp_matrix, mm, nlat = 3, nlon = 1)
+    expect_true(is.na(result[2, 1]))
+    expect_true(!is.na(result[1, 1]))
+    expect_true(!is.na(result[3, 1]))
+  })
+  
+  test_that("Extreme temperature handling", {
+    extreme_temps <- c(-50, -10, 0, 50, 60, 100)
+    
+    for (temp in extreme_temps) {
+      # Functions should not crash and should return finite or zero values
+      result_liu <- .liuhelmerssohn(temp)
+      result_mordecai <- .mordecai(temp)
+      result_wesolowski <- .wesolowski(temp)
+      
+      expect_true(is.finite(result_liu) || result_liu == 0)
+      expect_true(is.finite(result_mordecai) || result_mordecai == 0)  
+      expect_true(is.finite(result_wesolowski) || result_wesolowski == 0)
+    }
+  })
+  
+  test_that("Array dimension preservation", {
+    # Test that output dimensions match input dimensions
+    temp_array <- create_test_temp_data(c(time = 50, lat = 3, lon = 2))
+    
+    result <- RNoughtIndices(temp_array, method = "mordecai", mm = NULL,
+                             lon_dim = NULL, lat_dim = NULL)
+    
+    expect_equal(dim(result), dim(temp_array))
+  })
+  
+  test_that("Empirical method with edge values", {
+    # Test empirical method with temperatures at the edge of data range
+    # First, let's read the empirical data to understand its range
+    empirical_file <- system.file("index_empirical", "empirical.txt", 
+                                  package = "CSIndicators", mustWork = FALSE)
+    
+    if (file.exists(empirical_file)) {
+      empirical <- read.table(empirical_file)
+      
+      # Test with temperatures at the boundary
+      min_temp <- min(empirical$temp)
+      max_temp <- max(empirical$temp)
+      
+      result_min <- .empirical(min_temp, empirical)
+      result_max <- .empirical(max_temp, empirical)
+      
+      expect_true(is.finite(result_min) && result_min >= 0)
+      expect_true(is.finite(result_max) && result_max >= 0)
+      
+      # Test with temperatures outside range
+      result_below <- .empirical(min_temp - 1, empirical)
+      result_above <- .empirical(max_temp + 1, empirical)
+      
+      expect_true(is.na(result_below))
+      expect_true(is.na(result_above))
+    } else {
+      cat("⚠️  Empirical data file not found, skipping empirical edge case tests\n")
+    }
+  })
+  
+  cat("✅ Edge case tests completed\n")
+}
+
+# ============================================================================
+# INTEGRATION TEST
+# ============================================================================
+
+test_integration <- function() {
+  cat("\n🔧 TESTING INTEGRATION\n")
+  
+  test_that("Full CST_RNoughtIndices workflow", {
+    # Test complete workflow with s2dv_cube
+    temp_cube <- create_test_s2dv_cube()
+    
+    # Test different methods
+    methods <- c("mordecai", "liuhelmerssohn", "wesolowski")
+    
+    for (method in methods) {
+      result <- CST_RNoughtIndices(temp_cube, method = method, mm = NULL,
+                                   lon_dim = NULL, lat_dim = NULL)
+      
+      # Should return s2dv_cube
+      expect_s3_class(result, "s2dv_cube")
+      
+      # Should have same structure as input
+      expect_equal(names(result), names(temp_cube))
+      
+      # Variable name should be updated
+      expect_equal(result$attrs$Variable$varName, "RNoughtIndex")
+      
+      # Data should be numeric and finite (where not NA)
+      finite_data <- result$data[is.finite(result$data)]
+      expect_true(all(finite_data >= 0))
+    }
+  })
+  
+  test_that("Caminade method integration", {
+    # Create compatible data for Caminade method
+    temp_dims <- c(time = 50, lat = 5, lon = 4)
+    temp_data <- create_test_temp_data(temp_dims)
+    mm <- create_test_kramer_matrix(lat_dim = 5, lon_dim = 4)
+    
+    result <- RNoughtIndices(temp_data, method = "caminade", mm = mm,
+                             lon_dim = "lon", lat_dim = "lat")
+    
+    # Should preserve input dimensions
+    expect_equal(dim(result), dim(temp_data))
+    
+    # Should produce non-negative finite values
+    finite_values <- result[is.finite(result)]
+    expect_true(all(finite_values >= 0))
+  })
+  
+  test_that("Temporal subsetting functionality", {
+    temp_cube <- create_test_s2dv_cube()
+    
+    # Test with start/end parameters
+    result <- CST_RNoughtIndices(temp_cube, method = "mordecai", mm = NULL,
+                                 lon_dim = NULL, lat_dim = NULL)
+    
+    # Should return s2dv_cube with potentially reduced time dimension
+    expect_s3_class(result, "s2dv_cube")
+    expect_true(is.array(result$data))
+  })
+  
+  cat("✅ Integration tests completed\n")
+}
+
+# ============================================================================
+# PERFORMANCE TEST
+# ============================================================================
+
+test_performance <- function() {
+  cat("\n⚡ TESTING PERFORMANCE\n")
+  
+  test_that("Performance with large datasets", {
+    # Test with larger arrays to ensure reasonable performance
+    large_dims <- c(time = 365, lat = 20, lon = 30)  # Realistic climate data size
+    
+    cat("   Testing with", prod(large_dims), "data points...\n")
+    
+    start_time <- Sys.time()
+    temp_data <- create_test_temp_data(large_dims)
+    
+    result <- RNoughtIndices(temp_data, method = "mordecai", mm = NULL,
+                             lon_dim = NULL, lat_dim = NULL)
+    
+    end_time <- Sys.time()
+    computation_time <- as.numeric(end_time - start_time, units = "secs")
+    
+    cat("   Computation time:", round(computation_time, 3), "seconds\n")
+    
+    # Should complete in reasonable time (less than 30 seconds for this size)
+    expect_true(computation_time < 30)
+    
+    # Should produce correct output dimensions
+    expect_equal(dim(result), dim(temp_data))
+  })
+  
+  test_that("Memory efficiency", {
+    # Test that the function doesn't create excessive intermediate objects
+    temp_data <- create_test_temp_data(c(time = 100, lat = 10, lon = 10))
+    
+    # Get initial memory usage
+    initial_memory <- gc()
+    
+    result <- RNoughtIndices(temp_data, method = "mordecai", mm = NULL,
+                             lon_dim = NULL, lat_dim = NULL)
+    
+    # Force garbage collection
+    final_memory <- gc()
+    
+    # Function should not dramatically increase memory usage
+    # (allowing for some overhead)
+    memory_increase <- final_memory[2, 2] - initial_memory[2, 2]
+    expect_true(memory_increase < 100)  # Less than 100MB increase
+  })
+  
+  cat("✅ Performance tests completed\n")
+}
+
+# ============================================================================
+# VALIDATION TEST
+# ============================================================================
+
+test_output_validation <- function() {
+  cat("\n✅ TESTING OUTPUT VALIDATION\n")
+  
+  test_that("R0 values are in reasonable range", {
+    # Test with realistic temperature data (20-30°C)
+    realistic_temps <- seq(20, 30, by = 0.5)
+    
+    methods <- c("mordecai", "liuhelmerssohn", "wesolowski")
+    
+    for (method in methods) {
+      cat("   Testing", method, "method...\n")
+      
+      for (temp in realistic_temps) {
+        temp_array <- array(temp, dim = c(time = 1))
+        
+        result <- RNoughtIndices(temp_array, method = method, mm = NULL,
+                                 lon_dim = NULL, lat_dim = NULL)
+        
+        # R0 values should typically be between 0 and 10 for realistic temperatures
+        expect_true(result >= 0 && result <= 10,
+                    info = paste("Unrealistic R0 value", result, "for temp", temp, "method", method))
+      }
+    }
+  })
+  
+  test_that("Temperature response makes biological sense", {
+    # R0 should generally increase with temperature up to an optimum, then decrease
+    temps <- c(15, 20, 25, 30, 35)
+    temp_array <- array(temps, dim = c(time = length(temps)))
+    
+    result_mordecai <- RNoughtIndices(temp_array, method = "mordecai", mm = NULL,
+                                      lon_dim = NULL, lat_dim = NULL)
+    
+    # At very low and high temperatures, R0 should be low
+    expect_true(result_mordecai[1] < max(result_mordecai))  # 15°C < peak
+    expect_true(result_mordecai[5] < max(result_mordecai))  # 35°C < peak
+    
+    # Peak should be somewhere in the middle range
+    peak_index <- which.max(result_mordecai)
+    expect_true(peak_index >= 2 && peak_index <= 4)  # Peak between 20-30°C
+  })
+  
+  cat("✅ Output validation tests completed\n")
+}
+
+# ============================================================================
+# MAIN TEST RUNNER
+# ============================================================================
+
+run_all_tests <- function() {
+  
+  # Run all test suites
+  tryCatch({
+    test_parameter_validation()
+    test_mathematical_models() 
+    test_edge_cases()
+    test_integration()
+    test_performance()
+    test_output_validation()
+    
+    cat("\n🎉 ALL TESTS COMPLETED SUCCESSFULLY!\n")
+    
+    return(TRUE)
+    
+  }, error = function(e) {
+    cat("\n❌ TEST FAILED:\n")
+    cat("Error:", e$message, "\n")
+    
+    return(FALSE)
+  })
+}
\ No newline at end of file