//This simulates the TASEP with periodic BC using random sequential update
//The generated output file contains profiles at equidistant times
//written (in bad style) by Stefan Grosskinsky


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <sys/types.h>
#include <time.h>

#define l 200
#define rhol 0.2
#define rhor 0.8

FILE *output_ptr;
char *output_file;
char dataname[80];

int e[l];     //# of cars, configuration array
//int c;           //current counter
double t;        //time counter


void init();
void update();
void file();
void output();

int main(void)
{
  int k, delta;
  int tmeas;
  time_t t1;
  (void) time(&t1);
  srand48((long) t1);  //random seed with system time t1
  printf("This is a very good traffic model for Q2.3\n OUTPUT FILE CONTAINS CONFIGURATIONS!\n\n");
  file();              //creat output file
  //  tequi=l;             //set equilibration time
  tmeas=l;             //set waiting time between outputs
  delta=200;            //set number of outputs
  init();            //initialze system
  output();          //output initial config
  for(k=1; k<delta; k++)
  {  
    t=0;
    //    while(t<tequi) update();   //equilibrate
    //    t=0; c=0;
    while(t<tmeas) update();   //update
    output(); // Output the configuration
  }
  return 0;
}


/*void init()    //uniformly distribute n particles on l sites
{
  int i, j, x[l];
  for(i=0; i<l; i++)
  {
    x[i]=i;
    e[i]=0;
  }
  for(i=0; i<n; i++)
  {
    j=lrand48()%(l-i);
    e[x[j]]++;
    x[j]=x[l-i-1];
  }
}
*/

void init()    //distribute particles indep. on left and right half
{
  int i, j, x[l];
  double r, rho1, rho2;
  rho1=rhol;
  rho2=rhor;
  for(i=0; i<(l/2); i++)
  {
    r=drand48();
    if( r<rho1 ) e[i]=1;
    else e[i]=0;
  }
  for(i=(l/2); i<l; i++)
  {
    r=drand48();
    if( r<rho2 ) e[i]=1;
    else e[i]=0;
  }
}


void update()
{
  int i, j;
  double r;
  i=lrand48()%(l+1);
  if( i==l )            //left boundary event
  {
    if( e[0]==0 )
    {
      r=drand48();
      if( r<rhol ) e[0]=1;
    }
  }
  else
  {
    if( i==(l-1) )    //right boundary event
    {
      if( e[l-1]==1 )
      {
	r=drand48();
	if( r>rhor ) e[l-1]=0;
      }
    }
    else             //bulk event
    {
      j=(i+1);
      if( e[i]*(1-e[j])==1 )
      {
	e[i]=0;      //move particle with rate 1
	e[j]=1;
	//	c++;         //update current counter
      }
    }
  }
  t+=1./(l+1.);
}


void file()    //create the output file
{
  sprintf(dataname,"traffic_config_l%i.dat",l);
  output_file = dataname;
  output_ptr = fopen(output_file,"w");
  fclose(output_ptr);
}

void output() // Output the current configuration
{
  output_ptr = fopen(output_file,"a"); 
  int i;
  for(i = 0; i < l; i++ ) {
    fprintf(output_ptr,"%i ",e[i]);
  }
  fprintf(output_ptr,"\n");
  fclose(output_ptr);

}