calc_bulk_hyphil.cpp

#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
#include <cassert>
#include <cstdlib>
#include <vector>
#include <sstream>
 
#define INLINE __attribute__((always_inline))
 
//T should be double or float
template<typename T>
struct vector3{
  T x,y,z;
  
  //constructors
  INLINE vector3<T>(void){
    x = y = z = 0.0;
  }
  INLINE vector3<T>(const vector3<T> &rhs) : x(rhs.x), y(rhs.y), z(rhs.z) {}
  INLINE vector3<T>(const T &r){
    x = y = z = r;
  }
  INLINE vector3<T>(const T *p){
    x = p[0]; y = p[1]; z = p[2];
  }
  INLINE vector3<T>(const T &x_,const T &y_,const T &z_){
    x = x_; y = y_; z = z_;
  }
 
  //operators
  INLINE const vector3<T> operator+(const vector3<T>& obj) const {
    vector3<T> ret(x+obj.x,y+obj.y,z+obj.z);
    return ret;
  }
  INLINE const vector3<T> operator-(const vector3<T>& obj) const {
    vector3<T> ret(x-obj.x,y-obj.y,z-obj.z);
    return ret;
  }
  INLINE const vector3<T>& operator+=(const vector3<T>& obj){
    x += obj.x; y += obj.y; z += obj.z;
    return *this;
  }
  INLINE const vector3<T>& operator-=(const vector3<T>& obj){
    x -= obj.x; y -= obj.y; z -= obj.z;
    return *this;
  }
  INLINE const vector3<T> operator*(const T& cf) const {
    vector3<T> ret(x*cf,y*cf,z*cf);
    return ret;
  }
  INLINE const T operator*(const vector3<T>& obj) const {
    const T ret = x*obj.x + y*obj.y + z*obj.z;
    return ret;
  }
  INLINE friend const vector3<T> operator*(const vector3<T>& obj , const T& cf) {
    vector3<T> ret(obj.x*cf,obj.y*cf,obj.z*cf);
    return ret;
  }
  INLINE const vector3<T>& operator*=(const T& cf){
    x *= cf; y *= cf; z *= cf;
    return *this;
  }
  INLINE const vector3<T>& operator/=(const T& cf){
    x /= cf; y /= cf; z /= cf;
    return *this;
  }
 
  //inverse square root
  //for IEEE754
  INLINE float m_rsqrt(const float& r) const {
    const float rHalf = 0.5 * r;
    const int   tmp   = 0x5F3759DF - ( *(int*)&r >> 1 );
    float rRes        = *(float*)&tmp;
    rRes *= ( 1.5 - ( rHalf * rRes * rRes ) );
    rRes *= ( 1.5 - ( rHalf * rRes * rRes ) );
    return rRes;
  }
  //for IEEE754
  INLINE double m_rsqrt(const double& r) const {
    const double         rHalf = 0.5 * r;
    const long long int  tmp   = 0x5FE6EB50C7B537AAl - ( *(long long int*)&r >> 1);//initial guess
    double         rRes        = *(double*)&tmp;
    rRes *= ( 1.5 - ( rHalf * rRes * rRes ) );
    rRes *= ( 1.5 - ( rHalf * rRes * rRes ) );
    return rRes;
  }
  
  //ret norm
  INLINE T invnorm2(void) const {
    const T dr2 = (*this)*(*this);
    return m_rsqrt(dr2);
  }
  INLINE T norm2(void) const {
    const T dr2 = (*this)*(*this);
    return sqrt(dr2);
  }
  INLINE T fastnorm2(void) const {
    const T dr2 = (*this)*(*this);
    return m_rsqrt(dr2) * dr2;
  }
  INLINE T dist2(void) const {
    const T dr2 = (*this)*(*this);
    return dr2;
  }
 
  INLINE void clear(void) {
    x = y = z = (T)0.0;
  }
    
  void elemprint(void) const {
    std::cout << "(" << x << "," << y << "," << z << ")" << std::endl;
  }
};
 
#ifndef __CUDACC__
typedef vector3<double> double3;
typedef vector3<float > float3;
#else
typedef vector3<double> dv3;
typedef vector3<float > df3;
#endif
 
#undef INLINE

int main(int argc,char* argv[]){
  using namespace std;
  assert(argc == 4);
  const string c_dir = argv[1];
  const int b_time = atoi(argv[2]);
  const double half_plane = atof(argv[3]);
  
  string buf_str;
  buf_str = c_dir + "/particle_data.txt";
  ifstream fin1(buf_str.c_str());
  buf_str = c_dir + "/macro_data.txt";
  ifstream fin2(buf_str.c_str());
  
  int wN,hbN;
  int all_time,time_step;
  double buf;
  double L[3];
  fin2 >> wN >> hbN >> L[0] >> buf >> all_time >> time_step;
  L[1] = L[0];
  L[2] = L[0];
  const int syssize = wN + hbN;

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

  const string tar_file = c_dir + "/bulk_concent.txt";
  ofstream fout(tar_file.c_str());
  fout << "#up     down    all" << endl;
  
  for(int time=0; time < all_time; time += time_step){
    for(int i=0; i<syssize; i++){
      double3 buf; int bprp; bool bchm;
      fin1 >> buf.x >> buf.y >> buf.z >> bprp >> bchm;
      if(bprp == 1 && bchm == false){
	pos.push_back(buf); prp.push_back(bprp); chm.push_back(bchm);
      }
    }
    
    if(time >= b_time){
      const size_t hyphil_size = pos.size();
      size_t count = 0;
      for(size_t i=0; i<hyphil_size; i++){
	if(pos[i].z > half_plane) count++;
      }
      fout << count << " " << hyphil_size - count << " " << hyphil_size << endl;
    }
    
    pos.clear();
    prp.clear();
    chm.clear();

  }
  
  return 0;
}