/*** mc.c:  simple Ising simulation   ***/

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

/** macros **/
/* for allocating memory in a nice way: **/
#define TFNE_MALLOC_AH(POINTER,ANZAHL,TYPE,STRING) \
{ POINTER = (TYPE *) malloc((ANZAHL)*sizeof(TYPE)); \
  if(POINTER == NULL) \
  {  fprintf(stderr, "Achtung: Kein Speicher fuer %s \n", STRING); exit(1);}}   
  
/* for using this macros the variables dim,next,j have to be       */     
/* (locally) defined !                                             */ 
/* use the macro 'NEXT' to access the neighbours of a spin         */
/* NEXT(t,r) gives the neighbour of spin with id 't'               */
/* in direction 'r', with r=0: +x direction, r=1: -x direction     */
/* r=2: +y direction, r=3: -y direction                            */
/* The subroutine mf_init_short takes care of the initalization.   */
/* You can use it as a black box, i.e. you don't have to under-    */
/* stand that subroutine.                                          */
/* An example how NEXT is used is give in the subroutine 'energy'. */
/* You have to read and to understand that subroutine!             */
#define INDEX(t, r, d) ((t)*(2*(d)+2)+r)    
#define NEXT(t, r) next[INDEX(t, r, dim)]   
/******************* mf_init_arr() **********************/    
/** Allocates and initialises 'n[]', 'N', 'next[][]',  **/    
/** 't_feld[]'                                         **/    
/** PARAMETERS: (*)= return-paramter                   **/    
/** dim: dimension of system                           **/    
/**   l: size of system in x,y,..                      **/    
/**    period: with pbc ?                              **/    
/**   (*) n_p: p. tosize of 0,1,.. dimensional plane   **/    
/**   (*) N_p: p. to number of spins = inner nodes     **/    
/**(*) next_p: gives neighbours next[0..N][0..2*dim+1] **/    
/**     in +x,-x,+y,-y,...,N+1 (sink), 0(source)       **/    
/** RETURNS:                                           **/    
/**     0-> OK                                         **/    
/********************************************************/    
int mf_init_short(int dim, int *l, int period, int **n_p, int *N_p ,
int **next_p)
{ 
  int t, t2, N;     
  int *n, *next;    
  
  
  TFNE_MALLOC_AH(n, dim+1, int, "n-Feld");  
  *n_p = n;   
  n[0] = 1;   
  for(t=1; t<dim+1; t++)  
    n[t] = l[t-1]*n[t-1]; 
  N = n[dim]; 
  *N_p = N;   
  TFNE_MALLOC_AH(next, (N+1)*(2*dim+2), int, "next-Array");   
  *next_p = next;   
  for(t=1; t<N+1; t++)    
  {     
    NEXT(t, 2*dim) = N+1; 
    NEXT(t, 2*dim+1) = 0; 
    for(t2=0; t2<dim; t2++)     
    {   
      if(!period)   
      { 
	if((((t-1)%n[t2+1])/n[t2])==0)    
	{    /* "linker" Rand */    
	  NEXT(t, 2*t2) = t+n[t2];  
	  NEXT(t, 2*t2+1) = 0;
	}   
	else if((((t-1)%n[t2+1])/n[t2])==(l[t2]-1))   
	{   /* "rechter" Rand */    
	  NEXT(t, 2*t2) = 0;  
	  NEXT(t, 2*t2+1) = t-n[t2];
	}    
	else
	{   
	  NEXT(t, 2*t2) = t+n[t2];  
	  NEXT(t, 2*t2+1) = t-n[t2];
	}   
      } 
      else    
      { 
	if((((t-1)%n[t2+1])/n[t2])==0)    
	{    /* "linker" Rand */    
	  NEXT(t, 2*t2)= t+n[t2]; /* recht. Nachf. */ 
	  NEXT(t, 2*t2+1) = t+(l[t2]-1)*n[t2];  
	}   
	else if((((t-1)%n[t2+1])/n[t2])==(l[t2]-1))   
	{   /* "rechter" Rand */    
	  NEXT(t, 2*t2) = t-(l[t2]-1)*n[t2];    
	  NEXT(t, 2*t2+1) = t-n[t2]; /* linker Nachf. */    
	}
	else
	{   
	  NEXT(t, 2*t2) = t+n[t2];  
	  NEXT(t, 2*t2+1) = t-n[t2];
	}   
      } 
    }   
  }     
  return(0);  
}  

/************************ energy() **********************/    
/** Calculates energy of 'dim'-dimensional Ising system**/
/** 'spin' haveing 'N' spins                           **/    
/** PARAMETERS: (*)= return-paramter                   **/    
/**   dim: dimension of system                         **/    
/**     N: number of spins                             **/
/**  next: neighbour relation, use NEXT() macro!       **/
/**  spin: [1]..[N]: values of spins, [0] is unused!   **/
/** RETURNS:                                           **/    
/**     energy per spin                                **/    
/********************************************************/    
double energy(int dim, int N, int *next, short int *spin)
{
  double energy = 0.0;
  int t,r;

  for(t=1; t<=N; t++)                                /* loop over all spins */
    for(r=0; r<2*dim; r++)                      /* loop over all neighbours */
      energy += -spin[t] * spin[NEXT(t,r)];          /* energy contribution */

  return(energy/(2*N));                       /* each bond is counted twice */
}


int main(int argc, char **argv)
{
  int  *next;              /* gives neighbours of spin i. Use macro NEXT ! */
  int t, r;                                               /* loop counters */ 
  int dim;                                         /* dimension of systems */ 
  int L;                                           /* length of system box */ 
  int *l;          /* array with sizes of system-box for x,y,... dimension */
  int N;                                      /* number of spins in system */
  int *n;                                                         /* dummy */
  short int *spin;                                      /* values of spins */
  int do_random;                            /* initialise spins randomly ? */

  do_random = atoi(argv[1]);                            /* read parameters */

  L = 30;   
  dim = 2;                                                /* initalisation */
  TFNE_MALLOC_AH(l, dim, int, "l");     
  for(r=0; r<dim; r++)      
    l[r] = L;   
  mf_init_short(dim, l, 1, &(n), &N, &(next));   

  TFNE_MALLOC_AH(spin, N+1, short int, "spin");          /* allocate spins */
  if(!do_random)
    for(t=1;t<=N; t++)                          /* initialise spins all up */
      spin[t] = 1;
  else
  {
    for(t=1; t<=N; t++)                   /* initialise spins all randomly */
      if(drand48() < 0.5)
	spin[t] = 1;
      else
	spin[t] = -1;
  }

  printf("%f\n",  energy(dim, N, next, spin));  
  printf("\n#");
  for(t=1; t<=N; t++)                               /* print configuration */
  {
    if(spin[t]>0)
	printf("+");
    else
	printf("-");
    if(t%L == 0)
	printf("\n#");
  }
  printf("\n");
  free(n); free(l); free(next); free(spin);                /* free memory */ 
  
  return(0);
}