#####Preparing the export file
write.table(rbind(c("identifiant","Simulation","undrying.rate","growthrate.max","rain.nb.max","rain.radius.max","pop.initiale","nb.refuges","nb.migration","rain.k","0","]0-1]","]1-10]","]10-25]","]25-50]","]50-100]","]100-1000]","]1000-10000]","]>10000","year","rain", "Moy.pop","Pop.total","Capacite.total","Nb.colonise","Statut")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")
#####Loop for replicates (here 2000 replicats)
for(r in 1:5000){
  
  rm(list=ls())
  gc()
  
  #####Variables
  capacite.initial.desert <- 10
  capacite.initial.refuges.multiplicateur <- 1000
  capacitemin.refuges <- 100
  undrying.rate <- 0.35
  growthrate.max  <- 5
  rain.nb.max <- 5
  rain.radius.max <- 1
  rain.k <- 1000
  pop.initiale <- 2
  pop.refuges.initial <- 20
  nb.refuges <- 5
  nb.migration <- 3
  pop.min.refuges <- 20
  year <- 0
  simulation <- sample(1:10000000,1)
  identifiant <- c("ninanew",simulation)
  mat.simpson.sans.refuges <- matrix(capacite.initial.desert, nrow = 100, ncol = 100)
  i <- c(-1,0,1,10,25,50,100,1000,10000,10000000000)
  Variables <- rbind(c(identifiant,undrying.rate,growthrate.max,rain.nb.max,rain.radius.max,pop.initiale,nb.refuges,nb.migration,rain.k))
  zero <- rbind(c(NA,NA,NA,NA))
  quatorze <- rbind(c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA))
  neuf <- rbind(c(NA,NA,NA,NA,NA,NA,NA,NA,NA))
  
  rain.dataframe <- NA
  
  #####Simulation
  #clearing matrices and preparing tested variables 
  Rainpattern <- sample(c(5,12,25,37,50),1)
  maigre <- 6
  if(Rainpattern==5) grasselow <- 10
  if(Rainpattern==5) grassehigh <- 10
  if(Rainpattern==12) grasselow <- 29
  if(Rainpattern==12) grassehigh <- 49
  if(Rainpattern==25) grasselow <- 85
  if(Rainpattern==25) grassehigh <- 105
  if(Rainpattern==37) grasselow <- 141
  if(Rainpattern==37) grassehigh <- 161
  if(Rainpattern==50) grasselow <- 198
  if(Rainpattern==50) grassehigh <- 218
  
  rain.nb.max <- Rainpattern
  nb.refuges <- sample(c(1,10,50,100,200,350),1)
  nb.migration <- sample(c(1,2,3,4,5,8),1)
  simulation <- simulation+1
  identifiant <- c("ninanew",simulation)
  Variables <- rbind(c(identifiant,undrying.rate,growthrate.max,rain.nb.max,rain.radius.max,pop.initiale,nb.refuges,nb.migration,rain.k))
  write.table(rbind(c(Variables,quatorze,"DEBUT")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")
  
  mat.N <- matrix(pop.initiale, nrow = 100, ncol = 100)
  
  
  
  ##### Random instalation of the refuges
  mat.refuges <- matrix(1, nrow = 100, ncol = 100)
  for(o in 1:nb.refuges) {
    y <- round(runif(1, 1, 100))
    x <- round(runif(1, 1, 100))
    mat.refuges[x,y] <- capacite.initial.refuges.multiplicateur
  }
  mat.N[which(mat.refuges==capacite.initial.refuges.multiplicateur)] <- pop.refuges.initial
  ##### initial K of desert cell
  mat.simpson.K <- mat.simpson.sans.refuges*mat.refuges
  
  ####BIBLE TIME#######
  for(r in 1:9) { 
    #####  El Nino
    for(r in 1:7) {
      ##### Rain     
      year <- year+1
      mat.pluies <- matrix(0, nrow = 100, ncol = 100)
      rain <- sample(c(sample(1:maigre,1)),1)
      for(r in 1:rain) {
        y <- round(runif(1, 3, 97))
        x <- round(runif(1, 3, 97))
        mat.pluies[(x-rain.radius.max):(x+rain.radius.max),(y-rain.radius.max):(y+rain.radius.max)] <- rain.k;
      }
      mat.simpson.K <- mat.simpson.K+mat.pluies;
      ##### Growth    
      mat.N <- mat.N+(mat.N*growthrate.max*(1-(mat.N/mat.simpson.K)))
      #mat.N <- mat.N+(mat.N*growthrate.max)
      ##### Local extension
      mat.N[which(mat.N<1*runif(1))] <- 0
      ##### Migration
      for(m in 1:nb.migration){
        mat.emigre <- mat.N-mat.simpson.K
        mat.emigre[which(mat.emigre<0)] <- 0
        mat.N <- mat.N-mat.emigre
        pourcent <- mat.N*0.1
        pourcent[which(mat.refuges==capacite.initial.refuges.multiplicateur)] <- ifelse(mat.N[which(mat.refuges==capacite.initial.refuges.multiplicateur)]<pop.min.refuges, 0, pourcent[which(mat.refuges==capacite.initial.refuges.multiplicateur)])
        mat.emigre <- (mat.emigre+pourcent)/4
        mat.N <- mat.N-pourcent
        mat.emigre <- rbind(mat.emigre[100,],mat.emigre,mat.emigre[1,])
        mat.emigre <- cbind(mat.emigre[,100],mat.emigre,mat.emigre[,1])
        
        mat.imigre <- matrix(NA, nrow = 100, ncol = 100)
        for(y in 1:100) 
        {
          for (x in 1:100){
            mat.imigre[x,y] <- sum(mat.emigre[(x+1),(y+1+1)],mat.emigre[(x+1+1),(y+1)],mat.emigre[(x+1),(y+1-1)],mat.emigre[(x+1-1),(y+1)])
            
          }}
        mat.N <- mat.N+mat.imigre
      }
      ##### Elimination over K
      mat.mort <- mat.N-mat.simpson.K
      mat.mort[which(mat.mort<0)] <- 0
      mat.N <- mat.N-mat.mort
      
      ##### Ressource depletion
      mat.simpson.K <- mat.simpson.K*undrying.rate
      mat.simpson.K[which(mat.refuges==capacite.initial.refuges.multiplicateur)] <- ifelse(mat.simpson.K[which(mat.refuges==capacite.initial.refuges.multiplicateur)]<capacitemin.refuges, capacitemin.refuges, mat.simpson.K[which(mat.refuges==capacite.initial.refuges.multiplicateur)])
      
      Moy.pop.total <- mean(mat.N)
      freq.N.sans <- rbind(table(cut(mat.N[which(mat.refuges<capacite.initial.refuges.multiplicateur)], i)))
      data.N.sans <- as.data.frame(freq.N.sans)
      Nb.colonise <- length(mat.N[which(mat.N[which(mat.refuges<capacite.initial.refuges.multiplicateur)]>1)])
      Capacite.total <- sum(mat.simpson.K[which(mat.refuges<capacite.initial.refuges.multiplicateur)])  
      mat.N.sans <- mat.N
      mat.N.sans[which(mat.refuges==capacite.initial.refuges.multiplicateur)] <- 0
      mat.N.sans[which(mat.N.sans<1)] <- 0
      Pop.total <- sum(mat.N.sans)
      if(Nb.colonise>2000) { write.table(rbind(c(Variables,data.N.sans,year,rain,Moy.pop.total,Pop.total,Capacite.total,Nb.colonise,"pic")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")}
      if(Nb.colonise>0 & Nb.colonise<2000) { write.table(rbind(c(Variables,data.N.sans,year,rain,Moy.pop.total,Pop.total,Capacite.total,Nb.colonise,"faible")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")}
      if(Nb.colonise<1) { write.table(rbind(c(Variables,data.N.sans,year,rain,Moy.pop.total,Pop.total,Capacite.total,Nb.colonise,"mort")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")}
      
}
    #####  El Nina
    
    for(r in 1:2) {
      ##### Rain     
      year <- year+1
      mat.pluies <- matrix(0, nrow = 100, ncol = 100)
      rain <- sample(c(sample(grasselow:grassehigh,1)),1)
      for(r in 1:rain) {
        y <- round(runif(1, 3, 97))
        x <- round(runif(1, 3, 97))
        mat.pluies[(x-rain.radius.max):(x+rain.radius.max),(y-rain.radius.max):(y+rain.radius.max)] <- rain.k;
      }
      mat.simpson.K <- mat.simpson.K+mat.pluies;
      ##### Growth    
      mat.N <- mat.N+(mat.N*growthrate.max*(1-(mat.N/mat.simpson.K)))
      #mat.N <- mat.N+(mat.N*growthrate.max)
      ##### Local extension
      mat.N[which(mat.N<1*runif(1))] <- 0
      ##### Migration
      for(m in 1:nb.migration){
        mat.emigre <- mat.N-mat.simpson.K
        mat.emigre[which(mat.emigre<0)] <- 0
        mat.N <- mat.N-mat.emigre
        pourcent <- mat.N*0.1
        pourcent[which(mat.refuges==capacite.initial.refuges.multiplicateur)] <- ifelse(mat.N[which(mat.refuges==capacite.initial.refuges.multiplicateur)]<pop.min.refuges, 0, pourcent[which(mat.refuges==capacite.initial.refuges.multiplicateur)])
        mat.emigre <- (mat.emigre+pourcent)/4
        mat.N <- mat.N-pourcent
        mat.emigre <- rbind(mat.emigre[100,],mat.emigre,mat.emigre[1,])
        mat.emigre <- cbind(mat.emigre[,100],mat.emigre,mat.emigre[,1])
        
        mat.imigre <- matrix(NA, nrow = 100, ncol = 100)
        for(y in 1:100) 
        {
          for (x in 1:100){
            mat.imigre[x,y] <- sum(mat.emigre[(x+1),(y+1+1)],mat.emigre[(x+1+1),(y+1)],mat.emigre[(x+1),(y+1-1)],mat.emigre[(x+1-1),(y+1)])
            
          }}
        mat.N <- mat.N+mat.imigre
      }
      ##### Elimination over K
      mat.mort <- mat.N-mat.simpson.K
      mat.mort[which(mat.mort<0)] <- 0
      mat.N <- mat.N-mat.mort
      
      ##### Ressource depletion
      mat.simpson.K <- mat.simpson.K*undrying.rate
      mat.simpson.K[which(mat.refuges==capacite.initial.refuges.multiplicateur)] <- ifelse(mat.simpson.K[which(mat.refuges==capacite.initial.refuges.multiplicateur)]<capacitemin.refuges, capacitemin.refuges, mat.simpson.K[which(mat.refuges==capacite.initial.refuges.multiplicateur)])
      
      ##### Year results    
      Moy.pop.total <- mean(mat.N)
      freq.N.sans <- rbind(table(cut(mat.N[which(mat.refuges<capacite.initial.refuges.multiplicateur)], i)))
      data.N.sans <- as.data.frame(freq.N.sans)
      Nb.colonise <- length(mat.N[which(mat.N[which(mat.refuges<capacite.initial.refuges.multiplicateur)]>1)])
      Capacite.total <- sum(mat.simpson.K[which(mat.refuges<capacite.initial.refuges.multiplicateur)])  
      mat.N.sans <- mat.N
      mat.N.sans[which(mat.refuges==capacite.initial.refuges.multiplicateur)] <- 0
      mat.N.sans[which(mat.N.sans<1)] <- 0
      Pop.total <- sum(mat.N.sans)
      if(Nb.colonise>2000) { write.table(rbind(c(Variables,data.N.sans,year,rain,Moy.pop.total,Pop.total,Capacite.total,Nb.colonise,"pic")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")}
      if(Nb.colonise>0 & Nb.colonise<2000) { write.table(rbind(c(Variables,data.N.sans,year,rain,Moy.pop.total,Pop.total,Capacite.total,Nb.colonise,"faible")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")}
      if(Nb.colonise<1) { write.table(rbind(c(Variables,data.N.sans,year,rain,Moy.pop.total,Pop.total,Capacite.total,Nb.colonise,"mort")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")}
      
}
  }
  ##### Data
  
  if(Nb.colonise<1){write.table(rbind(c(Variables,quatorze,"EXTINCTION")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")}
  if(Nb.colonise>0 & Nb.colonise<2000) {write.table(rbind(c(Variables,quatorze,"SURVIE")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")}
  if(Nb.colonise>2000) {write.table(rbind(c(Variables,quatorze,"COLONISATION")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")}
  if(Moy.pop.total==0){write.table(rbind(c(Variables,quatorze,"ANNIHILATION")), file="data.rain.nina.csv", row.names=F, col.names=F, append=T, sep=",")}
  }
}