# Probability of extinction until tMax:

se <- 0.15               # Per-step standard deviation related to environmental stochasticity
tMax <- 10000            # Observation time (T in paper)
nSample <- 10000         # Number of trials for computation of survival probability

# Data frame to record survival probability depending on initial population size:
means <- data.frame(nStart=round(10^seq(0,5,0.2)),pSurvive=NA)

# Simulate survival probability for each value of nStart
for(i in seq(nrow(means))){
  samples <- data.frame(tex=rep(NA,nSample),tlim=NA,n=means$nStart[i])
  w <- seq(nSample)        # indices of surviving species
  t <- 0                   # time
  while(length(w)>0){
    # Too large n break rpois, so we split the sample up:
    w1 <- w[samples$n[w]<=1e6]
    w2 <- w[samples$n[w]>1e6]
    # Simulate one time step:
    samples$n[w1] <- rpois(length(w1),exp(rnorm(length(w1),0,se))*samples$n[w1])
    samples$n[w2] <- exp(rnorm(length(w2),0,se))*samples$n[w2] + rnorm(length(w2),0,sqrt(samples$n[w2]))
    t <- t+1
    if(t>tMax) break;
    w <- which(samples$n>0) # update indices of surviving species
  }
  means$pSurvive[i]<- length(w)/nSample # compute survivial probability
  print(means[i,])
}

ve <- se^2  # ve from paper
vd <- 1     # vd from paper; it's a constant since we use a Poisson distributions

doPDF <- F # generate PDF? doPDF <- T

if(doPDF) pdf(file = "pExt.pdf", width = 5, height = 5)

par(mfrow=c(2,1),mar=c(5,4,0.5,0.2))

xrange <- range(log(1+means$nStart*(ve/vd)))
xvalues <- seq(xrange[1],xrange[2],length.out=100)
xvalues <- log(1+means$nStart*(ve/vd))

plot(log(means$nStart),means$pSurvive,ylab="Survival probability",xlab=expression(log(italic(N)[0]))) # plot results
lines(log(means$nStart),xvalues*sqrt(2/pi/tMax/se^2),lty=2) # plot theoretical prediction
legend("topleft", legend=" ", title="(a)", bty='n')

plot(log(1+means$nStart*(ve/vd)),means$pSurvive,ylab="Survival probability",xlab=expression(italic(x)[0] == log(1+italic(N)[0]~italic(v)[e]/italic(v)[d])),xlim=xrange) # plot results
lines(xrange,xrange*sqrt(2/pi/tMax/se^2),lty=2) # plot theoretical prediction
legend("topleft", legend=" ", title="(b)", bty='n')

if(doPDF) dev.off()