#Vectorized code with correct dispersal array and correct order of events#
#setwd("~/Documents/PhD-project Laura/Spatial Model/Code/")
#rm(list=ls())

#correct the dt things#


###Run model#####
total_time = 0
tol <- dt/100
for(y in 1:year) {
  print(paste("year:", y))
  index_dt = 1
  index_t = 1
  for(t in seq(from = 1, to = num_days, by = dt)) {
    #print(t)
    #print(max(Ad_ill_new))
    ####Functions#####
    # Calculate the current day for output
    current_day <- floor(t)
    
    #FUNCTIONS
    #TEMPERATURE change within a year
    Temp <- Ta*cos((t - td) * 2 * pi / 365) + Tc
    #Temperature correction for rates
    tc <- (Temp-4)/(22-4)
    
    #FLORAL RESOURCES  #*
    #Functions flower food availability per habitat
    B1 <- B1l + B1m * dweibull((t-T1s)/T1d, k1, 1) 
    B4 <- B4l + B4m * dweibull((t-T4s)/T4d, ks, 1) #annual
    #B4 <- B4l + B4m * dweibull((t-T4s)/T4d, ks, 1) #perennial
    
    ###Update flowers that vary over time####
    B[woody_habitat] <- B1
    B[flower_margin] <- B4
    
    #APHID DYNAMICS 
    # Aphid resource availability and carrying capacity (Weibull, reverse Weibull)
    n1 <- n1s * dweibull((t-A1s)/A1d,ks,1) + n1a*dweibull((t-A1as)/A1ad,ks,1)
    K2 <- Km * dweibull((A2f-t)/A2d,ks,1)+phi
    K3 <- Km * dweibull((A3f-t)/A3d,ks,1)+phi
    n2 <- ifelse(K2>0.01*Km, 1, 0)
    n3 <- ifelse(K3>0.01*Km, 1, 0)
    
    #Aphid population growth rate
    r1 <- n1 * rN1 * tc * dt
    r2 <- n2 * rN2 * tc * dt
    r3 <- n3 * rN3 * tc * dt
    
    #time delay representing inactive stages
    tau <- tauh/tc  #tau niet delen door dt, want t is in eenheden van 1 dag  
    #needs to be rounded because array only every day. 
    
    #HOVERFLY HIBERNATION
    #Release from hibernation (around day H0=30) h0
    h0 <- dnorm(t, H0, sigma0) * dt
    
    #Induction into hibernation (around day H1=183) h1  
    h1 <- dnorm(t, H1, sigma1) * dt
    
    
    ####POPULATION DYNAMICS####
    
    ###Aphids####
    
    #aphids infest in other habitats 
    if(abs(t - T2i) < tol) {
      N_prev[early_crop] <- N2i
    }
    
    if(abs(t - T3i)< tol) {
      N_prev[late_crop] <- N3i
    }
    
    
    #Woody habitat
    #note that the function f() is already taking care of temperature (tc) and timestep (dt)

    if(t<=Wt1){
    ###Aphids
    N_susceptible <- pmax(N_prev[woody_habitat] - Nr, 0) #individuals that can die 
    N_save <- pmin(Nr, N_prev[woody_habitat]) #individuals that are hidden and cannot die.
    
      N_new[woody_habitat] <- 
      N_prev[woody_habitat] * r1 * (1-N_prev[woody_habitat]/K1) + #growth
    (N_susceptible - pmin(Juv_prev[woody_habitat]*f(N_susceptible), N_susceptible))*(1-mu1*tc*dt) + #mortality (predation followed by background)
        N_save #hidden aphids. 
        
    

    
    ###Juveniles hoverflies####
    Juv_new[woody_habitat] <- Juv_prev[woody_habitat] * (1-muj(N_susceptible)) * 
                              (1-e(N_susceptible)) +
                              G1(N_susceptible) * Ad_well_prev[woody_habitat] * (1-muB * tc * dt)
      if(min(Juv_new[woody_habitat]) < 0) {print(paste("Stop at year ", y, " and day ", t))} 
      }
    
    #winter mortality
    if(abs(t - Wt1) < tol) {
      N_new[woody_habitat] <- WN * N_new[woody_habitat]
      Juv_new[woody_habitat] <- 0
    } 


    #Early crop
    if(T2i <= t & t <= Wt2){
    #Aphids H2
    N_new[early_crop] <- N_prev[early_crop] * r2 * (1-N_prev[early_crop]/K2) + 
    (N_prev[early_crop] - pmin(Juv_prev[early_crop]*f(N_prev[early_crop]), N_prev[early_crop]))*(1-mu2*tc*dt)
     
    
    Juv_new[early_crop] <- Juv_prev[early_crop] * (1-muj(N_prev[early_crop])) * 
                              (1-e(N_prev[early_crop])) +
                              G2(N_prev[early_crop]) * Ad_well_prev[early_crop] * (1-muB * tc * dt)
    }


    #Harvest mortality
    if(abs(t - Wt2) < tol) {
      Juv_new[early_crop] <- 0
      N_new[early_crop] <- 0
    }

    
    #Late crop
    if(T3i <= t & t <= Wt3){
    #Aphids H3
    N_new[late_crop] <- N_prev[late_crop] * r3 * (1-N_prev[late_crop]/K3) + 
    (N_prev[late_crop] - pmin(Juv_prev[late_crop]*f(N_prev[late_crop]),N_prev[late_crop]))*(1-mu3*tc*dt)
    
    

    #Juveniles
    Juv_new[late_crop] <- Juv_prev[late_crop] * (1-muj(N_prev[late_crop])) * 
                              (1-e(N_prev[late_crop])) +
                              G2(N_prev[late_crop]) * Ad_well_prev[late_crop] * (1-muB * tc * dt)
      }

    if(abs(t - Wt3) < tol) { #mortality
      Juv_new[late_crop] <- 0
      N_new[late_crop] <- 0
    }


    
    ###Adult hoverflies new####
    #1) Survival
    Ad_ill_new <- Ad_ill_prev * (1-(muB+muS) * tc * dt)        #survival of ill-fed
    Ad_well_new <- Ad_well_prev * (1-muB * tc * dt)                 #survival of well-fed
    
    #2) Reproduction of adults into juveniles (no change in adults)
    #calculated under juveniles
    
    #3) change of feeding status
    Ad_well_save = Ad_well_new #necessary to save in order to not change status multiple times
    Ad_ill_save = Ad_ill_new
    
    Ad_ill_new  <- Ad_ill_new  + (Ad_well_save * a * tc * dt) - (Ad_ill_save * b * B * tc * dt)
    Ad_well_new <- Ad_well_new + (Ad_ill_save * b * B * tc * dt) - (Ad_well_save * a * tc *dt)
    
    
    #4 Juveniles develop into ill-fed adults
    tau_rounded <- round(tau / dt) * dt
    index <- round((t - tau_rounded - 1) / dt)
    if(t >= (tau + 2)) {
      susceptible <- pmax(N_array_days[, , index] - Nr, 0)
      Ad_ill_new <- Ad_ill_new + Juv_array_days[, , index] *
        (1-muj(susceptible)) *
        e(susceptible)
    }
    
    #5) End of hibernation 
    Ad_hib_new[woody_habitat] <- Ad_hib_prev[woody_habitat]  * (1 - h0) #this could be done without woody_habitat since only in woody habitat anyway
    Ad_ill_new[woody_habitat] <- Ad_ill_new[woody_habitat] + Ad_hib_prev[woody_habitat] * h0 #adult hoverflies leaving hibernation
    
    #6) Start of hibernation 
    Ad_hib_new[woody_habitat] <- Ad_hib_new[woody_habitat] * (1-mu0*dt*tc) + (Ad_well_new[woody_habitat])*h1
    Ad_well_new[woody_habitat] <- Ad_well_new[woody_habitat] * (1 - h1) #well-fed going into hibernation around day h1=183
    
   #7) winter mortality of adults
    #Eerst was Wt1 dag 209, maar ik heb er nu de laatste dag 210 van gemaakt
    if(abs(t - Wt1) < tol) {
      Ad_hib_new <- Ad_hib_new * WP
      Ad_well_new <- 0
      Ad_ill_new <- 0
    } 
    
    ####Dispersal of adults####
    #DISPERSAL####
    # Proportion of (ill-fed) hoverflies that stay in each cell due to flowers (zB)
    zB <- B / (Bd + B)  # Using Michaelis-Menten equation #*
    # Proportion of (well-fed) hoverflies that stay in each cell due to aphids (zN)
    zN <- N_new / (Nd + N_new)  # Using Michaelis-Menten equation
    
    #APPLYING DISPERSAL TO Populations
    #update dispersal not always necessary, maybe only every day and not every dt day
    #dispersal itself should happen each time step I think 
    if (round(t) %in% Update_dispersal && abs(t - round(t)) < tol) {
      #print(paste('At time', t, "dispersal arrays are updated"))
      dispersal_array_ill <- calculate_dispersal_array(B, distance_array, B_localsearch, distance_local, distance_global)
      dispersal_array_well <- calculate_dispersal_array(N_new, distance_array, N_localsearch, distance_local, distance_global)
    }
      Ad_ill_new <- apply_dispersal(Ad_ill_new, dispersal_array_ill, zB)  #*
      Ad_well_new <- apply_dispersal(Ad_well_new, dispersal_array_well, zN)  #*
     
      #print(paste("Arrays at time", t, " and index ", index))
      N_array_days[,,index_dt] <- N_new 
      Juv_array_days[,,index_dt] <- Juv_new
      #print(Juv_array_days[1,1,index])
      index_dt = index_dt + 1
    
    #DATA STORAGE#### 
    #Collect some data#
    
    sumN1 <- sumN1 + sum(N_new[woody_habitat])
    sumJuv1 <- sumJuv1 + sum(Juv_new[woody_habitat])
    sumN2 <- sumN2 + sum(N_new[early_crop])
    sumJuv2 <- sumJuv2 + sum(Juv_new[early_crop])
    sumN3 <- sumN3 + sum(N_new[late_crop])
    sumJuv3 <- sumJuv3 + sum(Juv_new[late_crop])

    if (abs(t - current_day) < tol) {
      total_time = floor(t) + (y - 1) * num_days
      #Creating array for all cells on all days of 1 year
      #N_array_days[,,t] <- N_new 
      #Juv_array_days[,,t] <- Juv_new
      Ad_ill_array_days[,,index_t] <- Ad_ill_new
      Ad_well_array_days[,,index_t] <- Ad_well_new
      Ad_hib_array_days[,,index_t] <- Ad_hib_new
      B_array_days[,,index_t] <- B
      
      # Creating a list with population averages for all days in 1 year
      N1_mean[[total_time]] <- sumN1*dt/ncellsH1
      sumN1<-0
      N2_mean[[total_time]] <- sumN2*dt/ncellsH2
      sumN2<-0
      N3_mean[[total_time]] <- sumN3*dt/ncellsH3
      sumN3<-0
      Juv1_mean[[total_time]] <- sumJuv1*dt/ncellsH1
      sumJuv1<-0
      Juv2_mean[[total_time]] <- sumJuv2*dt/ncellsH2
      sumJuv2<-0  
      Juv3_mean[[total_time]] <- sumJuv3*dt/ncellsH3
      sumJuv3<-0

      
      Ad_ill_mean[[total_time]] <- mean(Ad_ill_new)
      Ad_well_mean[[total_time]] <- mean(Ad_well_new)
      Ad_hib_mean[[total_time]] <- mean(Ad_hib_new)
      index_t = index_t + 1
      ###Visualise data####
      #PLOTTING distributions
      if (y %in% plotyears & t %% 5 == 0 & resplot == T) {
        
        
        
        
        #Aphids
        # Convert matrix to a dataframe for ggplot
        N_df <- melt(t(N_new))
        colnames(N_df) <- c("Row", "Column", "Value")
        
        Aphid_title <- paste0("Aphid density (Year: ", y, ", Day: ", t, ")")
        
        # Plot using ggplot with a gradient color scale
        N_plot = ggplot(N_df, aes(x = Column, y = Row, fill = Value)) +
          geom_tile() +
          scale_fill_gradientn(colors = c("blue", "green", "yellow", "red"),
                               name = "Density",
                               trans = "log",  # Apply log scale transformation
                               limits = c(0.01, 1000),  # Avoid log(0), so set min to 1
                               breaks = c(0.1, 1, 10, 100, 1000),  # Log scale breaks
                               labels = scales::comma) +
          
          theme_minimal() +
          labs(title = Aphid_title, fill = "Value") +
          scale_y_reverse() +  # Flip y-axis so (1,1) is top-left
          theme(axis.text.x = element_blank(),
                axis.text.y = element_blank(),
                axis.ticks = element_blank(),
                axis.title = element_blank(),
                panel.grid = element_blank())
        
        # plot(N_plot)
        
        # Store plot in a list with dynamic naming
        N_plot_list[[paste0("plot_day", t)]] <- N_plot
        
        
        ##Well-fed
        Ad_well_df <- melt(t(Ad_well_new))
        colnames(Ad_well_df) <- c("Row", "Column", "Value")
        
        Ad_well_title <- paste0("Well-fed hoverfly density (Year: ", y, ", Day: ", t, ")")
        
        # Plot using ggplot with a gradient color scale
        Ad_well_plot = ggplot(Ad_well_df, aes(x = Column, y = Row, fill = Value)) +
          geom_tile() +
          scale_fill_gradientn(colors = c("blue", "green", "yellow", "red"),
                               name = "Density",
                               trans = "log",  # Apply log scale transformation
                               limits = c(0.001, 10),  # Avoid log(0), so set min to 1
                               breaks = c(0.001, 0.01, 0.1, 1, 3),  # Log scale breaks
                               labels = scales::comma) +
          
          theme_minimal() +
          labs(title = Ad_well_title, fill = "Value") +
          scale_y_reverse() +  # Flip y-axis so (1,1) is top-left
          theme(axis.text.x = element_blank(),
                axis.text.y = element_blank(),
                axis.ticks = element_blank(),
                axis.title = element_blank(),
                panel.grid = element_blank())
        
        # plot(Ad_well_plot)
        
        Ad_well_plot_list[[paste0("plot_day", t)]] <- Ad_well_plot
        
        ##Ill-fed
        Ad_ill_df <- melt(t(Ad_ill_new))
        colnames(Ad_ill_df) <- c("Row", "Column", "Value")
        
        Ad_ill_title <- paste0("ill-fed hoverfly density (Year: ", y, ", Day: ", t, ")")
        
        # Plot using ggplot with a gradient color scale
        Ad_ill_plot = ggplot(Ad_ill_df, aes(x = Column, y = Row, fill = Value)) +
          geom_tile() +
          scale_fill_gradientn(colors = c("blue", "green", "yellow", "red"),
                               name = "Density",
                               trans = "log",  # Apply log scale transformation
                               limits = c(0.001, 10),  # Avoid log(0), so set min to 1
                               breaks = c(0.001, 0.01, 0.1, 1, 3),  # Log scale breaks
                               labels = scales::comma) +
          
          theme_minimal() +
          labs(title = Ad_ill_title, fill = "Value") +
          scale_y_reverse() +  # Flip y-axis so (1,1) is top-left
          theme(axis.text.x = element_blank(),
                axis.text.y = element_blank(),
                axis.ticks = element_blank(),
                axis.title = element_blank(),
                panel.grid = element_blank())
        
        # plot(Ad_ill_plot)
        
        Ad_ill_plot_list[[paste0("plot_day", t)]] <- Ad_ill_plot
  
      }
      
      
      
      
      
    } #end t<-current-day loop
    
    ####End, set new values to old values for the next time step####
    #and make sure there are no zero's and negatives
    N_new <- pmax(N_new, MinVal)
    N_prev <- N_new
    Juv_new <- pmax(Juv_new, MinVal)
    Juv_prev <- Juv_new
    Ad_ill_new <- pmax(Ad_ill_new, MinVal)
    Ad_ill_prev <- Ad_ill_new
    Ad_well_new <- pmax(Ad_well_new, MinVal)
    Ad_well_prev <- Ad_well_new
    Ad_hib_new  <- pmax(Ad_hib_new, MinVal)
    Ad_hib_prev <- Ad_hib_new
    
  } # end TIME-loop
}

###Save data####
time = seq(1, num_days * year)


NewData <- data.frame(time = time, 
                      well = Ad_well_mean[time], 
                      ill = Ad_ill_mean[time], 
                      hib = Ad_hib_mean[time],
                      Aphids_N1 = N1_mean[time], 
                      Aphids_N2 = N2_mean[time],
                      Aphids_N3 = N3_mean[time],
                      Juv_N1 = Juv1_mean[time],
                      Juv_N2 = Juv2_mean[time],
                      Juv_N3 = Juv3_mean[time]
                      )
NewData$Day = (NewData$time - 1) %% num_days + 1
NewData$Year = (NewData$time -1) %/% num_days + 1


