| 46 |  |  | 
| 47 |  | namespace oopse { | 
| 48 |  |  | 
| 49 | < | GofRAngle::GofRAngle(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2) | 
| 50 | < | : RadialDistrFunc(info, filename, sele1, sele2){ | 
| 49 | > | GofRAngle::GofRAngle(SimInfo* info, const std::string& filename, const std::string& sele1, | 
| 50 | > | const std::string& sele2, double len, int nrbins, int nangleBins) | 
| 51 | > | : RadialDistrFunc(info, filename, sele1, sele2), len_(len), nRBins_(nrbins), nAngleBins_(nangleBins){ | 
| 52 |  |  | 
| 53 | + | deltaR_ = len_ /nRBins_; | 
| 54 | + | deltaCosAngle_ = 2.0 / nAngleBins_; | 
| 55 | + |  | 
| 56 | + | histogram_.resize(nRBins_); | 
| 57 | + | avgGofr_.resize(nRBins_); | 
| 58 | + | for (int i = 0 ; i < nRBins_; ++i) { | 
| 59 | + | histogram_[i].resize(nAngleBins_); | 
| 60 | + | avgGofr_[i].resize(nAngleBins_); | 
| 61 | + | } | 
| 62 |  | } | 
| 63 |  |  | 
| 64 |  |  | 
| 78 |  |  | 
| 79 |  | void GofRAngle::processHistogram() { | 
| 80 |  |  | 
| 81 | + | int nPairs = getNPairs(); | 
| 82 |  | double volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume(); | 
| 83 | < | double pairDensity = npairs_ /volume; | 
| 84 | < | double pairConstant = ( 4.0 * PI * pairDensity ) / 3.0; | 
| 83 | > | double pairDensity = nPairs /volume; | 
| 84 | > | double pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0; | 
| 85 |  |  | 
| 86 |  | for(int i = 0 ; i < histogram_.size(); ++i){ | 
| 87 |  |  | 
| 111 |  | double distance = r12.length(); | 
| 112 |  | int whichRBin = distance / deltaR_; | 
| 113 |  |  | 
| 114 | < |  | 
| 115 | < | double cosAngle = evaluateAngle(sd1, sd2); | 
| 116 | < | double halfBin = (nAngleBins_ - 1) * 0.5; | 
| 117 | < | int whichThetaBin = halfBin * (cosAngle + 1.0) | 
| 118 | < | ++histogram_[whichRBin][whichThetaBin]; | 
| 119 | < |  | 
| 120 | < | ++npairs_; | 
| 114 | > | if (distance <= len_) { | 
| 115 | > | double cosAngle = evaluateAngle(sd1, sd2); | 
| 116 | > | double halfBin = (nAngleBins_ - 1) * 0.5; | 
| 117 | > | int whichThetaBin = halfBin * (cosAngle + 1.0); | 
| 118 | > | ++histogram_[whichRBin][whichThetaBin]; | 
| 119 | > |  | 
| 120 | > | ++npairs_; | 
| 121 | > | } | 
| 122 |  | } | 
| 123 |  |  | 
| 124 |  | void GofRAngle::writeRdf() { | 
| 151 |  | Vector3d r12 = pos1 - pos2; | 
| 152 |  | currentSnapshot_->wrapVector(r12); | 
| 153 |  | r12.normalize(); | 
| 154 | < | Vector3d dipole = sd1->getElectroFrame().getColumn(2)£» | 
| 154 | > | Vector3d dipole = sd1->getElectroFrame().getColumn(2); | 
| 155 |  | dipole.normalize(); | 
| 156 | < | return dot(); | 
| 156 | > | return dot(r12, dipole); | 
| 157 |  | } | 
| 158 |  |  | 
| 159 |  | double GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) { | 
| 160 |  | Vector3d v1 = sd1->getElectroFrame().getColumn(2); | 
| 161 |  | Vector3d v2 = sd1->getElectroFrame().getColumn(2); | 
| 162 | + | v1.normalize(); | 
| 163 | + | v2.normalize(); | 
| 164 | + | return dot(v1, v2); | 
| 165 |  | } | 
| 166 |  |  | 
| 167 |  |  |