library(readxl)  # Load package to read Excel files

###Initialization of all the vectors etc####

# DEFINING HABITATS H <- 1-4 in the LANDSCAPE####

#Load new landscape in script Run_spatialmodel
# Convert loaded landscape to matrix (since data frames behave differently)
H <- as.matrix(H_loaded)

# Ensure it's numeric
H <- apply(H, c(1,2), as.numeric)

#Defining nrow and ncol from habitat matrix
dimsL=dim(H)
nrowL=dimsL[1]
ncolL=dimsL[2]
ncells=length(H)

#Calculate proportion of each habitat in the landscape
cell_counts <- table(H)

ncellsH0=as.numeric(cell_counts["0"])
ncellsH1=as.numeric(cell_counts["1"])
ncellsH2=as.numeric(cell_counts["2"])
ncellsH3=as.numeric(cell_counts["3"])
ncellsH4=as.numeric(cell_counts["4"])



propH0=ncellsH0/ncells
propH1=ncellsH1/ncells
propH2=ncellsH2/ncells
propH3=ncellsH3/ncells
propH4=ncellsH4/ncells



#Set up habitats#
no_habitat <- (H==0)
woody_habitat <- (H == 1)
early_crop <- (H == 2)
late_crop <- (H == 3)
flower_margin <- (H == 4)




# SETTING UP ARRAYS representing populations in SPACE ##########

#Aphids
N_prev <- array(0, dim = c(nrow = nrowL, ncol = ncolL))
N_new <- array(0, dim = c(nrow = nrowL, ncol = ncolL))

#Juvenile hoverflies
Juv_prev <- array(0, dim = c(nrow = nrowL, ncol = ncolL))
Juv_new <- array(0, dim = c(nrow = nrowL, ncol = ncolL))

#Ill-fed adult hoverflies
Ad_ill_prev <- array(0, dim = c(nrow = nrowL, ncol = ncolL))
Ad_ill_new <- array(0, dim = c(nrow = nrowL, ncol = ncolL))

#Well-fed adult hoverflies
Ad_well_prev <- array(0, dim = c(nrow = nrowL, ncol = ncolL))
Ad_well_new <- array(0, dim = c(nrow = nrowL, ncol = ncolL))

#Hibernating adult hoverflies
Ad_hib_prev <- array(0, dim = c(nrow = nrowL, ncol = ncolL))
Ad_hib_new <- array(0, dim = c(nrow = nrowL, ncol = ncolL))

#Flowers (HtB: In crop fields right? )
B <- array(B23l, dim = c(nrow = nrowL, ncol = ncolL))

# Proportion of hoverflies that do not disperse from cell
zB  <-  array(0, dim  =  c(nrowL, ncolL))
zN  <-  array(0, dim  =  c(nrowL, ncolL))


# Dispersal arrays to store dispersal probabilities, separate for ill and well-fed adults #*
dispersal_array_ill <- array(0, dim  =  c(nrowL, ncolL, nrowL, ncolL))  #*
dispersal_array_well <- array(0, dim  =  c(nrowL, ncolL, nrowL, ncolL)) #*

#SPACExTIME arrays to store simulation results####
N_array_days       <- array(dim = c(nrow = nrowL, ncol = ncolL, num_days / dt))
Juv_array_days     <- array(dim = c(nrow = nrowL, ncol = ncolL, num_days / dt))
Ad_ill_array_days  <- array(dim = c(nrow = nrowL, ncol = ncolL, num_days))
Ad_well_array_days <- array(dim = c(nrow = nrowL, ncol = ncolL, num_days))
Ad_hib_array_days  <- array(dim = c(nrow = nrowL, ncol = ncolL, num_days))
B_array_days       <- array(dim = c(nrow = nrowL, ncol = ncolL, num_days))

#VECTORS to store means per habitat through time
sumN1  <-  0
sumN2  <-  0
sumN3  <-  0
sumN5  <-  0  
sumJuv1  <-  0
sumJuv2  <-  0
sumJuv3  <-  0
sumJuv5  <-  0
N1_mean  <-  array(dim = num_days)
N2_mean  <-  array(dim = num_days)
N3_mean  <-  array(dim = num_days)
N5_mean  <-  array(dim = num_days)
Juv1_mean  <-  array(dim = num_days)
Juv2_mean  <-  array(dim = num_days)
Juv3_mean  <-  array(dim = num_days)
Juv5_mean  <-  array(dim = num_days)
Ad_ill_mean  <-  array(dim = num_days)
Ad_well_mean  <-  array(dim = num_days)
Ad_hib_mean  <-  array(dim = num_days)


###Calculate the distance array (non-wrapped)####
# distance_array <- array(0, dim = c(nrowL, ncolL, nrowL, ncolL))
# for (i in 1:nrowL) {
#   for (j in 1:ncolL) {
#     for (k in 1:nrowL) {
#       for (l in 1:ncolL) {
#         distance_array[i, j, k, l] <- sqrt((i - k)^2 + (j - l)^2)
#       }
#     }
#   }
# }


#alternative distance array (wrapped)####
distance_array <- array(0, dim = c(nrowL, ncolL, nrowL, ncolL))
for (i in 1:nrowL) {
  for (j in 1:ncolL) {
    for (k in 1:nrowL) {
      for (l in 1:ncolL) {
        # Calculate the shortest distance in each dimension, accounting for wrapping
        dy <- min(abs(i - k), nrowL - abs(i - k))
        dx <- min(abs(j - l), ncolL - abs(j - l))
        # Calculate Euclidean distance using the wrapped distances
        distance_array[i, j, k, l] <- sqrt(dx^2 + dy^2)
      }
    }
  }
}
distance_NA <- (distance_array > max_distance | distance_array == 0) #where it is not relevant




###Initialization of the grid######
N_prev <- pmax(N_new, MinVal)
Juv_prev <- pmax(Juv_new, MinVal)
Ad_ill_prev <- pmax(Ad_ill_new, MinVal)
Ad_well_prev <- pmax(Ad_well_new, MinVal)
Ad_hib_prev <- pmax(Ad_hib_new, MinVal)

#Only for first year at t = 0
N_prev[woody_habitat]  <-  N1i
Ad_hib_prev[woody_habitat] <- Pi/propH1



#Initialisation for visualisation
N_plot_list <- list()
Ad_well_plot_list <- list()
Ad_ill_plot_list <- list()
hib_plot_list <- list()
Juv_plot_list <- list()