rm(list=ls());

#Paramètres
du=0.03;  Du=0.02;  au=0.08;  bu=-0.08; cu=0.04;
dv=0.08;  Dv=0.5;   av=0.1;   bv=0;     cv=-0.15;
Fmax=0.2;
Gmax=0.5;

# Pas de temps et d'espace
Dx<-1 #Dy sera identique
Dt<-0.1
Niter=6000 #Nombre de pas de temps
dim=100; #dimension du système

#Initialisation des matrices avec des valeurs de concentration autour de 1
Mu<- matrix(data=rnorm(dim*dim,1,0.01),nrow=dim,ncol=dim)  
Mv<- matrix(data=rnorm(dim*dim,1,0.01),nrow=dim,ncol=dim)

#filled.contour(z = Mu, nlevels = 100, color = terrain.colors)
setwd("D:/BOULOT/Enseignement/UERB35 - Parcours Science Santé - Ecole INSERM/R/Turing")

ipng=0

#Boucle d'évolution temporelle.
for (i in seq(1, Niter+1))
{
  #"./" permet de sauver dans le repertoire courant. Pour connaitre le répertoire 
  #courant, il faut utiliser la commande "print(getwd())"
  #Pour changer le répertoire courant utiliser setwd("chemin")
  #Sauvegarde d'une image tous les 100 pas de calcul
  if(i%%100 == 1)
  {
    #si moins de 100 images pour avoir une numérotation de 01 à 99.
    if(ipng<10) {
      out=paste0("./out0",toString(ipng),".png")
    } else { 
      out=paste0("./out",toString(ipng),".png")
    }
    png(file = out, width = 800, height = 800)
    #filled.contour(z = cbind(log(Mu),matrix(data=rep(0,5*dim),nrow=dim,ncol=5),log(Mv)), nlevels = 100, color = terrain.colors)
    filled.contour(z = Mu, nlevels = 100, color = terrain.colors)
    dev.off()
    ipng=ipng+1
  }
  
  #Remise à zero des matrices pour les calculs intermédiaires
  Fuv<- matrix(data=rep(0,dim*dim),nrow=dim,ncol=dim)
  Guv<- matrix(data=rep(0,dim*dim),nrow=dim,ncol=dim)
  Diffu<-matrix(data=rep(0,dim*dim),nrow=dim,ncol=dim)
  Diffv<-matrix(data=rep(0,dim*dim),nrow=dim,ncol=dim)

  #Calcul diffusion
  Diffu=Du*(-4*Mu
        +cbind(Mu[,2:dim],Mu[,1])
        +cbind(Mu[,dim],Mu[,1:dim-1])
        +rbind(Mu[2:dim,],Mu[1,])
        +rbind(Mu[dim,],Mu[1:dim-1,])
        )/(Dx^2)
  Diffv=Dv*(-4*Mv
            +cbind(Mv[,2:dim],Mv[,1])
            +cbind(Mv[,dim],Mv[,1:dim-1])
            +rbind(Mv[2:dim,],Mv[1,])
            +rbind(Mv[dim,],Mv[1:dim-1,])
            )/(Dx^2)
  #Calcul reaction
  Fuv=au*Mu+bu*Mv+cu; #calcul des valeurs de Fuv
  Fuv=Fuv*(Fuv>=0); #condition Fuv>=0
  Fuv<-Fuv*(Fuv<=Fmax)+(Fuv>Fmax)*matrix(data=rep(Fmax,dim*dim),nrow=dim,ncol=dim); #condition Fuv<Fmax
  
  Guv=av*Mu+bv*Mv+cv; #calcul des valeurs de Guv
  Guv=Guv*(Guv>=0); #condition Guv>=0
  Guv<-Guv*(Guv<=Gmax)+(Guv>Gmax)*matrix(data=rep(Gmax,dim*dim),nrow=dim,ncol=dim); #condition Fuv<Fmax
  
  #Calcul reaction
  Mu<-Mu+Dt*(Fuv-du*Mu+Diffu)
  Mv<-Mv+Dt*(Guv-dv*Mv+Diffv)
}

filled.contour(z = Mu, nlevels = 100, color = terrain.colors)


#Pour convertir les images en une animation
#Il faut installer la librairie: install.packages("av")
library(av)
# Liste les fichiers PNG (assurez-vous qu'ils sont dans le bon ordre)
files <- list.files(pattern = "out.*\\.png", full.names = TRUE)
files <- sort(files)  # S'assure que les images sont dans le bon ordre

# Créer la vidéo avec 10 fps
av_encode_video(files, output = "output.mp4", framerate = 10)