// Version 2 - Fermions


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


class fermidist
{
  int Emax;

  int *edist;
  double *edistavg;
  double *tdist;

  int Etot;
  int Ntot;
  double B0;
  double beta;
  double Ef;

  int plotopt;

  GPlot plotwin;

  public:

  fermidist(int sEmax);
  ~fermidist();

  int shuffle();
  int cool();
  void computeavgdist(int N);
  void plot();
  void toggleplot();
  void show();

  void calctots();
  double Etotfn(double betaarg);
  int findbeta();
  void calctdist();

  int iter(int niter); 
  int cooliter(int niter);
  int coolby(int deltaE);
};



fermidist::fermidist(int sEmax)
{
  int i, j;

  Emax = sEmax;
  edist = (int *)malloc(Emax*sizeof(int));
  edistavg = (double *)malloc(Emax*sizeof(double));
  tdist = (double *)malloc(Emax*sizeof(double));

  if (!edist || !edistavg || !tdist) exit(1);

  for (j = 0; j < Emax; j++) edist[j] = 0;
  for (i = 0; i < 1100; i++)
  {
    if (rand() % 10 > 7) edist[i] = 1;
    else edist[i] = 2;
  }

  calctots();
  plotwin = newGPlot();

  calctdist();

  computeavgdist(50);
  plotopt = 2;
  plot();
}

fermidist::~fermidist()
{
  deleteGPlot(plotwin);
  free(edist); free(edistavg); free(tdist);
}

int fermidist::shuffle()
{
  int i1, i2, i1p, i2p, Eexc;
  i1 = rand() % Emax;
  i2 = rand() % Emax;

  Eexc = 1 + rand()%100;

  i1p = i1 - Eexc;
  i2p = i2 + Eexc;

  if (i1p < 0 || i2p > Emax || 
      edist[i1] <= 0 || edist[i2] <= 0 || 
      edist[i1p] >= 2 || edist[i2p] >= 2 ||
      i1 == i2 || i1p == i2p) return 0;
  else
  {
    edist[i1]--; edist[i1p]++;
    edist[i2]--; edist[i2p]++;
    return 1;
  }
}

int fermidist::cool()
{
  int i1, i1p, Eexc;
  i1 = rand() % Emax;

  Eexc = 1 + rand()%100;

  i1p = i1 - Eexc;

  if (i1p < 0 ||  
      edist[i1] <= 0 || 
      edist[i1p] >= 2) return 0;
  else
  {
    edist[i1]--; edist[i1p]++;
    Etot -= Eexc;
    return 1;
  }
}



void fermidist::computeavgdist(int N)
{
  //convolves n

  int i, j;
  int j1, j2;
  int sum;
  for (i = 0; i < Emax; i++)
  {
    j1 = i-N > 0 ? i-N : 0;
    j2 = i+N < Emax ? i+N : Emax-1; 
    sum = 0;
    for (j = j1; j <= j2; j++)
      sum += edist[j];
    edistavg[i] = (double)sum / (j2 - j1 + 1.0);
  }

  return;
}

void fermidist::show()
{
  printf("   Etot = %d\n", Etot);
  printf("   Ntot = %d\n", Ntot);
  printf("   B0 = %.10f\n", B0);
  printf("   beta = %f\n", beta);
  printf("   Ef = %f\n", Ef);
}

void fermidist::plot()
{
  int i;
  double *x;  
  double *y1;
  double *y2;
 
  x  = (double *)malloc(Emax * sizeof(double));
  y1 = (double *)malloc(Emax * sizeof(double));
  y2 = (double *)malloc(Emax * sizeof(double));

  if (!x || !y1 || !y2) exit(1);

  for (i = 0; i < Emax; i++) 
  {
    x[i]  = (double)i;
    y1[i] = (double)edist[i];
    y2[i] = (edistavg[i]); 
  }


  useGPlot(plotwin);

  for (i = Emax-1; i > 0; i--)
    if (tdist[i] > 0.001 || edist[i] > 0) break;
  xrange(0,i+1);
  xlabel("E");
  ylabel("N(E)");
  yrange(0,2);
 
  plot_line(Emax * (plotopt!=1), x, y1);
  add_line(Emax * (plotopt>=1), x, y2);
  add_line(Emax, x, tdist);

  free(x); free(y1); free(y2);
}

void fermidist::toggleplot()
{
  plotopt = (plotopt + 1) % 3;
  plot();
  return ;
}


int fermidist::iter(int niter)
{
  int c = 0;

  if (niter <= 0 || niter > 1000000) 
  {
    printf("Bad iteration #\n");
    return 0;
  }

  while (c < niter) c += shuffle();
  computeavgdist(50);
  plot();
  return 1;
}

int fermidist::cooliter(int niter)
{
  int c = 0;
  int Elimit = (int)(Ntot*Ntot/4.0+0.5);

  if (niter <= 0 || niter > 1000000) 
  {
    printf("Bad iteration #\n");
    return 0;
  }

  while (c < niter && Etot > Elimit) c += cool();
  computeavgdist(50);
  calctdist();
  plot();   
  return 1;  
}

int fermidist::coolby(int deltaE)
{
  int targetE = Etot - deltaE;
  int Elimit = (int)(Ntot*Ntot/4.0+0.5);
  if (targetE <= 0 || targetE >= Etot) 
  {
    printf("Bad deltaE #\n");
    return 0;
  }
  if (targetE < Elimit) targetE = Elimit;

  while (Etot > targetE) cool();
  computeavgdist(50);
  calctdist();
  plot();   
  return 1;  
}



void fermidist::calctots()
{
  int i = 0;
  Etot = 0;
  Ntot = 0;
  
  for (i = 0; i < Emax; i++)
  {
    Etot += i * edist[i];
    Ntot += edist[i];
  }

  return;
}

double fermidist::Etotfn(double argbeta)
{
  //uses the argument beta, not the internal beta, for calc
  int i;
  double rEtot = 0.0;

  beta = argbeta;
  B0 = (1.0 - exp(beta*(Ntot/2.0-Emax))) / (exp(beta * Ntot / 2.0) - 1.0);
  for (i = 0; i < Emax; i++)
    tdist[i] = 2.0 / (1.0 + B0 * exp(beta * (double)i));
  
  for (i = 0; i < Emax; i++)
    rEtot += (double)i * tdist[i];

  return rEtot; 
}

int fermidist::findbeta()
{
  int i;
  double beta0 = Ntot / (Etot - Ntot*Ntot/8.0);
  double betalo = 0.3 * beta0;
  double betahi = 10.0 * beta0;
  double betamid;

  double vallo = Etotfn(betalo) - (double)Etot;
  double valhi = Etotfn(betahi) - (double)Etot;
  double valmid;


  if (vallo * valhi > 0.0) return 0;

  for (i = 0; i < 30; i++)
  {
    betamid = (betalo + betahi)/2.0;
    valmid = Etotfn(betamid) - (double)Etot;

    if ((vallo < 0.0) ? (valmid >= 0.0) : (valmid < 0.0)) //opposite signs
      betahi = betamid;  
    else                                         //same signs
      betalo = betamid;
  }

  beta = (betalo + betahi)/2.0;

  return 1;
}

void fermidist::calctdist()
{
  int i;
  
  findbeta();
  Etotfn(beta);

  Ef = -log(B0)/beta;

  return;
}


int main(int argc, char **argv)
{
  int i;
  int c;
  int niter;
  char command[20];  

  srand(10);

  fermidist ee(5000);
     
  printf("\nCommands:\n");
  printf("> i# to iterate for # iterations\n");
  printf("> c# to cool for # iterations\n");
  printf("> b# to cool by approx. # energy\n");
  printf("> r to print summary\n");
  printf("> p to print graph\n");
  printf("> s* to save graph\n");
  printf("> m to toggle smoothed/exact plots\n");
  printf("> x to exit\n\n");


  while (1)
  {
    printf("> ");
    scanf("%s",command);

    if (command[0] == 'i') ee.iter(atoi(command+1));
    if (command[0] == 'b') ee.coolby(atoi(command+1));
    if (command[0] == 'c') ee.cooliter(atoi(command+1));
    if (command[0] == 'r') ee.show();
    if (command[0] == 's') 
    {
      if (command[1] == 0) save_plot("plot.ps");
      else save_plot(command+1);
    }
    if (command[0] == 'p') print_plot();
    if (command[0] == 'm') ee.toggleplot();
    if (command[0] == 'x') break; 
  }
  
  return 0;
}