--- trunk/src/applications/staticProps/GofRAngle.cpp 2005/02/10 18:14:03 310 +++ trunk/src/applications/staticProps/GofRAngle.cpp 2005/02/16 21:00:39 356 @@ -46,9 +46,19 @@ namespace oopse { namespace oopse { -GofRAngle::GofRAngle(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2) - : RadialDistrFunc(info, filename, sele1, sele2){ +GofRAngle::GofRAngle(SimInfo* info, const std::string& filename, const std::string& sele1, + const std::string& sele2, double len, int nrbins, int nangleBins) + : RadialDistrFunc(info, filename, sele1, sele2), len_(len), nRBins_(nrbins), nAngleBins_(nangleBins){ + deltaR_ = len_ /nRBins_; + deltaCosAngle_ = 2.0 / nAngleBins_; + + histogram_.resize(nRBins_); + avgGofr_.resize(nRBins_); + for (int i = 0 ; i < nRBins_; ++i) { + histogram_[i].resize(nAngleBins_); + avgGofr_[i].resize(nAngleBins_); + } } @@ -68,9 +78,10 @@ void GofRAngle::processHistogram() { void GofRAngle::processHistogram() { + int nPairs = getNPairs(); double volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume(); - double pairDensity = npairs_ /volume; - double pairConstant = ( 4.0 * PI * pairDensity ) / 3.0; + double pairDensity = nPairs /volume; + double pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0; for(int i = 0 ; i < histogram_.size(); ++i){ @@ -100,13 +111,14 @@ void GofRAngle::collectHistogram(StuntDouble* sd1, Stu double distance = r12.length(); int whichRBin = distance / deltaR_; - - double cosAngle = evaluateAngle(sd1, sd2); - double halfBin = (nAngleBins_ - 1) * 0.5; - int whichThetaBin = halfBin * (cosAngle + 1.0) - ++histogram_[whichRBin][whichThetaBin]; - - ++npairs_; + if (distance <= len_) { + double cosAngle = evaluateAngle(sd1, sd2); + double halfBin = (nAngleBins_ - 1) * 0.5; + int whichThetaBin = halfBin * (cosAngle + 1.0); + ++histogram_[whichRBin][whichThetaBin]; + + ++npairs_; + } } void GofRAngle::writeRdf() { @@ -120,7 +132,7 @@ void GofRAngle::writeRdf() { double r = deltaR_ * (i + 0.5); for(int j = 0; j < avgGofr_[i].size(); ++j) { - double cosAngle = -1.0 + (i + 0.5)*deltaCosAngle_; + double cosAngle = -1.0 + (j + 0.5)*deltaCosAngle_; rdfStream << r << "\t" << cosAngle << "\t" << avgGofr_[i][j]/nProcessed_ << "\n"; } } @@ -139,14 +151,17 @@ double GofRTheta::evaluateAngle(StuntDouble* sd1, Stun Vector3d r12 = pos1 - pos2; currentSnapshot_->wrapVector(r12); r12.normalize(); - Vector3d dipole = sd1->getElectroFrame().getColumn(2)£» + Vector3d dipole = sd1->getElectroFrame().getColumn(2); dipole.normalize(); - return dot(); + return dot(r12, dipole); } double GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) { Vector3d v1 = sd1->getElectroFrame().getColumn(2); Vector3d v2 = sd1->getElectroFrame().getColumn(2); + v1.normalize(); + v2.normalize(); + return dot(v1, v2); }