/* Algorithm written by Yeng-Long Chen in Winter, 1999 */


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


void Rg(double *R, double *R2, double *Ix, double *Iy, double *Iz, double *Ixy, double *Ixz, double *Iyz);
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 R, R2, Ix, Iy, Iz, Ixy, Ixz, Iyz;
  double sumR, sumRx, sumRy, sumRz, sumRxy, sumRxz, sumRyz;
  double sumR2;
  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);
       
       sumRx = 0.;
       sumRy = 0.;
       sumRz = 0.;
       sumRxy = 0.;
       sumRxz = 0.;
       sumRyz = 0.;

       for(n=1; n<=nconfig; n++)
	 {
	   fscanf(stream, "%s %d %d", dummy1, &CN, &chainl);

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

	   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 periodic boundaries */
	   
	   unwrap(box);

	   /* Calculates Rg for each configuration */

	   Rg(&R, &R2, &Ix, &Iy, &Iz, &Ixy, &Ixz, &Iyz);

	   sumR = sumR +R;
	   sumR2 = sumR2 + R2;

	   /* Calculating the Moments of Inertia Ixx, Iyy, Izz, Ixy, Ixz, Iyz */

	   sumRx = sumRx +Ix;
	   sumRy = sumRy +Iy;
	   sumRz = sumRz +Iz;
	   sumRxy = sumRxy + Ixy;
	   sumRxz = sumRxz + Ixz;
	   sumRyz = sumRyz + Iyz;

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

  fclose(stream);

  R = sumR / nconfig;
  R2 = sumR2 / nconfig;

  Ix = sumRx / nconfig;
  Iy = sumRy / nconfig;
  Iz = sumRz / nconfig;
  Ixy = sumRxy / nconfig;
  Ixz = sumRxz / nconfig;
  Iyz = sumRyz / nconfig;

  printf("%d %d %lf %lf %lf %lf %lf\n", CN, chainl, R, sqrt((R2 - R*R)/(nconfig-1)), Ix, Iy, Iz);
}


void Rg(double *R, double *R2, double *Ix, double *Iy, double *Iz, double *Ixy, double *Ixz, double *Iyz)
{
  double Xcm, Ycm, Zcm;
  double dx, dy, dz, Ri2;
  double dxx, dyy, dzz, dxy, dxz, dyz; 
  double sumR, sumRx, sumRy, sumRz, sumRxy, sumRxz, sumRyz;
  double sumR2;
  int i, j;

  sumR = 0.;
  sumR2 = 0.;

  sumRx = 0.;
  sumRy = 0.;
  sumRz = 0.;
  sumRxy = 0.;
  sumRxz = 0.;
  sumRyz = 0.;
  
  for(i=1; i<=CN; i++)
    {
      Ri2 = 0;

      dxx = 0.;
      dyy = 0.;
      dzz = 0.;
      dxy = 0.;
      dxz = 0.;
      dyz = 0.;
      
      Xcm = 0.;
      Ycm = 0.;
      Zcm = 0.;

      /* Finding the Center-o-Mass of a given polymer*/

      for(j=1; j<=chainl; j++)
	{
	  Xcm = X[i][j] + Xcm;
	  Ycm = Y[i][j] + Ycm;
	  Zcm = Z[i][j] + Zcm;
	}

      Xcm = Xcm / chainl;
      Ycm = Ycm / chainl;
      Zcm = Zcm / chainl;

      /* Calculate Rg of a given polymer */

      for(j=1; j<=chainl; j++)
	{
	  dx = (X[i][j] - Xcm)*(X[i][j] - Xcm);
	  dy = (Y[i][j] - Ycm)*(Y[i][j] - Ycm);
	  dz = (Z[i][j] - Zcm)*(Z[i][j] - Zcm);

	  Ri2 = Ri2 + dx + dy + dz;

	  dxx = dxx + (Y[i][j] - Ycm)*(Y[i][j] - Ycm)+(Z[i][j] - Zcm)*(Z[i][j] - Zcm);
	  dyy = dyy + (X[i][j] - Xcm)*(X[i][j] - Xcm)+(Z[i][j] - Zcm)*(Z[i][j] - Zcm);
	  dzz = dzz + (X[i][j] - Xcm)*(X[i][j] - Xcm)+(Y[i][j] - Ycm)*(Y[i][j] - Ycm);
	  dxy = dxy + (X[i][j] - Xcm)*(Y[i][j] - Ycm);
	  dxz = dxz + (X[i][j] - Xcm)*(Z[i][j] - Zcm);
	  dyz = dyz + (Y[i][j] - Ycm)*(Z[i][j] - Zcm);
	}
      
      sumR = sumR + sqrt(Ri2/chainl);
      sumR2 = sumR2 + Ri2/chainl;

      sumRx = sumRx + dxx;
      sumRy = sumRy + dyy;
      sumRz = sumRz + dzz;
      sumRxy = sumRxy + dxy;
      sumRxz = sumRxz + dxz;
      sumRyz = sumRyz + dyz;
    }

  *R = sumR / CN;
  *R2 = sumR2 / CN;

  *Ix = sumRx / CN;
  *Iy = sumRy / CN;
  *Iz = sumRz / CN;
  *Ixy = sumRxy / CN;
  *Ixz = sumRxz / CN;
  *Iyz = sumRyz / CN;
}


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, "%lf %lf %lf\n", 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(*dx < (-box/2))
    *dz = *dz + box;
}