diff --git a/R/recovery.FUN.R b/R/recovery.FUN.R
index 60f9200de45e50283797353bfe7c281f22fc72d3..eadc963886a13d83c25b0c90ccc913e76249834f 100644
--- a/R/recovery.FUN.R
+++ b/R/recovery.FUN.R
@@ -1,3 +1,51 @@
+# internal function --------------------------------------------------------------
+update_status <- function(rate, rand, active, status, label, state, at, prog, expTime = NULL, prog.dist.scale = NULL, prog.dist.shape = NULL){
+  
+  smp_sz <- 0
+  at_idx <- NULL
+  
+  idsElig <- which(active == 1 & status == label)
+  nElig <- length(idsElig)
+  
+  if (nElig > 0) {
+    gElig <- rep(1, nElig)
+    
+    if (rand) {
+      vecProg <- which(stats::rbinom(nElig, 1, rate) == 1)
+      if (length(vecProg) > 0) {
+        idsProg <- idsElig[vecProg]
+        smp_sz <- length(idsProg)
+        status[idsProg] <- state
+        at_idx <- idsProg
+      }
+    }else{
+      do_sample <- TRUE
+      if(!prog){
+        smp_sz <- min(round(sum(rates[gElig == 1])), sum(gElig == 1))
+        smp_prob <- NULL
+      }else{
+        vecTimeSinceExp <- at - expTime[idsElig]
+        gammaRatesElig <- stats::pweibull(vecTimeSinceExp, prog.dist.shape, scale=prog.dist.scale)
+        smp_sz <- round(sum(gammaRatesElig[gElig == 1], na.rm=TRUE)) 
+        smp_prob <- gammaRatesElig[gElig == 1]
+        if(smp_sz <= 0) do_sample <- FALSE
+      }
+      
+      if(do_sample){
+        ids2bProg <- EpiModel::ssample(idsElig[gElig == 1], smp_sz, prob = smp_prob)
+        status[ids2bProg] <- state
+        at_idx <- ids2bProg
+      }
+      
+      
+    }
+    
+  }
+  
+  list(smp_sz, at_idx, status)
+}
+
+
 #' Progress icm
 #'
 #' Function to get progress of icms
@@ -21,56 +69,6 @@ recovery.FUN <- function(dat, at, seed = NULL) {
     return(dat)
   }
   
-  
-  # internal function --------------------------------------------------------------
-  update_status <- function(rate, rand, active, status, label, state, at, prog, expTime = NULL, prog.dist.scale = NULL, prog.dist.shape = NULL){
-    
-    smp_sz <- 0
-    at_idx <- NULL
-    
-    idsElig <- which(active == 1 & status == label)
-    nElig <- length(idsElig)
-    
-    if (nElig > 0) {
-      gElig <- rep(1, nElig)
-      
-      if (rand) {
-        vecProg <- which(stats::rbinom(nElig, 1, rate) == 1)
-        if (length(vecProg) > 0) {
-          idsProg <- idsElig[vecProg]
-          smp_sz <- length(idsProg)
-          status[idsProg] <<- state
-          at_idx <- idsProg
-        }
-      }else{
-        do_sample <- TRUE
-        if(!prog){
-          smp_sz <- min(round(sum(rates[gElig == 1])), sum(gElig == 1))
-          smp_prob <- NULL
-        }else{
-          vecTimeSinceExp <- at - expTime[idsElig]
-          gammaRatesElig <- stats::pweibull(vecTimeSinceExp, prog.dist.shape, scale=prog.dist.scale)
-          smp_sz <- round(sum(gammaRatesElig[gElig == 1], na.rm=TRUE)) 
-          smp_prob <- gammaRatesElig[gElig == 1]
-          if(smp_sz <= 0) do_sample <- FALSE
-        }
-        
-        if(do_sample){
-          ids2bProg <- EpiModel::ssample(idsElig[gElig == 1], smp_sz, prob = smp_prob)
-          status[ids2bProg] <<- state
-          at_idx <- ids2bProg
-        }
-        
-        
-      }
-      
-      dat$attr$status <<- status
-    }
-    
-    list(smp_sz, at_idx)
-  }
-  
-  
   # Variables ---------------------------------------------------------------
   active <- dat$attr$active
   status <- dat$attr$status
@@ -96,7 +94,8 @@ recovery.FUN <- function(dat, at, seed = NULL) {
   
   nProg <- res[[1]]
   if(!is.null(res[[2]])) dat$attr$infTime[res[[2]]] <- at
-  
+  status <- res[[3]]
+  dat$attr$status <- status
   
   if (type %in% c("SEIQHR", "SEIQHRF")) {
     
@@ -112,7 +111,8 @@ recovery.FUN <- function(dat, at, seed = NULL) {
                          state = "q", at, prog = FALSE)
     nQuar <- res[[1]]
     if(!is.null(res[[2]])) dat$attr$quarTime[res[[2]]] <- at
-    
+    status <- res[[3]]
+    dat$attr$status <- status
     
     # ----- need to be hospitalised -------
     hosp.rand <- dat$control$hosp.rand