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

#define n 100         //system size
#define runs 100      //number of averages
#define steps 200     //number of outputs in the file
#define tmax 1000     //total (real) simulation time
#define l 1.7         //lambda

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

int x[n], sum;           //array of individuals and sum
double t, pr, pi;  


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

int main(void)
{
  int k,i, j,imax;
  time_t t1;
  (void) time(&t1);
  srand48((long) t1);  //random seed with system time t1
  printf("This is a very good program for the contact process...\n");
  pr=1./(2.*l);     // recovery probability for l>=1
  pi=0.5;           // infection prob for l>=1
  imax=tmax*n*2*l/steps;
  file();
  output_ptr = fopen(output_file,"a");
  for(j=0; j<runs; j++)
  {
    init();
    for(k=0; k<steps; k++)
    {
      for(i=0; i<imax; i++) update();
      //      printf("%i %g\n",sum,t);
      fprintf(output_ptr,"%i %g\n",sum,t);
    }
  }
  fclose(output_ptr);    
  return 0;
}


void init()
{
  int i;
  for(i=0; i<n; i++) x[i]=1;
  t=0;
  sum=n;
}


void update()
{
  int i, j;
  double p,r;
  i=lrand48()%n;
  p=x[i]*pr + (1-x[i])*pi*(x[(i+1)%n]+x[(i-1+n)%n]);
  r=drand48();
  if(r<p)
    {
      x[i]=1-x[i];
      sum+=2*x[i]-1;
    }
  t+=1./(n*2*l);
}

void file()
{
  sprintf(dataname,"contact_n%i_l%g.dat",n,l*10);
  output_file = dataname;
  output_ptr = fopen(output_file,"w");
  fclose(output_ptr);
}