/****************************************************************/
/*** Simulates 1d SOS models                                  ***/
/*** A.K. Hartmann June 2004                                  ***/
/*** Rechnergestütztes Wissenschaftliches Rechnen             ***/
/*** University of Goettingen, Germany 2004                   ***/
/****************************************************************/

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

/** data structure for global system **/
typedef struct
{
  int          L;                      /* size of system */
  int num_sweeps;        /* number of Monte Carlo sweeps */
  double       T;              /* temperature of surface */
} system_t;


/****************** random_deposition() **************/
/** Deposes particles randomly on the 1d surface    **/
/** given by 'h'. One sweep is performed            **/
/** PARAMETERS: (*)= return-parameter               **/
/**       glob: global data                         **/
/**     (*)  h: height field                        **/
/** RETURNS:                                        **/
/**    nothing                                      **/
/*****************************************************/
void random_deposition(system_t *glob, int *h)
{
  int pos;          /* position on x-axis of surface */
  int step;

  for(step=0; step<glob->L; step++)
  {
    pos = (int) floor(glob->L*drand48());
    h[pos]++;
  }
}


/***************** ballistic_deposition() ************/
/** Deposes particles ballistically on the 1d surf. **/
/** given by 'h'. One sweep is performed            **/
/** Uses periodic boundary conditions (pbc)         **/
/** PARAMETERS: (*)= return-parameter               **/
/**       glob: global data                         **/
/**     (*)  h: height field                        **/
/** RETURNS:                                        **/
/**    nothing                                      **/
/*****************************************************/
void ballistic_deposition(system_t *glob, int *h)
{
  int pos;          /* position on x-axis of surface */
  int neighb;       /* position of neighbour         */
  int step;

  for(step=0; step<glob->L; step++)
  {
    pos = (int) floor(glob->L*drand48());
    h[pos]++;
    if(pos==0)            /* set left neighbour, pbc */
      neighb = glob->L - 1;
    else
      neighb = pos - 1;
    if(h[neighb] > h[pos]) /* stick part. to neighb. */
       h[pos] = h[neighb];
    if(pos==(glob->L-1)) /* set right neighbour, pbc */
      neighb = 0;
    else
      neighb = pos + 1;
    if(h[neighb] > h[pos]) /* stick part. to neighb. */
       h[pos] = h[neighb];
  }
}



/********************* diffusion() *******************/
/** Performs one Monte Carlo sweep for diffusion    **/
/** on 1d surface given by 'h'                      **/
/** Uses periodic boundary conditions (pbc)         **/
/** Activated diffusion is used, see M. Siegert and **/
/** M. Plischke, Phys. Rev. E 50, 917 (1994)        **/
/** PARAMETERS: (*)= return-parameter               **/
/**       glob: global data                         **/
/**     (*)  h: height field                        **/
/** RETURNS:                                        **/
/**    nothing                                      **/
/*****************************************************/
void diffusion(system_t *glob, int *h)
{
  int pos;          /* position on x-axis of surface */
  int neighb;       /* position of neighbour         */
  int step;
  int energy_old, energy_new;
  double prob;            /* Monte Carlo probability */

  for(step=0; step<glob->L; step++)     /* one sweep */
  {
    pos = (int) floor(glob->L*drand48());
    if(drand48() < 0.5)           /* choose randomly */
    {
      if(pos==0)              /* left neighbour, pbc */
        neighb = glob->L - 1;
      else
	neighb = pos - 1;
    }
    else
    {
      if(pos==(glob->L-1))   /* right neighbour, pbc */
      neighb = 0;
    else
      neighb = pos + 1;
    }
    energy_old = (h[pos]-h[neighb])*(h[pos]-h[neighb]);
    energy_new = (h[pos]-h[neighb]-2)*(h[pos]-h[neighb]-2);
    prob = 1.0/(1.0+exp((energy_new-energy_old)/glob->T));
    if(drand48() < prob)            /* accept move ? */
    {
      h[pos] -= 1;      /* particle jumps from 'pos' */
      h[neighb] += 1;                 /* to 'neighb' */
    }
  }
}


/****************** reset_surface() ******************/
/** sets all heights of 1d surface to 0             **/
/** PARAMETERS: (*)= return-parameter               **/
/**       glob: global data                         **/
/**      (*) h: height field                        **/
/** RETURNS:                                        **/
/**    nothing                                      **/
/*****************************************************/
void reset_surface(system_t *glob, int *h)
{
  int pos;          /* position on x-axis of surface */

  for(pos=0; pos<glob->L; pos++)
    h[pos]=0;
}


/****************** print_surface() ******************/
/** prints heights of 1d surface                    **/
/** PARAMETERS: (*)= return-parameter               **/
/**       glob: global data                         **/
/**          h: height field                        **/
/** RETURNS:                                        **/
/**    nothing                                      **/
/*****************************************************/
void print_surface(system_t *glob, int *h)
{
  int pos;          /* position on x-axis of surface */

  for(pos=0; pos<glob->L; pos++)
    printf("%d %d\n", pos, h[pos]);
}

/********************** roughness() ******************/
/** Calculates roughness of 1d surface              **/
/** PARAMETERS: (*)= return-parameter               **/
/**       glob: global data                         **/
/**          h: height field                        **/
/** RETURNS:                                        **/
/**    roughness                                    **/
/*****************************************************/
double roughness(system_t *glob, int *h)
{
  int pos;          /* position on x-axis of surface */
  double sum, sum2;       /* sum of heights, squared */

  sum = 0.0;
  sum2 = 0.0;
  for(pos=0; pos<glob->L; pos++)
  {
    sum += h[pos];
    sum2 += h[pos]*h[pos];
  }

  return(sqrt(sum2/glob->L - (sum*sum)/(glob->L*glob->L)));
}



int main(int argc, char *argv[])
{
  system_t glob;    /* contains global system data */
  int *h;                          /* height field */
  int num_runs;      /* number of independent runs */
  int run, sweep, sweep2;         /* loop counters */
  double *sum_rough, *sum_rough2; /* roughn. at  t */
  double rough;

  glob.T = 0.2;
  glob.L = atoi(argv[1]);       /* read parameters */
  glob.num_sweeps = atoi(argv[2]);
  num_runs = atoi(argv[3]);
  h = (int *) malloc(glob.L*sizeof(int));  /* init */
  sum_rough = (double *) malloc((glob.num_sweeps+1)*sizeof(double));
  sum_rough2 = (double *) malloc((glob.num_sweeps+1)*sizeof(double));
  for(sweep=0; sweep<=glob.num_sweeps; sweep++)
  {
    sum_rough[sweep] = 0;
    sum_rough2[sweep] = 0;
  }
  srand48(1000);

  for(run=0; run<num_runs; run++)
  {
    reset_surface(&glob, h);
    /** main loop: **/
    for(sweep=1; sweep<=glob.num_sweeps; sweep++)
    {
      random_deposition(&glob, h);
      diffusion(&glob, h);
      rough = roughness(&glob, h);
      sum_rough[sweep] += rough;
      sum_rough2[sweep] += rough*rough;
    }
    /*print_surface(&glob, h);*/
  }

  for(sweep=1; sweep<=glob.num_sweeps; sweep++)
    printf("%d %f %f\n", sweep, sum_rough[sweep]/num_runs,
	   sqrt((sum_rough2[sweep]/num_runs-pow(sum_rough[sweep]/num_runs, 2))/
		(num_runs-1)));

  free(h);
  free(sum_rough);
  return(0);
}