#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <mpi.h>

#define MAX  5000                        /* number of iterations  */
#define SIZE 6001                         /* size of grid array  */
#define PI 3.1415926535897932385E0
#define COLOR 16                         /* num of colors */

int grid[SIZE][SIZE];

int main() {
  int ipe, npe; 
  int hit[256], ii, jj;
  double angle;
  int i, j, k, x, y, dist, dir, step, trav;
  int gauss_ran();			/* gaussian random number */
  MPI_Status status;
  MPI_Init(NULL, NULL);
  MPI_Comm_rank(MPI_COMM_WORLD,&ipe);
  MPI_Comm_size(MPI_COMM_WORLD,&npe);

  for(i = 0;i < SIZE; i ++){
    for(j=0; j<SIZE; j++) grid[i][j] = 0;
  }
  grid[SIZE/2][SIZE/2] = 1;			/* one particle at the center */ 
  srand48((unsigned)time( NULL ));
  for(k=0; k<COLOR; k++){ 
    if(ipe == 0){
      printf("\nLT%d\n",k);     
    }
    for(i=0; i<MAX; i++){
      MPI_Barrier(MPI_COMM_WORLD);
      hit[ipe] = -1;
      angle = (2*PI*drand48());		/* random angle */
      x     = (double)(1+cos(angle))*SIZE/2;      
      y     = (double)(1+sin(angle))*SIZE/2;
      
      dist = gauss_ran(); 	/* random number gaussian dist. */
      if (dist<0){
        step=-1;	/* move forwards or backwards */
      } else {
        step=1;
      }
      trav=0;
      while((hit[ipe]==-1)&&(x<SIZE-2)&&(x>1)&&(y<SIZE-2)&&(y>1)&&(trav<abs(dist))){
        if(grid[x+1][y]+grid[x-1][y]+grid[x][y+1]+grid[x][y-1]>=1){
          hit[ipe] = x*SIZE + y;
          grid[x][y] = 1;        
        } else if(drand48() < 0.5){
          x+=step;      /* move horizontally */
        } else {
          y+=step;      /* move vertically */
        } 
        trav++;
      }  
      for(ii = 0;ii < npe; ii++){ 
        MPI_Bcast(hit + ii, 1, MPI_INT, ii, MPI_COMM_WORLD);
        if(hit[ii] != -1){
          MPI_Bcast(&(grid[hit[ii]/SIZE][hit[ii]%SIZE]), 1, MPI_INT, ii, MPI_COMM_WORLD);
        }
      }
      if(ipe == 0){
        for(ii = 0; ii < npe; ii++){
          for(jj = 0; jj < ii; jj ++){
            if(hit[ii] == hit[jj]){
              hit[ii] = -1;
              break;
            }
          }
          if(hit[ii] != -1){
            printf("%d\t%d\n",hit[ii]/SIZE,hit[ii]%SIZE);
          }
        }
      }
    }     
  }
  return 0;
}
/*--------------------------end of main program-----------------------*/

/* generates random numbers with gaussian distribution using the */
/* Box-Mueller method */
int gauss_ran()
{
   double fac, rr, r1, r2;
   static int old=0;		/* have to be static so information */
   static int mem;		/* survives between function calls */
   
   if (old==0)			/* no random number left from */
   {				/* previous function call */
      do
      {
         r1= 2.0*drand48()-1.0;      /* choose random point in */	 
         r2= 2.0*drand48()-1.0;	     /* the unit circle */
         rr= r1*r1+r2*r2; 			 
      }while ((rr>=1)||(rr==0));    
   
      fac=sqrt(-2*log(rr)/rr);
      mem=5000*r1*fac;		     /* save for next call */
      old=1;			     /* set flag */
      
      return((int)(5000*r2*fac));	
   }
   else				     /* return second number */
   {				     /* from last call */
      old=0;			     /* unset flag */
      return mem;   		     /* return number from last call */
   }
}