Simulate prey fields with a spectral power law of spatial frequency to the –1.5

prey_field.R

library(raster)
library(scales)
library(tidyverse)

powpow <- function(n, a, b) {
  # based on answers at: https://dsp.stackexchange.com/questions/47640/generating-a-timeseries-with-an-arbitrary-power-spectrum
  
  # Thank you Sam!
  
  if (n %% 2 != 0)
    stop("n must be even")
  
  n0 <- n / 2
  
  # the frequencies for the FFT of a time serie with n samples
  freqs <- 0:(n0 - 1) / n
  freqs <- c(freqs, -n0 / n, -rev(freqs[2:n0]))
  
  pow <- a * abs(freqs)^b # the power-law power spectrum we want
  pow[1] <- 1 # make the 0-frequency term finite
  
  phase <- runif(n0, 0, 2 * pi)
  phase <- c(phase, -phase)
  
  xft <- sqrt(pow) * exp(phase * 1i)
  Re(fft(xft, inverse = TRUE))
}

# simulate a 2d prey field
sim_prey <- function(n, a, b, r_sz, xmn, xmx, ymn, ymx) {
  prey <- tibble(x = powpow(n, a, b),
                 y = powpow(n, a, b)) %>% 
    mutate(x = rescale(x, to = c(xmn, xmx)),
           y = rescale(y, to = c(ymn, ymx)))
  # rasterize to get density
  r <- raster(matrix(0, r_sz, r_sz), 
              xmn = xmn, xmx = xmx,
              ymn = ymn, ymx = ymx)
  counts <- table(cellFromXY(r, as.matrix(prey)))
  r[as.numeric(names(counts))] <- counts
  r
}
prey_fields <- map(1:9, ~ sim_prey(2^22, 0.05, -1.5, 2.5e2, 0, 1e2, 0, 1e2))
par(mfrow = c(3, 3))
for (i in seq_along(prey_fields)) {
  plot(prey_fields[[i]], main = paste("Prey field", i))
}
par(mfrow = c(1, 1))