#include "RadGyr.h" 

double RadGyrFromStructure(file_records *f)
{
  int 
    an = f->atomnum ;
  NRVectP 
    X = dvector(1,an),
    Y = dvector(1,an),
    Z = dvector(1,an);
  int n = FileSelectToXYZ(f,X,Y,Z);
  AllocateNRVect3 v,u;
  double RadGyr ;
  v[1] = NRVect_Mean(X,n);
  v[2] = NRVect_Mean(Y,n);
  v[3] = NRVect_Mean(Z,n);

  u[1] = NRVect_SumSq(X,n) / n;
  u[2] = NRVect_SumSq(Y,n) / n;
  u[3] = NRVect_SumSq(Z,n) / n;
  
  RadGyr = sqrt(u[1] + u[2] + u[3] - (v[1]*v[1] + v[2]*v[2] + v[3]*v[3]));

  free_dvector(X,1,an);
  free_dvector(Y,1,an);
  free_dvector(Z,1,an);
  return RadGyr;
}

void DistScalDFromCentroidToSrfV (file_records *f)
{
  double RadGyr;
  AllocateNRVect3 v;
  int ii;
  
  CentroidFromStructure(f,v);
  RadGyr = RadGyrFromStructure(f);

  fprintf(stderr,"Centroid: %f %f %f\n",v[1],v[2],v[3]);
  fprintf(stderr,"RadGyr: %f\n",RadGyr);

  for (ii = 0;ii<f->atomnum; ii++) {
    atom_record * atom = f->atoms[ii];
    double dx = atom->x - v[1] ;
    double dy = atom->y - v[2] ;
    double dz = atom->z - v[3] ;
    atom->surface = sqrt(dx*dx + dy*dy + dz*dz);
    dx /= RadGyr ;
    dy /= RadGyr ;
    dz /= RadGyr ;
    atom->volume = sqrt(dx*dx + dy*dy + dz*dz);
  }
}