/*
    ising.c

    Monte-Carlo-Simulation fuer das Ising-Modell.

    Grundprogramm in C von Fabian Braun (09.01.1998)

*/


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


/* Hier kommen zunaechst ein paar systemspezifische Angaben.
   Fuer ein Spin verwenden wir der Einfachheit halber ein Byte
   Speicherplatz                                                       */

/* Was ist ein Byte? */
#define  BYTE  char

/* Wie wird ein Spin dargestellt? */
#define  SPIN  BYTE 
#define  UP    +1
#define  DOWN  -1

#define  sqr(x)  ((x)*(x))


/* MonteCarlo-Simulationen benoetigen ZUFALLSZAHLEN.
   Haeufig sind vom Compiler oder System zur Verfuegung gestellte Zufallszahlen
   unbrauchbar, da sie Korrelationen aufweisen. Wir verwenden daher einen 
   einfachen brauchbaren generator. 

   Der Generator ist entnommen aus:
   Numerical Recipes in C - The Art of Scientiffic Computing
   Second Edition
   W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery
   Cambridge University Press 1992

   Dort finden sich auch noch komplexere Generatoren.
*/


long    idum;

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define MASK 123459876

double ran0(long *idum)
{
	long   k;
	double ans;

	*idum ^= MASK;
	k=(*idum)/IQ;
	*idum=IA*(*idum-k*IQ)-IR*k;
	if (*idum < 0) *idum += IM;
	ans=AM*(*idum);
	*idum ^= MASK;
	return ans;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef MASK


/* Hier geht's endlich los: */
/* Zunaechst definieren wir einige Programmroutinen: */

/* randomize erzeugt ein Feld mit zufaelliger Spinverteilung. 
   Dabei wird die Anzahl der Up-Spins festgestellt. */

void randomize(int Nspins, SPIN *spin)
{
  int j;

  for(j=0; j<Nspins; j++)
    if(ran0(&idum) < 0.5)   /* Fifty-Fifty */
      spin[j]=UP;
    else
      spin[j]=DOWN;
}


/* innerEnergy berechnet die innere Energie H fuer eine vorgegebene 
   Spinkonfiguration. Dies wird nur einmal zur Berechnung der inneren 
   Energie der zufaellig erzeugten Konfiguration benoetigt.           */

double innerEnergy(int Nspins, SPIN *spin, double K, double hz, int *NspinsUp)
{
  double U;
  int    j;

  U = 0; *NspinsUp = 0;

  for(j=0; j<Nspins-1; j++)      /* Spin-Spin-Wechselwirkung */
    U -= K*spin[j]*spin[j+1]; 
  U -= K*spin[Nspins-1]*spin[0]; /* periodische Randbedingungen */
 
  for(j=0; j<Nspins; j++)        /* aeusseres Magnetfeld */
    U -= hz*spin[j];

  for(j=0; j<Nspins; j++)        /* zaehle Anzahl der Up-Spins */
    if(spin[j]==UP) (*NspinsUp)++;

  return U;
}


/* MCsteps fuehrt Nsteps Monte-Carlo-Schritte aus. */

void MCsteps(int Nsteps, int Nspins, SPIN *spin, 
	     double K, double hz, double beta, 
	     double *U, int *NspinsUp)
     /* Nsteps:   Anzahl der Monte-Carlo Schritte
        Nspins:   Systemgroesse 
        spin:     Zeiger auf Array von SPINs
	K, hz:    Kopplungen (siehe Aufgabe)
	beta:     inverse Temperatur
        U:        innere Energie
        NspinsUp: Anzahl der up-spins im System; daraus erhalten wir mit
                  Nspins die Magnetisierung                            */

{
  int    selectedSpin, prevSpin, nextSpin;
  double DeltaU, prob;
  int    DeltaSpinsUp;

  int    j;

  for(j=0; j<Nsteps; j++)
    {

      /* Schritt 1:  Bestimme einen zufaelligen Spin. */
      selectedSpin = Nspins*ran0(&idum);
      
      /* Finde Nachbarspins (unter Beruecksichtigung periodischer 
                             Randbedingungen). */
      prevSpin   = ( selectedSpin > 0 ? selectedSpin-1 : Nspins-1 ); 
      nextSpin   = ( selectedSpin < Nspins-1 ? selectedSpin+1 : 0 ); 

      /* Flippe Spin   */
      spin[selectedSpin] *= -1; 

      
      /* Schritt 2: Berechne die Aenderungen in Energie und NspinsUp */
      DeltaU       = -2*K*(spin[prevSpin]*spin[selectedSpin]+
			   spin[selectedSpin]*spin[nextSpin])
	             -2*hz*spin[selectedSpin];
      DeltaSpinsUp = spin[selectedSpin];
  
      /* Schritt 3: Energiebilanz */
      if(DeltaU<0)
	{ /* Schritt akzepieren */
	  *U        += DeltaU;
	  *NspinsUp += DeltaSpinsUp;
	} else {
	  /* Schritt 4: Nur mit Wahrscheinlichkeit akzeptieren */
	  prob = exp(-beta*DeltaU);
	  if(ran0(&idum)<prob)
	    { /* OK, Akzeptieren */
	      *U        += DeltaU;
	      *NspinsUp += DeltaSpinsUp;
	    } else /* Verwerfen: Flippe zurueck */
	      spin[selectedSpin] *= -1;
	}
    }
}


/* simulation fuehrt eine gesamte Monte-Carlo-Simultion durch.  

   Vorgaben: inverse Temperatur, Kopplungen K und hz, sowie die Systemgroesse 
   und die Anzahl der gewuenschten Messungen. Es wird die Magnetisierung und
   die spezifische Waerme gemessen. */

void simulation(int Nspins, SPIN *spin, double beta, double K, double hz, 
		int Nmessungen, double *magnetization, double *cV)
     /* Nspins:        Systemgroesse
	spin:          Vektor aus SPINs
        beta:          inverse Temperatur
        K, hz:         Kopplungen
        Nmessungen:    Anzahl der Messungen, ueber die gemittelt wird
        magnetization: liefert Magnetisierung
        cV:            liefert spezifische Waerme (pro Spin) */
{
  double U, ExpNUp, UMess, U2Mess;
  int    NspinsUp, j, NUpMess;

  /* Innere Energie vom Spinfeld ausrechnen und Up-Spins zaehlen: */
  U = innerEnergy(Nspins, spin, K, hz, &NspinsUp); 

  /* Ein erster Lauf ohne Messung, um das System aus der Anfangskonfiguration 
     ins Gleichgewicht zu bringen                                             */
  MCsteps(100*Nspins, Nspins, spin, K, hz, beta, &U, &NspinsUp);

  /* Jetzt werden "Messungen" gemacht */
  NUpMess = NspinsUp;             /* erste Messung ...        */
  UMess   = U;                    /* innere Energie PRO Spin  */
  U2Mess  = sqr(U);               /* Quadr. inn. En. PRO Spin */
  for(j=1; j<Nmessungen; j++)     /* Nmessungen-1 Schleifendurchlaeufe */
    {
      MCsteps(5*Nspins, Nspins, spin, K, hz, beta, &U, &NspinsUp);
      NUpMess += NspinsUp ;       /* ... und naechste Messung */
      UMess   += U;
      U2Mess  += sqr(U);
    }

  ExpNUp = ((double) NUpMess)/(double) Nmessungen; 

  *magnetization =  (2*ExpNUp-Nspins)/(double)(Nspins) ;
  *cV            =  sqr(beta)*(U2Mess/Nmessungen - sqr(UMess/Nmessungen))/
                    Nspins;
}


/* Hauptprogramm: 
   Gibt zwei Zahlentabellen aus: 
   a) Magnetisierung als Funktion von h
   b) spezifische Waerme als Funktion von T
*/


/* Systemgroesse */
#define NSPINS   100

int main(void)
{
  double beta, K, hz, hMax, magn, cV;
  int    steps, j;
  SPIN   spin[NSPINS+1];

  idum   = 24121917;       /* beliebige Zahl zur Initialisierung des 
                              Zufallszahlengenerator */
  steps  = 100;            /* Messpunkte auf den Ausgabekurven */
  steps--;

  /* Erste Tabelle: Magnetisierung als Funktion von h */

  beta = 1;                /* waehle 1/beta als Einheit */
  K    = beta;             /* waehle ein K            */

  /* Lasse nun h zwischen 0 und hMax laufen */
  hMax = 4*exp(-2*K)*beta;
  for(hz=0; hz<=hMax; hz+= hMax/steps)
    {
      /* Spinvektor wird initialisiert: */
      randomize(NSPINS, spin);
      /* MC-Simulation ausfuehren       */
      simulation(NSPINS, spin, beta, K, hz, 1000, &magn, &cV);
      printf("%lf\t%lf\n", hz, magn);
    }
  printf("*************************************************************\n");
  
  /* Um cV auszurechnen, fangen wir bei hohen Temperaturen an und 
     verringern diese dann langsam. Die Natur als Vorbild nehmend senken
     wir die Temperatur exponentiell langsam, was etwa dem Abkuehlen eines 
     Metalls entspricht (sog. "simulated annealing").                      */

  K    = 1;               /* waehle K als Einheit */
  hz   = 0;               /* kein Magnetfeld      */
  
  /* Spinvektor wird nur EINMAL initialisiert: */
  randomize(NSPINS, spin);
    
  for(j=-2;j<=steps-2;j++)
    {
      /* Senke die Temperatur exponentiell. Die Konstanten im folgenden
         Ausdruck laden zum herumprobieren ein...                       */
      beta = 1.0/(2.0*exp(-4.0*j/(steps+1)));            /* beta = 1/kT */
      simulation(NSPINS, spin, beta, K, hz, 1000, &magn, &cV);
      /* Die Konfiguration der alten Simulation wird aufbewahrt und als 
         Start fuer die naechste Konfiguration genommen. */

      printf("%lf\t%lf\n", 1/beta, cV);
    }
}