/* Calculate gr by choosing each monomer as origin and then average over sum of g_i(r), modified to read in multiple configurations. g(r) is defined as the par-pair correlation funcion between monomers belonging to different polymers*/

/* Algorithm written with the help of Allen & Tildesley, Dr. Duane Johnson, and Mr. Mark Dewing */

/* Yeng-Long Chen, Winter 1999 */

#define BUFFER 1
#define PI 3.14159
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>



void gr(double *GR, double box, double rho, double Dr);
void sk(double *SK, double *SK2, double Dk, double Kl);

void unwrap(double box);
double **matrix(int num, int size);
void free_matrix(double **M, int num);
double *dvector(int size);
void error_exit(char *error_message);
void adjust(double *dx, double *dy, double *dz, double box);

double **X, **Y, **Z;  

int CN, chainl;  /* CN = number of chains, chainl = #of monomers/chain */


void main()
{
  int nconfig, nskip;
  int i, j, n, C, C_n, dummy;
  int bin, rbin, kbin, Kl;
  double box, rho, dummy2; 
  double Dr, *GR, *sumGR;
  double Dk, *SK, *SK2, *sumSK, *sumSK2, errsq;
  char dummy1[10];

  FILE *stream;
  
  nconfig = 1;       

  /* Read in the configuration file that contains multiple configurations of polymer coordinates generated by cbmc */

  if((stream = fopen("CONF.crd", "r"))!=NULL)
     {
       fscanf(stream, "%s", dummy1);
       fscanf(stream, "%s %d", dummy1, &nconfig);
       fscanf(stream, "%s %lf %lf %lf", dummy1, &box, &dummy2, &dummy2);
       fscanf(stream, "%s %lf %lf %lf", dummy1, &dummy2, &dummy2, &dummy2);

       /* Define the grid size of g(r) and s(k) */

       Dr = 0.05;
       Kl = 20;
       Dk = 2*PI / 50;
       
       rbin = (int)(box/2/Dr) +1;
       kbin = (int)(Kl /Dk) +1;
       
       GR = dvector(rbin);
       sumGR = dvector(rbin);
       
       SK = dvector(kbin);
       SK2 = dvector(kbin);
       sumSK = dvector(kbin);
       sumSK2 = dvector(kbin);
       
       for(n=1; n<=nconfig; n++)
	 {
	   fscanf(stream, "%s %d %d", dummy1, &CN, &chainl);

	   rho = (CN*chainl) / (box*box*box);

	   printf("\nCN=%d chainl=%d box=%lf, rho=%lf, nconfig=%d", CN, chainl, box, rho, nconfig);

	   X = matrix(CN, chainl);
	   Y = matrix(CN, chainl);
	   Z = matrix(CN, chainl);	   

	   for(i=1; i<=CN; i++)
	     for(j=1; j<=chainl; j++)
	       fscanf(stream, "%lf %lf %lf", &X[i][j], &Y[i][j], &Z[i][j]);

	   /* unwrap the periodic boundary conditions */

	   unwrap(box);
	   
	   /* find g(r) and s(k) of a given configuration */
	   
	   gr(GR, box, rho, Dr);
	   sk(SK, SK2, Dk, Kl); 
	   
	   for(bin = 1; bin<= rbin; bin++)
	     sumGR[bin] = sumGR[bin] + GR[bin];	   
	   for(bin = 1; bin<= kbin; bin++)
	     {
	       sumSK[bin] = sumSK[bin] + SK[bin];
	       sumSK2[bin] = sumSK2[bin] +SK2[bin];
	     }

	   free_matrix(X, CN);
	   free_matrix(Y, CN);
	   free_matrix(Z, CN);
	 }
     }
  else
    error_exit("Configuration File Not Found!");

  fclose(stream);
  printf("\n");

  /* Average over g(r) and s(k) of all configurations */

  for(bin = 1; bin<=rbin; bin++)
    GR[bin] = sumGR[bin] / nconfig;

  for(bin = 1; bin<=kbin; bin++)
    {
      SK[bin] = sumSK[bin] / nconfig;
      SK2[bin] = sumSK2[bin] / nconfig;
    }

  stream = fopen("gr.dat", "w");
  for(bin=1; bin<rbin; bin++)
    fprintf(stream, "%lf %lf\n", bin*Dr, GR[bin]);
  fclose(stream);

  stream = fopen("sk.dat", "w");
  for(bin=4; bin<=kbin; bin++)
    {
      errsq = (SK[bin]*SK[bin]-SK2[bin])/(CN*chainl*nconfig);
      fprintf(stream, "%lf %lf %lf\n", bin*Dk, SK[bin], errsq);
    }
  fclose(stream);    	      

  free(SK);
  free(SK2);
  free(GR);
  free(sumGR);
  free(sumSK);
  free(sumSK2);
}


void gr(double *GR, double box, double rho, double Dr)
{
  int i, j, m, n, Iseed;
  int bin, maxbin;
  int C, C_n;

  double hbox;
  double Nideal, Const;
  
  double dx, dy, dz, Rij; 
  double *Hist;
  double Rlower, Rupper; /* defining the differential space r to r+dr */


  hbox = box/2;
  
  /* The distance correlation is only good up to half the box length because of the peiordic boundaris */

  maxbin  = (int)(hbox/Dr) + 1;
  
  printf("\nmax R_bin = %d\n", maxbin);
  
  /* assign bins to histogram */
  
  Hist = dvector(maxbin);

  for(bin=1; bin<=maxbin; bin++)
    GR[bin] = 0.;
  
  /* Calculate g(r) centering on a random monomer */
  
  /* Build Histogram */
  for(C = 1; C <= (CN-1); C++)
    for(C_n = 1; C_n <= chainl; C_n++)
      for(i=(C+1); i<=CN; i++)
	{
	  if(i==C)
	    {
	      for(j=(C_n+1); j <= chainl; j++)
		{
		  dx = X[i][j] - X[C][C_n];
		  dy = Y[i][j] - Y[C][C_n];
		  dz = Z[i][j] - Z[C][C_n];
		  
		  adjust(&dx, &dy, &dz, box);
		  
		  Rij = sqrt(dx*dx + dy*dy + dz*dz);
		  bin = (int)(Rij / Dr) + 1;
		  
		  if(bin <= maxbin)
		    Hist[bin] = Hist[bin] + 2.;
		}
	    }
	  else
	    {
	      for(j=1; j<= chainl; j++)
		{
		  dx = X[i][j] - X[C][C_n];
		  dy = Y[i][j] - Y[C][C_n];
		  dz = Z[i][j] - Z[C][C_n];
		  
		  adjust(&dx, &dy, &dz, box);
		  
		  Rij = sqrt(dx*dx + dy*dy + dz*dz);
		  bin = (int)(Rij / Dr) + 1;
		
		  if(bin <= maxbin)
		    Hist[bin] = Hist[bin] + 2.;		
		}
	    }
	}

  printf("npart = %lf\n", npart);

  /* Normalize Histogram by local average */

  Const = 4.0 * PI * rho /3.;

  for(bin = 1; bin<=maxbin; bin++)
    {
      Rlower = (bin-1) * Dr;
      Rupper = Rlower + Dr;
  
      Nideal = Const * ((Rupper*Rupper*Rupper) - (Rlower*Rlower*Rlower));

      GR[bin] = Hist[bin] / (CN*chainl)/ Nideal;
    }

  free(Hist);
}
  


void sk(double *SK, double *SK2, double Dk, double Kl)
{
  int n, i, j, k, maxbin;
  int C, C_n;
  
  double *kx, *ky, *kz;
  double K, KR, csum, ssum;
  double *count;
  
  maxbin = (int)(Kl / Dk) + 1;

  printf("\nmax K_bin = %d", maxbin);
  
  kx = dvector(maxbin);
  ky = dvector(maxbin);
  kz = dvector(maxbin);
  
  count = dvector(maxbin);
  
  for(i=1; i<=maxbin; i++)
    {
      kx[i] = Dk*i; 
      ky[i] = Dk*i; 
      kz[i] = Dk*i; 
      SK[i] = 0.;
      SK2[i] = 0.;

    }

  /* since k is a vector, s(k) is calculated from the average of all k's of the same magnitudes, count[] keeps track of how many times each |k| is summed for the average*/

  for(i=1; i<=maxbin; i+=19)
    for(j=1; j<=maxbin; j+=19)
      for(k=1; k<=maxbin; k+=19)
	{
	  csum = 0.;
	  ssum = 0.;

	  K = sqrt(kx[i]*kx[i] + ky[j]*ky[j] + kz[k]*kz[k]);
	  
	  if ( K <= Kl )
	    {
	      n = (int) (K / Dk) +1;
	      count[n] = count[n] + 1.;
	      
	      for(C = 1; C <= CN; C++)
		for(C_n = 1; C_n <= chainl; C_n++)
		  {
		    KR = kx[n] * X[C][C_n] + ky[n]*Y[C][C_n] + kz[n]*Z[C][C_n];
		    
		    csum = csum + cos(KR);
		    ssum = ssum + sin(KR);
		  }
	    }	
	  else n=0;

	  SK[n] = csum*csum + ssum *ssum + SK[n];
	  SK2[n] = (csum*csum + ssum * ssum)*(csum*csum + ssum*ssum) +SK2[n];
	}
  
  for(n=1; n<=maxbin; n++)
    {
      /* normalize by nconfig and number of particles */

      if(count[n] !=0)
	{
	  SK[n] = SK[n] / (count[n]*CN * chainl);  
	  SK2[n] = SK2[n] / (count[n]*CN *chainl); 
	}
    }
  free(kx);
  free(ky);
  free(kz);
}


void unwrap(double box)
{
  int i,j;
  double dx, dy, dz;
  double hbox;
  FILE *stream;

  hbox = box /2;

  for(i=1; i<=CN; i++)
    for(j=1; j<chainl; j++)
      {
	dx = X[i][j+1] - X[i][j];
	dy = Y[i][j+1] - Y[i][j];
	dz = Z[i][j+1] - Z[i][j];
 
	if(dx > hbox)
	  X[i][j+1] = X[i][j+1] - box;
	if(dx < (-1*hbox))
	  X[i][j+1] = X[i][j+1] + box;

	if(dy > hbox)
	  Y[i][j+1] = Y[i][j+1] - box;
	if(dy < (-1*hbox))
	  Y[i][j+1] = Y[i][j+1] + box;

	if(dz > hbox)
	  Z[i][j+1] = Z[i][j+1] - box;
	if(dz < (-1* hbox))
	  Z[i][j+1] = Z[i][j+1] + box;
      }

  stream = fopen("newcord.dat", "w");
  for(i=1; i<=CN; i++)
    for(j=1; j<chainl; j++)
      fprintf(stream, "\n(%d %d) %lf %lf %lf", i,j, X[i][j], Y[i][j], Z[i][j]);
  fclose(stream);
}


double *dvector(int size)
{
  double *v;
  int i;
  if((v=(double *)malloc((size+BUFFER)*sizeof(double)))==NULL)
    error_exit("Vector Allocation Failed");
  for(i=0; i<=size; i++)
    v[i]=0.0;
  return v;
}

void error_exit(char *error_message)
{
  printf("\n%s\nProgram exiting...\n", error_message);
  exit(1);
}

double **matrix(int num, int size)
{
  int i;
  double **M;
  if((M=(double **)malloc((BUFFER+num)*sizeof(double *)))==NULL)
    error_exit("Matrix Allocation Failed");
  for(i=0; i<=num; i++)
    if((M[i]=(double *)malloc((size+BUFFER)*sizeof(double)))==NULL)
      error_exit("Matrix Allocation Failed");
  return(M);
}


void free_matrix(double **M, int num)
{
  int i;
  for(i=0; i<num; i++)
    free(M[i]);
  free(M);
}

void adjust(double *dx, double *dy, double *dz, double box)
{
  if(*dx > (box/2))
    *dx = *dx - box;
  else if(*dx < (-box/2))
    *dx = *dx + box;
  if(*dy > (box/2))
    *dy = *dy - box;
  else if(*dy < (-box/2))
    *dy = *dy + box;
  if(*dz > (box/2))
    *dz = *dz - box;
  else if(*dz < (-box/2))
    *dz = *dz + box;
}