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

void PolymerStep(int N, double *PositionsX, double *PositionsY, double Size, double k,
                 double a, double *pR);
void PrintStep(int N, double *PositionsX, double *PositionsY);

int main(int argc, char **argv)
{
  int i;                /* loop variable           */
  int N;                /* number of monomers      */
  int Nrounds;          /* number of rounds        */
  double koverT;        /* k/T                     */
  double aoverT;        /* a/T                     */
  double Size;          /* size of simulation box  */
  double R;             /* total distance square   */
  double *PositionsX;   /* x positions of monomers */
  double *PositionsY;   /* y positions of monomers */
  FILE *gfp;            /* file pointer to gnuplot */

  /* interpret command line arguments */
  if (argc!=6) {
    fprintf(stderr,"polymer usage: polymer <N> <Nrounds> <k/T> <a/T> <box_size>\n");
    return(1);
  }
  N=atoi(argv[1]);
  Nrounds=atoi(argv[2]);
  koverT=atof(argv[3]);
  aoverT=atof(argv[4]);
  Size=atof(argv[5]);

  /* get memory for configuration */
  PositionsX=(double *)malloc(2*N*sizeof(double));
  if (!PositionsX) {
    fprintf(stderr,"Not enough memory !\n");
    return(2);
  }
  PositionsY=PositionsX+N;

  /* initialize positions */
  R=0.0;
  for(i=0;i<N;i++) {
    PositionsX[i]=Size/2*sin(2*3.1415*i/(double)N);
    PositionsY[i]=Size/2*cos(2*3.1415*i/(double)N);
    R+=PositionsX[i]*PositionsX[i]+PositionsY[i]*PositionsY[i];
  }

  /* get gnuplot file pointer */
  gfp=popen("gnuplot","w");
  if (!gfp) {
    fprintf(stderr,"Could not open gnuplot !\n");
    free((void *)PositionsX);
    return(3);
  }

  /* do simulation */
  for(i=0;i<Nrounds;i++) {
    PolymerStep(N,PositionsX,PositionsY,Size,koverT,aoverT,&R);
    if (!(i%1000)) {
      PrintStep(N,PositionsX,PositionsY);
      fprintf(gfp,"set title \"R=%8.2f %d Steps\"\n",sqrt(R),i);
      fprintf(gfp,"plot [-%d:%d] [-%d:%d] \"polymer.dat\" t \"\", \"polymer.dat\" t \"\" w l\n",((int)Size)+1,
              ((int)Size)+1,((int)Size)+1,((int)Size)+1);
      fflush(gfp);
    }
  }

  /* free memory */
  fclose(gfp);
  free((void *)PositionsX);

  return(0);
}

void PolymerStep(int N, double *PositionsX, double *PositionsY, double Size, double koverT,
                 double aoverT, double *pR)
{
  int n;          /* chosen monomer to move */
  double dX,dY;   /* change in position     */
  double dl;      /* change in length       */
  double dE;      /* change in energy       */

  /* choose a monomer */
  n=(int)(N*drand48());
  /* choose a displacement */
  dX=2*drand48()-1;
  dY=2*drand48()-1;
  /* don't do anything if outside box */
  if (fabs(PositionsX[n]+dX)>Size || fabs(PositionsY[n]+dY)>Size) return;
  /* calculate energy in original position and subtract energy in new position */
  if (fabs(PositionsX[n])<1.0 && fabs(PositionsY[n])<1.0)
    dE=aoverT;
  else
    dE=0.0;
  if (fabs(PositionsX[n]+dX)<1.0 && fabs(PositionsY[n]+dY)<1.0)
    dE-=aoverT;
  if (n>0) {
    dl=sqrt((PositionsX[n]-PositionsX[n-1])*(PositionsX[n]-PositionsX[n-1])+
            (PositionsY[n]-PositionsY[n-1])*(PositionsY[n]-PositionsY[n-1]))-1.0;
    dE-=koverT/2.0*dl*dl;
    dl=sqrt((PositionsX[n]+dX-PositionsX[n-1])*(PositionsX[n]+dX-PositionsX[n-1])+
            (PositionsY[n]+dY-PositionsY[n-1])*(PositionsY[n]+dY-PositionsY[n-1]))-1.0;
    dE+=koverT/2.0*dl*dl;
  }
  if (n<N-1) {
    dl=sqrt((PositionsX[n]-PositionsX[n+1])*(PositionsX[n]-PositionsX[n+1])+
            (PositionsY[n]-PositionsY[n+1])*(PositionsY[n]-PositionsY[n+1]))-1.0;
    dE-=koverT/2.0*dl*dl;
    dl=sqrt((PositionsX[n]+dX-PositionsX[n+1])*(PositionsX[n]+dX-PositionsX[n+1])+
            (PositionsY[n]+dY-PositionsY[n+1])*(PositionsY[n]+dY-PositionsY[n+1]))-1.0;
    dE+=koverT/2.0*dl*dl;
  }
  /* decide if we move */
  if (dE<0 || drand48()<exp(-dE)) {
    *pR-=PositionsX[n]*PositionsX[n]+PositionsY[n]*PositionsY[n];
    PositionsX[n]+=dX;
    PositionsY[n]+=dY;
    *pR+=PositionsX[n]*PositionsX[n]+PositionsY[n]*PositionsY[n];
  }
}

void PrintStep(int N, double *PositionsX, double *PositionsY)
{
  int i;     /* loop variable           */
  FILE *fp;  /* file pointer for output */

  fp=fopen("polymer.dat","w");
  if (fp) {
    for(i=0;i<N;i++)
      fprintf(fp,"%8.3f %8.3f\n",PositionsX[i],PositionsY[i]);
    fclose(fp);
  }
  for(i=0;i<10000000;i++);
}