// COVID-19 infestion simulator prototype Ver 0.0 8-9th Apr.2020 Y.Okamoto  ref. from the following site  

// https://rad-it21.com/%E3%82%B5%E3%82%A4%E3%82%A8%E3%83%B3%E3%82%B9/michikoshi-shugo_20200331/

int imax = 100;     // x grid max
int jmax = 100;     // y gridmax
int kmax = 100;     // time step max
int dmax = 100;     // initial infected persons

int m;             // x grid for initial seeds 
int n;             // y grid 

float D = 10;      // grid span
float R = 8;       // circle radius 

int [][] F = new int[imax+1][jmax+1];     // flag for infected 
int [][] S = new int[imax+1][jmax+1];     // flag for initial seed (not use this program)
int [][] T = new int[imax+1][jmax+1];     // flag for temporary infection
int [] NN = new int[kmax+1];            // array for infected
int [] NNN = new int[kmax+1];           // array for recovered
int [] IF = new int[kmax+1];            // array for NN - NNN

color rd = color(255, 0, 0);   // red
color bl = color(0, 0, 255);   // blue
color gr = color(0, 255, 0);   // green

float IP = 0.25;  // infection rate (alpha of reference)
float RP = 0.05;  // recovery rate (gamma of ref.)
float VP = 0.5;  // void rate (no human space, n_void of ref.)

String o_f_name = "IP"+IP+"_RP"+RP+"_VP"+VP ;    // outut file name (parameters)

PrintWriter output_ct;     

// init. set up
void setup() {
  size(1000,1000);  // window size
  randomSeed(0);     // set up for random seed 
  background(255);   // background color (white)
// array initialize  
  for (int j = 0; j < jmax; j++) {
    for (int i = 0; i < imax; i++) { 
      F[i][j] = 0;    // susceptible initialize  -> may be no need (definite in the next term)
      T[i][j] = 0;    // temporay array infected initialize
    }
  } 
  
// init.cond.of human and space  
  for (int j = 0; j < jmax; j++) {
    for (int i = 0; i < imax; i++) {
         if (VP < random(1)){      // recovery threshold
            F[i][j] = 0;           // recovery flag
            stroke(0);
            fill(255);             // susseptible person color white 
            ellipse(float(i)*D,float(j)*D, R, R);   
         }
         else {
            F[i][j] = -1;           // recovery flag
            stroke(0);
            fill(88);               // void color gray
            ellipse(float(i)*D,float(j)*D, R, R);             
         }
         
    }
  }
  
// initial infected seeds 
  int d = 0;  
  while( d < dmax){
    int m = int(random(imax));   //rundom set fire
    int n = int(random(jmax));
    if (F[m][n] == 0){
    F[m][n] = 1;    // seed flag
//    S[m][n] = 1;
    stroke(0);    // line color black
    fill(rd+bl);  // paint color purple           
    ellipse(float(m)*D,float(n)*D, R, R);   //circle plot for infected seeds               
    d = d+1;   // renewal d
    }
  }    
}

// contact process initialize 
int k = 0;   // counter for time step
int u = 0;   // counter for infected 
int v = 0;   // counter for recovered

// main routine

void draw(){
    
// contact process

  for (int j = 1; j < jmax-1; j++) {
    for (int i = 1; i < imax-1; i++) { 
      if (F[i-1][j] == 1 || F[i][j-1] == 1 || F[i+1][j] == 1 || F[i][j+1] == 1){      // for neighbor sites   
        if (F[i][j] == 0){    // no infected
          if (IP > random(1)){    // threshold for ingection
            T[i][j] = 1;        // infected temporary flag
            stroke(0);
            fill(bl);            
            ellipse(float(i)*D,float(j)*D, R, R);           
          }
        }
      }
    }
  }

// renewal for array F  

  for (int j = 0; j < jmax; j++) {
    for (int i = 0; i < imax; i++) {
      if (T[i][j] == 1){ 
        F[i][j] = T[i][j];     // infected flag at each time
        u = u + 1;    // counter for infected
        T[i][j] = 0;    // infection initialize
      }
    }
  }
  
  NN[k]=u;        // sumation of infected at each time 
  
// recovery process

  for (int j = 0; j < jmax; j++) {
    for (int i = 0; i < imax; i++) {
      if (F[i][j] == 1){
         if (RP > random(1)){      // recovery threshold
            F[i][j] = 2;           // recovery flag
            stroke(0);
            fill(gr);            
            ellipse(float(i)*D,float(j)*D, R, R);  
            v = v + 1;             // recovery counter 
         }
      }
    }
  }

  NNN[k]=v;       // sumation of recovered at each time 
  
  IF[k]= u-v;     // temporary infected = cumulative infected - recovered

  k=k+1;            // time step advancing

  // end of calc (output csv file)
  
  if (k>kmax){     
    output_ct = createWriter(o_f_name + "_tmct.csv");     // output csv file name   
    for(int n=0; n < kmax; n++){
      output_ct.print((n+1)+" ");                         // time  
      output_ct.print(((imax-2)*(jmax-2)-(NN[n]))+" ");   // susceptible
      output_ct.print(NNN[n]+" ");                        // recovered (cumulative)
      output_ct.println(dmax+IF[n]+" ");                  // infected (survival)  
    }
  output_ct.flush();   // buffer flush
  output_ct.close();   // output file close
  noLoop();             // end of drwaing
  }
}