/* (klassische) MC-Simulation eines 'Quantendots'
   Numerische Methoden der Vielteilchenphysik WS01/02
   Roland Winkler und Guenter Zwicknagel */

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

double upot (int np, double alpha, double gamma,
             double partx[], double party[]);
float ran1 (int *idum);

void dswap (double **x, double **y);

/*****************/
/* Hauptprogramm */
/*****************/

int main (int argc, char *argv[])
{
  int np, anzges, i, j, idum, schritt, akzep;
  double alpha, gamma, zeta, upot_alt, upot_neu, du, upot_sum, rskala, upot_sq;
  double temp = 1.0;
  FILE *daten_file, *pos_file;
  size_t nps;
  double *partx, *party, *xn, *yn, *gp1mc, *rhor;
  int startsamp = 20000, anzmess = 1000, gpanz = 2000;

  /* Das Programm erwartet vier Argumente:                        */
  /* Teilchenzahl, alpha, gamma und Zahl der Schritte             */
  /* Ein optionales fuenftes Argument ist die Temperatur, mit der */
  /* alpha und gamma skaliert werden (default: temp = 1.0)        */

  if (argc == 5 || argc == 6) {
    sscanf (argv [1], "%i",  &np);       /* Teilchenzahl      */
    sscanf (argv [2], "%lf", &alpha);    /* alpha             */
    sscanf (argv [3], "%lf", &gamma);    /* gamma             */
    sscanf (argv [4], "%i",  &anzges);   /* Zahl der Schritte */
    if (argc == 6) {
      sscanf (argv [5], "%lf",  &temp);  /* Temperatur        */
      alpha /= temp;
      gamma /= temp;
    }
  }
  else {
    printf ("%s: Wrong number of arguments\n", argv [0]); 
    exit (1);
  };

  /* Felder alloziieren (nutze nicht das Element [0]) */
  nps   = (np+1) * sizeof (double);
  partx = malloc (nps);
  party = malloc (nps);
  xn    = malloc (nps);
  yn    = malloc (nps);
  gp1mc = malloc ((gpanz+1) * sizeof (double));
  rhor  = malloc ((gpanz+1) * sizeof (double));

  /* Die charakteristische Laengenskala */
  rskala = 4.0 / sqrt (alpha) + 5.0 * pow (gamma/alpha, 1.0/3.0);

  /* Random-Generator initialisieren */
  idum = - fabs (rskala*np*temp + 123.) - 0.5;
  printf ("Initialisierung: %f\n", ran1 (&idum));

  /* Teilchenorte initialisieren */
  zeta = 0.3 * rskala; 
  for (i = 1; i <= np; i++) {
    partx [i] = zeta * (ran1 (&idum) - 0.5);
    party [i] = zeta * (ran1 (&idum) - 0.5);
  };
  for (i = 0; i <= gpanz; i++) {
    rhor  [i] = 0.0;
    gp1mc [i] = 0.0;
  };

  /* Berechnung des Potentials fuer die Anfangskonfiguration */
  upot_alt = upot (np, alpha, gamma, partx, party);
  printf ("Startpotential   %e\n", temp * upot_alt / np);

  /* Ausgabe initialisieren */
  daten_file = fopen ("datqdclass.dat", "w");
  fprintf (daten_file, "# Teilchenzahl:      %i\n", np);
  fprintf (daten_file, "# Zahl der Schritte: %i\n", anzges);
  fprintf (daten_file, "# alpha:             %f\n", alpha);
  fprintf (daten_file, "# gamma:             %f\n", gamma);

  pos_file = fopen ("posqdclass.dat", "w");

  upot_sum = 0.;
  upot_sq  = 0.;
  akzep    = 0;
  zeta     = 0.02 * rskala;

  /* Beginn des Schleifenprogramms */
  for (schritt = 1; schritt <= anzges; schritt++) {

    /* Erzeugen einer neuen Konfiguration */
    for (i = 1; i <= np; i++) {
      xn [i] = partx [i] + zeta * (ran1 (&idum) - 0.5);
      yn [i] = party [i] + zeta * (ran1 (&idum) - 0.5);
    }

    /* Berechnung der Aenderung der pot. Energie */
    upot_neu = upot (np, alpha, gamma, xn, yn);
    du = upot_neu - upot_alt;

    /* Befrage das Orakel */
    if (exp (-du) > ran1 (&idum)) {
      akzep ++;
      upot_alt  = upot_neu;
      dswap (&partx, &xn);
      dswap (&party, &yn);
    }

    if (schritt == startsamp) printf ("Sample Potential %e\n",
                                      temp * upot_alt / np);

    /* Observable bestimmen */
    if (schritt > startsamp) {
      /* Potentielle Energie */
      double fak = temp * upot_alt;
      upot_sum += fak;
      upot_sq  += fak * fak;

      /* Dichte rho */
      for (i = 1; i <= np; i++) {
        int kl = (int) (sqrt (partx [i] * partx [i] + party [i] * party [i])
                        * gpanz / (0.65 * rskala));
        kl = (kl < gpanz) ? kl : gpanz;
        rhor [kl] ++;
      }

      /* Paarkorrelationsfunktion */
      for (i = 1; i <= np-1; i++) {
        for (j = i+1; j <= np; j++) {
          double rx = partx [i] - partx [j];
          double ry = party [i] - party [j];
          int kl = (int) (gpanz * sqrt (rx*rx + ry*ry) / rskala);
          kl = (kl < gpanz) ? kl : gpanz;
          gp1mc [kl] ++;
        }
      }
    }

    if ((schritt % anzmess) == 0) {
      /* Datenausgabe */
      if (schritt > startsamp) {
        int nn = np * (schritt - startsamp);
        double u1 = upot_sum / nn;
        double u2 = upot_sq  / nn;
        fprintf (daten_file, "%i %13.6f  %13.6e  %13.6e  %13.6f\n",
                 schritt, u1, sqrt ((u2 - u1 * u1) / nn),
                 (double) akzep / anzmess, zeta);
        for (i = 1; i <= np; i++)
          fprintf (pos_file, "%13.6e  %13.6e\n", partx[i], party[i]);
      }

      /* Anpassung des zeta-Parameters */
      if ((double) akzep / anzmess > 0.5)
        zeta *= 1.05;
      else
        zeta *= 0.95;

      akzep = 0;
    }
  }

  printf ("Final Potential  %e\n", temp * upot_alt / np);

  fclose (daten_file);

  /* Ausgabe der Dichte */
  daten_file = fopen ("rhoqdclass.dat", "w");
  {
    double deltr = 0.65 * rskala / gpanz;
    double fak   = 1.0 / (deltr * (anzges-startsamp) * np);
    fprintf (daten_file, "# r,   rho(r)\n");
    for (i = 0; i <= gpanz; i++)
      fprintf (daten_file, "%13.6e  %13.6e\n",
               (i + 0.5) * deltr, fak * rhor [i]);
  }
  fclose (daten_file);

  /* Ausgabe der Paarkorrelationsfunktion */
  daten_file = fopen ("grqdclass.dat", "w");
  {
    double deltr = rskala / gpanz;
    double fak   = 1.0 / (deltr * (anzges-startsamp) * np * (np - 1));
    fprintf (daten_file, "# r,   g(r)\n");
    for (i = 0; i <= gpanz; i++)
      fprintf (daten_file, "%13.6e  %13.6e\n",
               (i + 0.5) * deltr, fak * gp1mc [i]);
  }
  fclose (daten_file);

  return 0;
}

/***********************/
/* Potentielle Energie */
/***********************/

double upot (int np, double alpha, double gamma,
             double partx[], double party[])
{
  double u = 0.;
  int i, j;

  /* harmonischer Anteil */
  for (i = 1; i <= np; i++)
    u += partx [i] * partx [i] + party [i] * party [i];
  u *= alpha;

  /* Coulomb-Anteil */
  for (i = 1; i <= np-1; i++) {
    for (j = i+1; j <= np; j++) {
      double rx = partx [i] - partx [j];
      double ry = party [i] - party [j];
      u += gamma / sqrt (rx*rx + ry*ry);
    }
  }
  return u;
}

/***************/
/* Swap arrays */
/***************/

void dswap (double **x, double **y)
{
  double *temp;
  temp = *x;
  *x   = *y;
  *y   = temp;
}

/**************************************************/
/* Zufallszahlengenerator aus "Numercial Recipes" */
/**************************************************/

#define M1 259200
#define IA1 7141
#define IC1 54773
#define RM1 (1.0/M1)
#define M2 134456
#define IA2 8121
#define IC2 28411
#define RM2 (1.0/M2)
#define M3 243000
#define IA3 4561
#define IC3 51349

float ran1(int *idum)
{
  static long ix1,ix2,ix3;
  static float r[98];
  float temp;
  static int iff=0;
  int j;

  if (*idum < 0 || iff == 0) {
    iff=1;
    ix1=(IC1-(*idum)) % M1;
    ix1=(IA1*ix1+IC1) % M1;
    ix2=ix1 % M2;
    ix1=(IA1*ix1+IC1) % M1;
    ix3=ix1 % M3;
    for (j=1;j<=97;j++) {
      ix1=(IA1*ix1+IC1) % M1;
      ix2=(IA2*ix2+IC2) % M2;
      r[j]=(ix1+ix2*RM2)*RM1;
    }
    *idum=1;
  }
  ix1=(IA1*ix1+IC1) % M1;
  ix2=(IA2*ix2+IC2) % M2;
  ix3=(IA3*ix3+IC3) % M3;
  j=1 + ((97*ix3)/M3);
  temp=r[j];
  r[j]=(ix1+ix2*RM2)*RM1;
  return temp;
}

#undef M1
#undef IA1
#undef IC1
#undef RM1
#undef M2
#undef IA2
#undef IC2
#undef RM2
#undef M3
#undef IA3
#undef IC3