ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/applications/staticProps/GofRAngle.cpp
(Generate patch)

Comparing trunk/src/applications/staticProps/GofRAngle.cpp (file contents):
Revision 309 by tim, Thu Feb 10 14:15:52 2005 UTC vs.
Revision 328 by tim, Sun Feb 13 20:36:24 2005 UTC

# Line 46 | Line 46 | namespace oopse {
46  
47   namespace oopse {
48  
49 < GofRAngle::GofRAngle(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2, double len)
50 <    : RadialDistrFunc(info, filename, sele1, sele2, len){
49 > GofRAngle::GofRAngle(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2)
50 >    : RadialDistrFunc(info, filename, sele1, sele2){
51  
52 –    histogram_.resize(nbins_);
53 –    avgGofr_.resize(nbins_);
52   }
53  
54  
55   void GofRAngle::preProcess() {
56 <    avgGofr_.resize(nbins_);
57 <    for (int i = 0; i < avgGofr_.size(); ++i)
56 >
57 >    for (int i = 0; i < avgGofr_.size(); ++i) {
58          std::fill(avgGofr_[i].begin(), avgGofr_[i].end(), 0);
59 +    }
60   }
61  
62   void GofRAngle::initalizeHistogram() {
# Line 71 | Line 70 | void GofRAngle::processHistogram() {
70  
71      double volume = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
72      double pairDensity = npairs_ /volume;
73 <    double pairConstant = ( 4.0 * PI * pairDensity ) / 3.0;
73 >    double pairConstant = ( 4.0 * NumericConstant::PI * pairDensity ) / 3.0;
74  
75      for(int i = 0 ; i < histogram_.size(); ++i){
76  
# Line 99 | Line 98 | void GofRAngle::collectHistogram(StuntDouble* sd1, Stu
98      currentSnapshot_->wrapVector(r12);
99  
100      double distance = r12.length();
101 <    int whichBin = distance / deltaR_;
101 >    int whichRBin = distance / deltaR_;
102  
103 <    
104 <    double cosAngle = evaluateAngle(sd1, sd2);
105 <    histogram_[whichBin] ++;
106 <    
107 <    npairs_++;
103 >    if (distance <= len_) {
104 >        double cosAngle = evaluateAngle(sd1, sd2);
105 >        double halfBin = (nAngleBins_ - 1) * 0.5;
106 >        int whichThetaBin = halfBin * (cosAngle + 1.0);
107 >        ++histogram_[whichRBin][whichThetaBin];
108 >        
109 >        ++npairs_;
110 >    }
111   }
112  
113   void GofRAngle::writeRdf() {
# Line 119 | Line 121 | void GofRAngle::writeRdf() {
121              double r = deltaR_ * (i + 0.5);
122  
123              for(int j = 0; j < avgGofr_[i].size(); ++j) {
124 <                double cosAngle = ;
124 >                double cosAngle = -1.0 + (i + 0.5)*deltaCosAngle_;
125                  rdfStream << r << "\t" << cosAngle << "\t" << avgGofr_[i][j]/nProcessed_ << "\n";
126              }
127          }
# Line 132 | Line 134 | void GofRAngle::writeRdf() {
134      rdfStream.close();
135   }
136  
135 –
136 –
137   double GofRTheta::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) {
138      Vector3d pos1 = sd1->getPos();
139      Vector3d pos2 = sd2->getPos();
140      Vector3d r12 = pos1 - pos2;
141      currentSnapshot_->wrapVector(r12);
142      r12.normalize();
143 <    Vector3d dipole = sd1->getElectroFrame().getColumn(2)£»
143 >    Vector3d dipole = sd1->getElectroFrame().getColumn(2);
144      dipole.normalize();    
145 <    return dot();
145 >    return dot(r12, dipole);
146   }
147  
148   double GofROmega::evaluateAngle(StuntDouble* sd1, StuntDouble* sd2) {
149      Vector3d v1 = sd1->getElectroFrame().getColumn(2);
150 <    Vector3d v2 = sd1->getElectroFrame().getColumn(2);
151 <
152 <    
150 >    Vector3d v2 = sd1->getElectroFrame().getColumn(2);    
151 >    v1.normalize();
152 >    v2.normalize();
153 >    return dot(v1, v2);
154   }
155  
156  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines