calcgyration.cpp

#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <cmath>
#include "Eigen/Core"
#include "Eigen/Eigenvalues"
#include <cstdlib>

//global val
double grid_leng = 2.0;
int grid_numb[3];
int all_grid_numb;
double L[3];

struct double3{
  double x,y,z;
  double3 operator+(const double3& obj){
    double3 ret;
    ret.x = x + obj.x;
    ret.y = y + obj.y;
    ret.z = z + obj.z;
    return ret;
  }
  double3 operator-(const double3& obj){
    double3 ret;
    ret.x = x - obj.x;
    ret.y = y - obj.y;
    ret.z = z - obj.z;
    return ret;
  }
  double3& operator*=(const double& cf){
    x *= cf;
    y *= cf;
    z *= cf;
    return *this;
  }
  double3& operator+=(const double3& obj){
    x += obj.x;
    y += obj.y;
    z += obj.z;
    return *this;
  }
  double3& operator-=(const double3& obj){
    x -= obj.x;
    y -= obj.y;
    z -= obj.z;
    return *this;
  }
  double3& operator/=(const double& cf){
    x /= cf;
    y /= cf;
    z /= cf;
    return *this;
  }
};

int GenHash(const double3& r){
  int i = (int)(r.x / grid_leng);
  int j = (int)(r.y / grid_leng);
  int k = (int)(r.z / grid_leng);
 
  if(i == grid_numb[0]) i--;
  if(j == grid_numb[1]) j--;
  if(k == grid_numb[2]) k--;
  
  const int hash = i + grid_numb[0] * (j + k * grid_numb[1]);
  assert(hash < all_grid_numb);
  return hash;
}
 
void MinImage(double3& a){
  a.x -= L[0] * (int) (2. * a.x / L[0]);
  a.y -= L[1] * (int) (2. * a.y / L[1]);
  a.z -= L[2] * (int) (2. * a.z / L[2]);
}

int main(int argc,char* argv[]){
  using namespace Eigen;
  using namespace std;

  const char* cur_dir = argv[1];
  int wN,hbN,sys_size;
  double buf;
  int all_time,time_step;

  //load system val  
  stringstream ss;
  ss  << cur_dir << "/macro_data.txt";
  ifstream fin(ss.str().c_str());
  fin >> wN >> hbN >> L[0] >> buf >> all_time >> time_step;
  sys_size = wN + hbN;
  L[1] = L[0];
  L[2] = L[0];

  grid_numb[0] = static_cast<int>(L[0] / grid_leng);
  grid_numb[1] = static_cast<int>(L[1] / grid_leng);
  grid_numb[2] = static_cast<int>(L[2] / grid_leng);
  all_grid_numb = grid_numb[0] * grid_numb[1] * grid_numb[2];
 
  grid_leng = L[0] / grid_numb[0];
  
  //for output
  ss.str("");
  ss << cur_dir << "/gyrationt.txt";
  ofstream fout(ss.str().c_str());

  //for load particle data
  ss.str("");
  ss << cur_dir << "/particle_data.txt";
  ifstream finp(ss.str().c_str());

  vector<double3> pos;
  vector<int>     prp;
  vector<bool>    chm;

  for(int time=0; time < all_time; time+=time_step){
    
    for(int i=0; i<sys_size; i++){
      double3 buf; int bprp; bool bchm;
      finp >> buf.x >> buf.y >> buf.z >> bprp >> bchm;
      if(bchm == 1){
	pos.push_back(buf); prp.push_back(bprp); chm.push_back(bchm);
      }
    }
    
    //get base pos
    double3 base_pos;
    const int base_sample = 6;
    int baseN = 1;
    int bef_hash=-1;
    for(size_t i=0; i<pos.size(); i++){
      double3 temp = pos[i];
      const int temp_hash = GenHash(temp);
      if(temp_hash == bef_hash){
	base_pos += pos[i];
	baseN++;
	if(baseN == base_sample) break;
      }else{
	baseN=1;
	base_pos *= 0.;
	base_pos += pos[i];
	bef_hash = temp_hash;
      }
    }
    base_pos *= 1./base_sample;
 
    //calc center of mass
    double3 cm_pos;
    for(size_t i=0; i<pos.size(); i++){
      double3 temp_pos = pos[i] - base_pos;
      MinImage(temp_pos);
      cm_pos += temp_pos;
    }
    cm_pos /= pos.size();
    
    double gyr[3][3] = {{0.,0.,0.},
			{0.,0.,0.},
			{0.,0.,0.}};
  
    for(size_t i=0; i<pos.size(); i++){
      double3 temp_pos = pos[i] - base_pos;
      MinImage(temp_pos);
      temp_pos -= cm_pos;
      gyr[0][0] += temp_pos.x * temp_pos.x;
      gyr[1][1] += temp_pos.y * temp_pos.y;
      gyr[2][2] += temp_pos.z * temp_pos.z;
      gyr[0][1] += temp_pos.x * temp_pos.y;
      gyr[0][2] += temp_pos.x * temp_pos.z;
      gyr[1][2] += temp_pos.y * temp_pos.z;
    }
  
    gyr[0][0] *= 1./pos.size();
    gyr[0][1] *= 1./pos.size();
    gyr[0][2] *= 1./pos.size();
    gyr[1][1] *= 1./pos.size();
    gyr[1][2] *= 1./pos.size();
    gyr[2][2] *= 1./pos.size();

    gyr[1][0] = gyr[0][1];
    gyr[2][0] = gyr[0][2];
    gyr[2][1] = gyr[1][2];
    
    Matrix3d gyrmat;
    for(int i=0; i<3; i++){
      for(int j=0; j<3; j++){
	gyrmat(i,j) = gyr[i][j];
      }
    }

    SelfAdjointEigenSolver<Matrix3d> es(gyrmat);
    if(es.info() != Success) std::abort();
    Vector3d tempvec = es.eigenvalues();
    double egvals[3] = {tempvec(0),tempvec(1),tempvec(2)};
    std::sort(egvals,egvals + 3);
  
    const double gyr_rad  = egvals[0] * egvals[0] + egvals[1] * egvals[1] + egvals[2] * egvals[2];
    const double aspheric = 1.5 * egvals[2] * egvals[2] - 0.5 * gyr_rad;
    const double acylind  = egvals[1] * egvals[1] - egvals[0] * egvals[0];
    const double shapeparam = (aspheric * aspheric + 0.75 * acylind * acylind) / (gyr_rad * gyr_rad);
    fout << aspheric << " " << acylind << " " << shapeparam << std::endl;

    //clear vector
    pos.clear();
    prp.clear();
    chm.clear();
  }
  return 0;
}