lambda0 <- 20
c <- 2

f <- function(x1, x2, c){
  a <- sqrt(x1 * x1 + x2 * x2)
  if(a * c <= 1){
    return(1 - a * c)
  }
  return(0)
}

# Randkorrektur: Simuliere den Elternprozess auf [-1/c, 1 + 1/c]^2
N <- rpois(1, lambda0 * (2 + c)^2 / c^2)

x0 <- runif(N) * (1 + 2/c) - 1/c
y0 <- runif(N) * (1 + 2/c) - 1/c

x1 <- c()
y1 <- c()

for(i in 1:N){
  x2 <- c()
  y2 <- c()
  
  N1 <- rpois(1, 4/c^2)
  
  x3 <- runif(N1) * 2/c - 1/c
  y3 <- runif(N1) * 2/c - 1/c
  u <- runif(N1)
  
  if(N1 > 0){
    for(j in 1:N1){
      if(f(x3[j], y3[j], c) <= u[j]){
        x2 <- c(x2, x3[j]) 
        y2 <- c(y2, y3[j]) 
      }
    }
    
    
    x2 <- x2 + x0[i]
    y2 <- y2 + y0[j]
    
    tmpx <- x2[0 <= x2 & x2 <= 1 & 0 <= y2 & y2 <= 1]
    tmpy <- y2[0 <= x2 & x2 <= 1 & 0 <= y2 & y2 <= 1]
  
    x2 <- tmpx
    y2 <- tmpy
    
    x1 <- c(x1, x2)
    y1 <- c(y1, y2)
  }
  
}

print(x1)
print(y1)
plot(x1, y1, xlim=c(0,1), ylim=c(0,1))