--- trunk/src/applications/staticProps/GofR.cpp 2005/02/09 17:08:22 306 +++ trunk/src/applications/staticProps/GofR.cpp 2005/02/13 20:36:24 328 @@ -42,19 +42,23 @@ #include #include #include "applications/staticProps/GofR.hpp" +#include "utils/simError.h" namespace oopse { -GofR::GofR(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2, double len) - : RadialDistrFunc(info, filename, sele1, sele2, len){ +GofR::GofR(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2) + : RadialDistrFunc(info, filename, sele1, sele2){ - histogram_.resize(nbins_); - avgGofr_.resize(nbins_); + deltaR_ = len_ /nRBins_; + + histogram_.resize(nRBins_); + avgGofr_.resize(nRBins_); + + setOutputName(getPrefix(filename) + ".gr"); } void GofR::preProcess() { - avgGofr_.resize(nbins_); std::fill(avgGofr_.begin(), avgGofr_.end(), 0.0); } @@ -68,12 +72,12 @@ void GofR::processHistogram() { double volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume(); double pairDensity = npairs_ /volume; - double pairConstant = ( 4.0 * PI * pairDensity ) / 3.0; + double pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0; for(int i = 0 ; i < histogram_.size(); ++i){ - double rLower = i * delta_; - double rUpper = rLower + delta_; + double rLower = i * deltaR_; + double rUpper = rLower + deltaR_; double volSlice = ( rUpper * rUpper * rUpper ) - ( rLower * rLower * rLower ); double nIdeal = volSlice * pairConstant; @@ -91,12 +95,15 @@ void GofR::collectHistogram(StuntDouble* sd1, StuntDou Vector3d pos1 = sd1->getPos(); Vector3d pos2 = sd2->getPos(); Vector3d r12 = pos1 - pos2; - - double distance = (pos1 - pos2).length(); + currentSnapshot_->wrapVector(r12); - int whichBin = distance / delta_; - histogram_[whichBin] ++; - npairs_++; + double distance = r12.length(); + + if (distance < len_) { + int whichBin = distance / deltaR_; + ++histogram_[whichBin]; + ++npairs_; + } } @@ -104,13 +111,20 @@ void GofR::writeRdf() { std::ofstream rdfStream(outputFilename_.c_str()); if (rdfStream.is_open()) { rdfStream << "#radial distribution function\n"; - rdfStream << "#selection1: " << selectionScript1_; - rdfStream << "#selection2: " << selectionScript2_; + rdfStream << "#selection1: (" << selectionScript1_ << ")\t"; + rdfStream << "selection2: (" << selectionScript2_ << ")\n"; + rdfStream << "#r\tcorrValue\n"; + for (int i = 0; i < avgGofr_.size(); ++i) { + double r = deltaR_ * (i + 0.5); + rdfStream << r << "\t" << avgGofr_[i]/nProcessed_ << "\n"; + } } else { } + + rdfStream.close(); } }