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

#define l 1000

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

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


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

int main(void)
{
  int k, delta;
  int tequi, tmeas;
  time_t t1;
  (void) time(&t1);
  srand48((long) t1);  //random seed with system time t1
  printf("This is a very good traffic model...\n");
  file();              //creat output file
  tequi=l;             //set equilibration time
  tmeas=l;             //set measurement time
  delta=20;            //set delta\rho
  for(k=1; k<delta; k++)
  {  
    n=k*l/delta;
    init();            //initialze system
    t=0;
    while(t<tequi) update();   //equilibrate
    t=0; c=0;
    while(t<tmeas) update();   //measure
    printf("%g\n",c/(l*t));    //output
    output_ptr = fopen(output_file,"a");
    fprintf(output_ptr,"%g\n",c/(l*t));
    fclose(output_ptr);
  }
  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 update()
{
  int i, j;
  i=lrand48()%l;
  j=(i+1)%l;
  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;
}


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