/****************************************************************/
/*** Simulates d-dim percolation                              ***/
/*** A.K. Hartmann June 2004                                  ***/
/*** Rechnergestütztes Wissenschaftliches Rechnen             ***/
/*** University of Goettingen, Germany 2004                   ***/
/****************************************************************/

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

#include "stacks.h"

#define INDEX(t, r, nn) ((t)*(nn)+r)
#define NEXT(t, r) next[INDEX(t, r, num_n)]

/******************* mf_init_arr() **********************/
/** Allocates and initialises 'next[][]' For simple    **/
/** hypercubic lattice. Use macro                      **/
/** 'NEXT' to access. NEXT(t,r) gives neighbour of     **/
/** element 't'=1,2,... in direction 'r' with          **/
/** 'r'=0,1,..,num_n-1 correspond to                   **/
/**  +x,-x,+y,-y,...                                   **/
/** PARAMETERS: (*)= return-paramter                   **/
/**       dim: dimension of system                     **/
/**         l: size of system in x,y,..                **/
/**   (*) n_p: ptr to number of sites                  **/
/** RETURNS:                                           **/
/**     pointer to neighbour array                     **/
/**      next: gives neighbours next[0..N][0..num_n-1] **/
/********************************************************/
int *percol_init_next(int dim, int *l, int *n_p)
{
  int t, t2;                            /* loop counter */
  int N;                             /* number of sites */
  int *n;    /* n[d]= # of sites in subsystem of dim. d */  
  int *next;        /* contains neighbours of all sites */
  int num_n;                    /* number of neighbours */

  num_n = 2*dim;
  n = (int *) malloc((dim+1)*sizeof(int));
  n[0] = 1;            /* calculates sizes of subsytems */
  for(t=1; t<dim+1; t++)
    n[t] = l[t-1]*n[t-1];
  N = n[dim];
  *n_p =N;
  next = malloc((N+1)*(2*dim)*sizeof(int));
  for(t=1; t<N+1; t++)           /* loop over all sites */
  {
    /*NEXT(t, 2*dim) = N+1;
      NEXT(t, 2*dim+1) = 0;*/
    for(t2=0; t2<dim; t2++)
    {
      if((((t-1)%n[t2+1])/n[t2])==0) 
      {                                /* "left" border */
	NEXT(t, 2*t2)= t+n[t2];       /* right neigbour */
	NEXT(t, 2*t2+1) = t+(l[t2]-1)*n[t2];
      }
      else if((((t-1)%n[t2+1])/n[t2])==(l[t2]-1))
      {                               /* "right" border */
	NEXT(t, 2*t2) = t-(l[t2]-1)*n[t2];  
	NEXT(t, 2*t2+1) = t-n[t2];    /* left neighbour */
      }
      else     
      {                                /* standard case */
	NEXT(t, 2*t2) = t+n[t2];
	NEXT(t, 2*t2+1) = t-n[t2];
      }
    }
  }
  free(n);
  return(next);
}




/******************* percol_cluster() *******************/
/** Calculates the connected clusters of the lattice   **/
/** 'site'. Occupied sites have 'site[i]=1' (i=1..N).  **/
/** Neigbours of occupied sites form clusters.         **/
/** For each site, in 'cluster[i]' the id of the       **/
/** cluster (starting at 0) is stored                  **/ 
/** PARAMETERS: (*)= return-paramter                   **/
/**       num_n: number of neighbours                  **/
/**           N: number of sites                       **/
/**        next: gives neighbours (0..N)x(0..num_n-1)  **/
/**              0 not used here. Use NEXT() to access **/
/**        site: tells whether site is occupied        **/
/** (*) cluster: id of clusters sites are contained in **/
/** RETURNS:                                           **/
/**     number of clusters                             **/
/********************************************************/
int percol_cluster(int num_n, int N, int *next, 
		   short int *site, int *cluster)
{
  int num_clusters=0;
  int t, r;       /* loop counters over sites, directions */
  int current, neighbour;                        /* sites */
  lstack_t *members;      /* stack of members for cluster */

  for(t=1; t<=N; t++)    /* initialise all clusters empty */
    cluster[t] = -1;
  members = lstack_new(N, sizeof(int));

  for(t=1; t<=N; t++)              /* loop over all sites */
  {
    if((site[t] == 1)&&(cluster[t]==-1)) /* new cluster ? */
    {
      lstack_push(members, &t);          /* start cluster */
      cluster[t] = num_clusters;
      while(!lstack_is_empty(members))  /* extend cluster */
      {
	lstack_pop(members, &current);
	for(r=0; r<num_n; r++)   /* loop over neighbours */
	{
	  neighbour = NEXT(current, r);
	  if((site[neighbour]==1)&&(cluster[neighbour]==-1))
	  {          /* neighbour belongs to same cluster */
	    lstack_push(members, &neighbour);
	    cluster[neighbour] = num_clusters;
	  }
	}
      }
      num_clusters++;
    }
  }

  lstack_delete(members);
  return(num_clusters);
}

/******************* percol_print() *********************/
/** Prints 'site' line by line. All occupied sites     **/
/** (site[i]==1, i=1..N) are indicated by printing X.  **/
/** PARAMETERS: (*)= return-paramter                   **/
/**           N: number of sites                       **/
/**           l: size of system in x,y,..              **/
/**        site: tells whether site is occupied        **/
/** RETURNS:                                           **/
/**     nothing                                        **/
/********************************************************/
void percol_print(int N, int *l, short int *site)
{
  int t;                                /* loop counter */
  
  for(t=1; t<=N; t++)            /* loop over all sites */
  {
    if(site[t])
      printf("X");
    else
      printf(" ");
    if( ((t-1)%l[0]) == (l[0]-1) )        /* end of line ? */
      printf("\n");
  }
}

/******************* percol_largest() *******************/
/** Picks largest cluster from all clusters            **/
/** PARAMETERS: (*)= return-paramter                   **/
/**           N: number of sites                       **/
/**      num_cl: number of clusters                    **/
/**     cluster: id of clusters sites are contained in **/
/** RETURNS:                                           **/
/**     size of largest clusters                       **/
/********************************************************/
int percol_largest(int N, int num_cl, int *cluster)
{
  int *size;             /* size[i] = size of cluster i */
  int t, cl;       /* loop counter over sites, clusters */
  int max_size;              /* size of largest cluster */

  size = (int *) malloc(num_cl*sizeof(int));    /* init */
  for(cl=0; cl<num_cl; cl++)
    size[cl] = 0;

  for(t=1; t<=N; t++)                /* determine sizes */
    if((cluster[t] >= 0) && (cluster[t] < num_cl))
      size[cluster[t]]++;

  max_size = size[0];              /* determine maximum */
  for(cl=1; cl<num_cl; cl++)
    if(size[cl] > max_size)
      max_size = size[cl];

  free(size);
  return(max_size);
}




/******************* percol_analysis() ******************/
/** Performs evaluation with given data array          **/
/** PARAMETERS: (*)= return-paramter                   **/
/**     num_pts: number of points                      **/
/**       value: data                                  **/
/** RETURNS:                                           **/
/**     Binder parameter                               **/
/********************************************************/
double percol_analysis(int num_pts, double *value)
{
  double sum, sum2, sum4;
  int t;

  sum = 0; sum2=0; sum4=0;
  for(t=0; t<num_pts; t++)
  {
    sum += value[t];
    sum2 += value[t]*value[t];
    sum4 += value[t]*value[t]*value[t]*value[t];
  }

  /*return(sum/num_pts);*/
  return(0.5*(3.0-(sum4/num_pts)/(sum2*sum2/(num_pts*num_pts))));
}


int main(int argc, char *argv[])
{
  int dim;                    /* dimension of system */
  int num_n;    /* number of neighbours of each site */
  int *l;                         /* sizes of system */
  int *next;           /* stores neighbours of sites */
  int t, r;                         /* loop counters */
  int N;                          /* number of sites */
  double prob;           /* probability of occ. site */
  short int *site;     /* site is occupied or empty? */
  int *cluster;           /* cluster id of each site */
  int num_clusters;             /* numb. of clusters */
  double *value;                   /* stores results */
  double mean;                       /* final result */
  int run, num_runs; /* cur., max no. of indep. runs */

  dim = 2;            /* initialization + parameters */
  num_n = 2*dim;                   /* simple lattice */
  l = (int *) malloc(dim*sizeof(int));
  l[0] = atoi(argv[1]);
  l[1] = l[0];
  next = percol_init_next(dim, l, &N);

  /*for(t=1; t<=N; t++)
    for(r=0; r<2*dim; r++)
    printf("NEXT(%d,%d)=%d\n",t, r, NEXT(t,r));
  exit(1);*/
  site = (short int *) malloc((N+1)*sizeof(short int));
  cluster = (int *) malloc((N+1)*sizeof(int));
  sscanf(argv[2], "%lf", &prob); 
  num_runs = atoi(argv[3]);
  value = (double *) malloc(num_runs*sizeof(double));

  for(run=0; run<num_runs; run++)
  {
    for(t=1; t<=N; t++)
      if(drand48() < prob)
	site[t] = 1;
      else
	site[t] = 0;

    num_clusters = percol_cluster(num_n, N, next, site, cluster);
    value[run] = (double )percol_largest(N, num_clusters, cluster)/N;

    /*for(t=1; t<=N; t++)          
    {
      if(size[t] == 1)
	printf("%3d", cluster[t]);
      else
	printf("   ");
      if( ((t-1)%l[0]) == (l[0]-1) ) 
	printf("\n");
	}*/

  }
  mean = percol_analysis(num_runs, value);
  printf("%f %f\n", prob, mean);

  

  free(value);
  free(cluster);
  free(site);
  free(next);
  free(l);
  return(0);
}