###Functions###
# Exponential decay function for distance
decay_function  <-  function(distance, mean_distance) {
  return(exp(-distance / mean_distance))
}


#function to calculate dispersal array#
#resource_array is the array with resource densities (either flowers for ill fed, or aphids for well-fed)
#distance array is the distance array (which can be changed to donut or something else)
#resource_local is a parameter, determining if individuals search locally or globally
#distance_local is a parameter, determining the local distance
#distance_global is a parameter, determining the global distance

calculate_dispersal_array <- function(resource_array, distance_array, resource_local, distance_local, distance_global) {
  # Initialize the dispersal array
  dispersal_array <- array(0, dim = dim(distance_array))
  
  # Determine the mean distance based on the resource levels
  mean_distance_array <- ifelse(resource_array > resource_local, distance_local, distance_global)

  #expand to 4d -> necessary because we want to have it similar to the distance array
  mean_distance_expanded <- array(mean_distance_array, dim = dim(distance_array))

  # Apply the decay function using vectorized operations
  dispersal_array <- exp(-distance_array / mean_distance_expanded) 
  dispersal_array[distance_NA] = 0 #where distance is either zero or too far
  #Normalize (total_dispersal sums over the last two dimension (l,k) = destination)
  #for each cell in i,j. alternatively, it does sum(dispersal_array[i,j, ,]) in a for loop over both i and j
  total_dispersal <- apply(dispersal_array, c(1, 2), sum) #=2 dim dim nrow ncol
  
  # Avoid division by zero by setting all zeros to 1 (cause x/1 = x)
  total_dispersal[total_dispersal == 0] <- 1
  
  # Normalize by dividing each slice by the corresponding total dispersal
  #here, the dispersal_array is divided by the total_dispersal,
  #c(1,2) tells R what to match, so it will make sure dispersal_array[i,j,k,l]
  #gets divided by total_dispersal[i,j] (and not total_dispersal[k,l] or [j,k] etc)
  dispersal_array <- sweep(dispersal_array, c(1, 2), total_dispersal, "/")
  
  return(dispersal_array)
}

# Function to apply dispersal for one time step 
#Ad_new = Array of adults (either well or ill fed)
#dispersal_array = dispersal probabilities
#stay_fraction = fraction not dispersing.
apply_dispersal <- function(Ad_new, dispersal_array, stay_fraction) {
  # Calculate the outflows from each cell
  outflows <- Ad_new * (1 - stay_fraction) * dt
  
  # Update the Ad_new matrix to account for the remaining population
  Ad_new <- Ad_new * (1 - dt + dt * stay_fraction)
  
  # Calculate the inflows by applying the dispersal array
  #slice gets one layer of the dispersal_array, in this case the destination [k, l]
  #than, it calculates how many individuals from the source [i,j] (=outflows)
  #will move to that destination
  #it would be similar as a for loop for each destination cell [k, l] compute 
  #how many individuals will move from the source cell [i,j]
  inflows <- apply(dispersal_array, c(3, 4), function(slice) {
    sum(outflows * slice)
  })
  
  # Add inflows to the updated Ad_new matrix
  Ad_new <- Ad_new + inflows
  
  return(Ad_new)
}



#HOVERFLY LIFE HISTORY
#Type 2 functional response hoverflies to aphids
f <- function(N){
  fm * (N / (N + Nf)) * tc * dt
}

#Juvenile developmental rate hoverflies
e <- function(N){
  em * (N/(N + Ne) + phi) * tc * dt
}

#Juvenile mortality rate hoverflies
muj <- function(N){
  mum * ((N + Ne * 0.5)/(N + epsilon)) * tc * dt
}

#Reproduction rate hoverflies (lower in H1)
G1 <- function(N){
  G1m * (N / ( N + Ng)) * tc * dt
}
G2 <- function(N){
  G2m * (N/(N + Ng)) * tc * dt
}



#plotting
is_multiple <- function(t, step, tol){
 abs(t / step - round(t / step)) < tol
}