diff --git a/R/get_pop.R b/R/get_pop.R
index e8c4e5c..ae17a4b 100644
--- a/R/get_pop.R
+++ b/R/get_pop.R
@@ -1,4 +1,31 @@
 # Helper function
+#' Download country-specific constrained world pop 100*100m raster (2015-2030) from the worldpop data repo: 
+#' @param country country iso code
+#' @param year the year of the population raster
+#' @param dest_dir folder name where the downloaded population raster will be saved
+download_worldpop_constrained <- function(country, year, dest_dir = "worldpop") {
+  if (!dir.exists(dest_dir)) dir.create(dest_dir)
+  
+  out_file <- file.path(
+    dest_dir,
+    paste0(tolower(country), "_pop_", year, "_CN_100m_R2025A_v1.tif")
+  )
+  
+  if (file.exists(out_file)) return(out_file)
+  
+  url <- paste0(
+    "https://data.worldpop.org/GIS/Population/Global_2015_2030/",
+    "R2025A/", year, "/", country,
+    "/v1/100m/constrained/",
+    tolower(country), "_pop_", year, "_CN_100m_R2025A_v1.tif"
+  )
+  
+  message("Downloading WorldPop constrained raster: ", url)
+  curl::curl_download(url, out_file)
+  
+  return(out_file)
+}
+
 #' Estimate Adjustment Factors for Population Data (for years >=2021, use the adjustment factors at 2020 instead)
 #'
 #' @param year numeric vector: Years for which population needs to be estimated.
@@ -24,7 +51,7 @@ estimate_adj_factors <- function(
   }
   
   # Adjust year if it exceeds the range
-  year <- if (year > 2020) 2020 else year
+  year <- if (year < 2015) 2015 else year
   
   # Validate year
   if (any(!(year %in% WPP2024$Time))) {
@@ -32,21 +59,12 @@ estimate_adj_factors <- function(
          paste0(range(WPP2024$Time), collapse = "-"))
   }
   
-  # Construct file path for the raster
-  pop_file_path <- file.path(raster_dir, paste0(tolower(country), "_ppp_", year, ".tif"))
-  
-  # Download raster if missing
-  if (!file.exists(pop_file_path)) {
-    message("Downloading raster for ", country, " (", year, ")...")
-    pop_file_path <- wpgpDownloadR::wpgpGetCountryDataset(
-      ISO3 = country,
-      covariate = paste0("ppp_", year),
-      destDir = dest_dir,
-      method = "curl"
-    )
-  }
-  
-  # Load raster and calculate population
+  # download and load population raster
+  pop_file_path <- download_worldpop_constrained(
+    country = country,
+    year = year,
+    dest_dir = raster_dir
+  )
   pop_raster <- raster::raster(pop_file_path)
   
   # Get the country-level shapefile
@@ -73,6 +91,7 @@ estimate_adj_factors <- function(
   return(adj_factors)
 }
 
+#' Estimate population function - updated in May 2026: use the updated population raster between 2015-2030
 #' @export
 #' @title get_pop
 #' @name get_pop
@@ -97,58 +116,37 @@ get_pop <- function(
     dir.create(pop_raster_path)
   }
   
-  pop_year <- if (year > 2020) 2020 else year
+  data("WPP2024", package = "OutbreakExtractR")
   
-  raster_file <- paste0(pop_raster_path,"/",country,"_ppp_",pop_year,'.tif')
-  if(!file.exists(raster_file)){
-    raster_file <- wpgpDownloadR::wpgpGetCountryDataset(ISO3 = country,
-                                                        covariate = paste0("ppp_",pop_year),
-                                                        destDir = pop_raster_path,
-                                                        method = "curl")
-  }
+  country <- toupper(country)
+  year <- as.integer(year)
+  
+  # Get the population raster
+  raster_file <- download_worldpop_constrained(
+    country = country,
+    year = year,
+    dest_dir = pop_raster_path
+  )
+  
+  pop_raster <- raster::raster(raster_file)
+  
+  # Calculate population data
+  pop <- exactextractr::exact_extract(
+    pop_raster,
+    shp$geometry,
+    "sum"
+  )
+  
+  # Extract country-year specific adjusting factors to align with UN pop 2024 version
+  adj_factors <- estimate_adj_factors(
+    country = country,
+    year = year,
+    raster_dir = pop_raster_path,
+    dest_dir = pop_raster_path,
+    country_shp = country_shp
+  )
   
-  if(year > 2020){
-    pop_raster_2020 <- raster::raster(raster_file)
-    
-    # Aligning worldpop estimates to the UN population estiamtes at the country level
-    adj_factors <- estimate_adj_factors(
-      country = country, 
-      year = 2020, 
-      raster_dir = pop_raster_path, 
-      dest_dir = pop_raster_path,
-      country_shp = country_shp
-    )
-    
-    pop_2020 <- exactextractr::exact_extract(pop_raster_2020, shp$geometry,'sum') * adj_factors
-    
-    pop_country_2020 <- WPP2024[WPP2024$ISO3_code == country & WPP2024$Time == 2020,]$PopTotal  * 1e3
-    pop_country_after_2020 <- WPP2024[WPP2024$ISO3_code == country & WPP2024$Time == year,]$PopTotal * 1e3
-    
-    if(length(pop_country_after_2020)>1){ # for years after 2022, there are multiple UN pop estimates and the average of them is used.
-      pop_country_after_2020 = mean(pop_country_after_2020,na.rm=T)
-    }
-    
-    pop_export <- pop_2020/pop_country_2020*pop_country_after_2020
-    
-  } else{
-    pop_raster <- raster::raster(raster_file)
-    pop <- exactextractr::exact_extract(pop_raster, shp$geometry,'sum')
-    
-    ## CA 1 Apr debug NER issue: Ensure single value for pop
-    if (length(pop) > 1) {
-      message("`exact_extract` returned multiple values for a single row")
-    }
-    
-    # Aligning worldpop estimates to the UN population estiamtes at the country level
-    adj_factors <- estimate_adj_factors(
-      country = country, 
-      year = year, 
-      raster_dir = pop_raster_path, 
-      dest_dir = pop_raster_path,
-      country_shp = country_shp
-    )
   pop_export <- pop * adj_factors
-  }
   
   # Check if pop is 0, if so, replace the pop with GHS population 
   if(pop_export == 0){
@@ -168,5 +166,8 @@ get_pop <- function(
     }
   }
   
-  return(pop_export)
+  return(list(
+    estimated_pop = pop_export,
+    adj_factors = adj_factors # also return the adjusting factors in the output
+  ))
 }
diff --git a/R/identify_epidemic_start.R b/R/identify_epidemic_start.R
index 522b1d3..7008d12 100644
--- a/R/identify_epidemic_start.R
+++ b/R/identify_epidemic_start.R
@@ -5,28 +5,91 @@
 #' @description per location period per start_day group: weekly cholera incidence
 identify_epidemic_start <- function(
     outbreak_data = outbreak_data, #consecutive weekly outbreak data
-    minimum_consecutive_reports=3
+    outbreak_start_definition = c("consecutive","dual_window"),
+    min_weeks_above = 2,
+    require_increasing_trend = FALSE,
+    window_weeks = 3,
+    use_cumulative_trigger=TRUE,
+    cumulative_trigger_type=c("cumulative_case_threshold","cumulative_case_threshold_and_min_cases",'cumulative_case_threshold_and_nonzero_weeks'),
+    cumulative_windows = 3,
+    cumulative_case_threshold_ratio = 1.5,
+    cumulative_min_cases = NULL,
+    nonzero_windows =NULL
 ) {
-  #identify consecutively increasing weeks
-  outbreak_data$diff_sCh=c(diff(outbreak_data$sCh),0)
-  outbreak_data_with_epidemic_start<-data.frame()
+  outbreak_data_with_epistart<- outbreak_data
+  outbreak_data_with_epistart$epidemic_start <- FALSE
 
-    outbreak_data_with_epistart<- outbreak_data
-    outbreak_data_with_epistart$consecutive_increase <- FALSE
-    outbreak_data_with_epistart$epidemic_start <- FALSE
-
-    if(nrow(outbreak_data_with_epistart[which(outbreak_data_with_epistart$risk=='high'),])>=minimum_consecutive_reports){
-      for (idx in 1:(nrow(outbreak_data_with_epistart)-minimum_consecutive_reports+1)) {
-        if(all(outbreak_data_with_epistart$diff_sCh[idx:(idx+minimum_consecutive_reports-2)]>0)){
-          outbreak_data_with_epistart$consecutive_increase[idx:(idx+minimum_consecutive_reports-1)]=TRUE
-          #the first week has to be exceeding the outbreak threshold
-          if(outbreak_data_with_epistart$risk[idx] == "high"){
-            outbreak_data_with_epistart$epidemic_start[idx] = TRUE
+  if(outbreak_start_definition == "consecutive"){
+    if(nrow(outbreak_data_with_epistart) < min_weeks_above){
+      return(outbreak_data_with_epistart)
+    }
+    
+    if(nrow(outbreak_data_with_epistart[which(outbreak_data_with_epistart$risk=='high'),])>=min_weeks_above){
+      if(require_increasing_trend){
+        outbreak_data$diff_sCh=c(diff(outbreak_data_with_epistart$sCh),0)
+        outbreak_data_with_epistart$consecutive_increase=FALSE
+        
+        for (idx in 1:(nrow(outbreak_data_with_epistart)-min_weeks_above+1)) {
+          if(all(outbreak_data_with_epistart$diff_sCh[idx:(idx+min_weeks_above-2)]>0)){
+            outbreak_data_with_epistart$consecutive_increase[idx:(idx+min_weeks_above-1)]=TRUE
+            #the first week has to be exceeding the outbreak threshold
+            if(outbreak_data_with_epistart$risk[idx] == "high"){
+              outbreak_data_with_epistart$epidemic_start[idx] = TRUE
+            }
           }
         }
+      } else {
+         for (idx in 1:(nrow(outbreak_data_with_epistart)-min_weeks_above+1)) {
+           if(all(outbreak_data_with_epistart$risk[idx:(idx+min_weeks_above-1)] == "high")){
+             outbreak_data_with_epistart$epidemic_start[idx] = TRUE
+           }
+        }
       }
-    }
-    outbreak_data_with_epidemic_start <- rbind(outbreak_data_with_epidemic_start, outbreak_data_with_epistart)
+  }
+} else if(outbreak_start_definition=="dual_window"){
+  
+  if (sum(outbreak_data_with_epistart$risk=="high")<min_weeks_above) {
+    return(outbreak_data_with_epistart)
+  }
   
-  return(outbreak_data_with_epidemic_start)
+  outbreak_data_with_epistart <- outbreak_data_with_epistart %>%
+    mutate(
+      risk_num= as.integer(risk == 'high'),
+      sCh_cum_thresh = threshold * pop * cumulative_case_threshold_ratio,
+      sCh_pos= as.integer(sCh > 0)
+    ) %>%
+    mutate(
+      d1   = slider::slide_dbl(risk_num, sum, .before = window_weeks - 1,.complete = TRUE), # calculate the number of high risk weeks among window weeks
+      d2   = slider::slide_dbl(sCh,sum, .before = cumulative_windows - 1,.complete = TRUE),# calculate the cumulative cases over the cumulative windows
+      d2_2 = slider::slide_dbl(sCh_pos,  sum, .before = cumulative_windows - 1,.complete = TRUE) # measure the positive case weeks over the cumulative windows
+    ) %>%
+    mutate(
+      d1= replace_na(d1, 0),
+      d2 = replace_na(d2, 0),
+      d2_2= replace_na(d2_2, 0),
+      
+      d2_thresh_met = d2 >= sCh_cum_thresh & sCh_cum_thresh != 0,
+      
+      cumulative_alert = if (!use_cumulative_trigger) {
+        rep(FALSE, n())
+      } else if (cumulative_trigger_type == "cumulative_case_threshold") {
+        d2_thresh_met
+      } else if (cumulative_trigger_type == "cumulative_case_threshold_and_min_cases") {
+        d2_thresh_met & d2 >= cumulative_min_cases
+      } else if (cumulative_trigger_type == "cumulative_case_threshold_and_nonzero_weeks") {
+        d2_thresh_met & d2_2 >= nonzero_windows
+      } else {
+        rep(FALSE, n())
+      },
+      row_id = row_number(),
+      alert = (d1 >= min_weeks_above) | cumulative_alert,
+      d1_start = if_else(d1 >= min_weeks_above, row_id - (window_weeks - 1), NA_integer_),
+      d2_start = if_else(cumulative_alert, row_id - (cumulative_windows - 1), NA_integer_),
+      d1_start = if_else(d1_start %in% row_id[risk == "high"], d1_start, NA_integer_),
+      epidemic_start = row_id %in% na.omit(c(d1_start, d2_start))
+      
+      )
+  }
+
+  return(outbreak_data_with_epistart[,colnames(outbreak_data_with_epistart)%in%c(colnames(outbreak_data),"epidemic_start")])
 }
diff --git a/R/identify_epidemic_tail.R b/R/identify_epidemic_tail.R
index 0542e9e..a130a3a 100644
--- a/R/identify_epidemic_tail.R
+++ b/R/identify_epidemic_tail.R
@@ -8,10 +8,10 @@ identify_epidemic_tail <- function (
 ){
   outbreak_data$epidemic_tail <- FALSE
   #identify tails of outbreaks
-  ## 1. below outbreak threshold. 2. for certain number of consecutive weeks
+  ## 1. below outbreak threshold. 2. for certain number of consecutive weeks. 3. not a outbreak start.
   
   consecutive_vector <- rep(TRUE,nrow(outbreak_data))
-  below_threshold_vector <- outbreak_data$risk == "low"
+  below_threshold_vector <- outbreak_data$risk == "low" & !outbreak_data$epidemic_start
   tail_vector <- consecutive_vector*below_threshold_vector
 
   tail_position <- rle(tail_vector)
diff --git a/R/identify_outbreaks.R b/R/identify_outbreaks.R
index 0f2a64a..41ec621 100644
--- a/R/identify_outbreaks.R
+++ b/R/identify_outbreaks.R
@@ -16,7 +16,18 @@ identify_outbreaks <- function(
     original_data,
     zero_case_assumption = T,
     customized_TL = NULL,
-    customized_TR = NULL
+    customized_TR = NULL,
+    outbreak_start_definition = c("consecutive","dual_window"), 
+    min_weeks_above = 2,
+    require_increasing_trend = FALSE,
+    window_weeks = window_weeks,
+    cumulative_windows = cumulative_windows,
+    cumulative_case_threshold_ratio = cumulative_case_threshold_ratio,
+    cumulative_trigger_type=cumulative_trigger_type,
+    use_cumulative_trigger=use_cumulative_trigger,
+    cumulative_min_cases=cumulative_min_cases,
+    nonzero_windows = nonzero_windows,
+    tail_period =6
     ){
   
   # Identify cholera outbreak thresholds
@@ -38,9 +49,18 @@ identify_outbreaks <- function(
     new_outbreak_by_location=data.frame()
     
     # extract outbreak start and end
-    preoutbreak_by_location_start<-OutbreakExtractR::identify_epidemic_start(outbreak_data = preoutbreak_by_location)
-    
-    preoutbreak_by_location_start_end_washout<-OutbreakExtractR::identify_epidemic_tail(outbreak_data = preoutbreak_by_location_start, tail_period = 6)
+    preoutbreak_by_location_start<-OutbreakExtractR::identify_epidemic_start(outbreak_data = preoutbreak_by_location,
+                                                                             outbreak_start_definition = outbreak_start_definition, 
+                                                                             require_increasing_trend = require_increasing_trend,
+                                                                             min_weeks_above = min_weeks_above,
+                                                                             window_weeks = window_weeks,
+                                                                             use_cumulative_trigger=use_cumulative_trigger,
+                                                                             cumulative_trigger_type=cumulative_trigger_type,
+                                                                             cumulative_windows = cumulative_windows,
+                                                                             cumulative_case_threshold_ratio = cumulative_case_threshold_ratio,
+                                                                             cumulative_min_cases=cumulative_min_cases,
+                                                                             nonzero_windows = nonzero_windows)
+    preoutbreak_by_location_start_end_washout<-OutbreakExtractR::identify_epidemic_tail(outbreak_data = preoutbreak_by_location_start, tail_period = tail_period)
     
     # get the row idx for epidemic start
     preoutbreak_by_location_start_end_washout$row_idx = rownames(preoutbreak_by_location_start_end_washout)
@@ -53,7 +73,7 @@ identify_outbreaks <- function(
         for (idx in seq(length(epidemic_start_row_idx)-1)) {
           
           data_between_epidemic_start = preoutbreak_by_location_start_end_washout[epidemic_start_row_idx[idx]:epidemic_start_row_idx[idx+1],]
-          if(any(data_between_epidemic_start$epidemic_tail)){
+          if(any(data_between_epidemic_start$epidemic_tail) & nrow(data_between_epidemic_start)>=tail_period+2){
             outbreak_end = min(as.numeric(data_between_epidemic_start[data_between_epidemic_start$epidemic_tail,]$row_idx))
             preoutbreak_by_location_start_end_washout[epidemic_start_row_idx[idx]:as.numeric(as.numeric(outbreak_end)+2-1),]$outbreak_number = 
               min(outbreak_number_idx,preoutbreak_by_location_start_end_washout[epidemic_start_row_idx[idx]:as.numeric(as.numeric(outbreak_end)+2-1),]$outbreak_number[preoutbreak_by_location_start_end_washout[epidemic_start_row_idx[idx]:as.numeric(as.numeric(outbreak_end)+2-1),]$outbreak_number>0])