| 45 |
|
|
| 46 |
|
namespace oopse { |
| 47 |
|
LegendreCorrFunc::LegendreCorrFunc(SimInfo* info, const std::string& filename, const std::string& sele1, const std::string& sele2, int order) |
| 48 |
< |
: ParticleTimeCorrFunc(info, filename, sele1, sele2, DataStorage::dslElectroFrame){ |
| 48 |
> |
: ParticleTimeCorrFunc(info, filename, sele1, sele2, DataStorage::dslAmat | DataStorage::dslElectroFrame){ |
| 49 |
|
|
| 50 |
|
setCorrFuncType("Legendre Correlation Function"); |
| 51 |
|
setOutputName(getPrefix(dumpFilename_) + ".lcorr"); |
| 52 |
+ |
histogram_.resize(nTimeBins_); |
| 53 |
+ |
count_.resize(nTimeBins_); |
| 54 |
+ |
nSelected_ = seleMan1_.getSelectionCount(); |
| 55 |
+ |
std::cout << "sele1 = " << sele1 << "\n"; |
| 56 |
+ |
std::cout << "sele2 = " << sele2 << "\n"; |
| 57 |
+ |
std::cout << "nSelected = " << nSelected_ << "\n"; |
| 58 |
+ |
assert( nSelected_ == seleMan2_.getSelectionCount()); |
| 59 |
+ |
|
| 60 |
|
LegendrePolynomial polynomial(order); |
| 61 |
|
legendre_ = polynomial.getLegendrePolynomial(order); |
| 62 |
< |
|
| 62 |
> |
} |
| 63 |
|
|
| 64 |
+ |
void LegendreCorrFunc::correlateFrames(int frame1, int frame2) { |
| 65 |
+ |
Snapshot* snapshot1 = bsMan_->getSnapshot(frame1); |
| 66 |
+ |
Snapshot* snapshot2 = bsMan_->getSnapshot(frame2); |
| 67 |
+ |
assert(snapshot1 && snapshot2); |
| 68 |
+ |
|
| 69 |
+ |
RealType time1 = snapshot1->getTime(); |
| 70 |
+ |
RealType time2 = snapshot2->getTime(); |
| 71 |
+ |
|
| 72 |
+ |
int timeBin = int ((time2 - time1) /deltaTime_ + 0.5); |
| 73 |
+ |
count_[timeBin] += nSelected_; |
| 74 |
+ |
|
| 75 |
+ |
int i; |
| 76 |
+ |
int j; |
| 77 |
+ |
StuntDouble* sd1; |
| 78 |
+ |
StuntDouble* sd2; |
| 79 |
+ |
|
| 80 |
+ |
for (sd1 = seleMan1_.beginSelected(i), sd2 = seleMan2_.beginSelected(j); |
| 81 |
+ |
sd1 != NULL && sd2 != NULL; |
| 82 |
+ |
sd1 = seleMan1_.nextSelected(i), sd2 = seleMan2_.nextSelected(j)) { |
| 83 |
+ |
|
| 84 |
+ |
Vector3d corrVals = calcCorrVals(frame1, frame2, sd1, sd2); |
| 85 |
+ |
histogram_[timeBin] += corrVals; |
| 86 |
|
} |
| 87 |
+ |
|
| 88 |
+ |
} |
| 89 |
|
|
| 90 |
< |
RealType LegendreCorrFunc::calcCorrVal(int frame1, int frame2, StuntDouble* sd1, StuntDouble* sd2) { |
| 91 |
< |
Vector3d v1 =sd1->getA(frame1).getColumn(2); |
| 92 |
< |
Vector3d v2 = sd2->getA(frame2).getColumn(2); |
| 90 |
> |
void LegendreCorrFunc::postCorrelate() { |
| 91 |
> |
for (int i =0 ; i < nTimeBins_; ++i) { |
| 92 |
> |
if (count_[i] > 0) { |
| 93 |
> |
histogram_[i] /= count_[i]; |
| 94 |
> |
} |
| 95 |
> |
} |
| 96 |
> |
} |
| 97 |
|
|
| 98 |
< |
return legendre_.evaluate(dot(v1, v2)/(v1.length()*v2.length())); |
| 98 |
> |
void LegendreCorrFunc::preCorrelate() { |
| 99 |
> |
// Fill the histogram with empty Vector3d: |
| 100 |
> |
std::fill(histogram_.begin(), histogram_.end(), Vector3d(0.0)); |
| 101 |
> |
// count array set to zero |
| 102 |
> |
std::fill(count_.begin(), count_.end(), 0); |
| 103 |
|
} |
| 104 |
|
|
| 105 |
|
|
| 106 |
+ |
|
| 107 |
+ |
Vector3d LegendreCorrFunc::calcCorrVals(int frame1, int frame2, StuntDouble* sd1, StuntDouble* sd2) { |
| 108 |
+ |
Vector3d v1x = sd1->getA(frame1).getColumn(0); |
| 109 |
+ |
Vector3d v2x = sd2->getA(frame2).getColumn(0); |
| 110 |
+ |
Vector3d v1y = sd1->getA(frame1).getColumn(1); |
| 111 |
+ |
Vector3d v2y = sd2->getA(frame2).getColumn(1); |
| 112 |
+ |
Vector3d v1z = sd1->getA(frame1).getColumn(2); |
| 113 |
+ |
Vector3d v2z = sd2->getA(frame2).getColumn(2); |
| 114 |
+ |
|
| 115 |
+ |
RealType uxprod = legendre_.evaluate(dot(v1x, v2x)/(v1x.length()*v2x.length())); |
| 116 |
+ |
RealType uyprod = legendre_.evaluate(dot(v1y, v2y)/(v1y.length()*v2y.length())); |
| 117 |
+ |
RealType uzprod = legendre_.evaluate(dot(v1z, v2z)/(v1z.length()*v2z.length())); |
| 118 |
+ |
|
| 119 |
+ |
return Vector3d(uxprod, uyprod, uzprod); |
| 120 |
+ |
|
| 121 |
+ |
} |
| 122 |
+ |
|
| 123 |
+ |
|
| 124 |
|
void LegendreCorrFunc::validateSelection(const SelectionManager& seleMan) { |
| 125 |
|
StuntDouble* sd; |
| 126 |
|
int i; |
| 127 |
|
for (sd = seleMan1_.beginSelected(i); sd != NULL; sd = seleMan1_.nextSelected(i)) { |
| 128 |
|
if (!sd->isDirectionalAtom()) { |
| 129 |
|
sprintf(painCave.errMsg, |
| 130 |
< |
"LegendreCorrFunc::validateSelection Error: selected atoms do not Directional\n"); |
| 130 |
> |
"LegendreCorrFunc::validateSelection Error: selected atoms are not Directional\n"); |
| 131 |
|
painCave.isFatal = 1; |
| 132 |
|
simError(); |
| 133 |
|
} |
| 135 |
|
|
| 136 |
|
} |
| 137 |
|
|
| 138 |
< |
} |
| 138 |
> |
void LegendreCorrFunc::writeCorrelate() { |
| 139 |
> |
std::ofstream ofs(getOutputFileName().c_str()); |
| 140 |
|
|
| 141 |
+ |
if (ofs.is_open()) { |
| 142 |
|
|
| 143 |
+ |
ofs << "#" << getCorrFuncType() << "\n"; |
| 144 |
+ |
ofs << "#time\tPn(costheta_x)\tPn(costheta_y)\tPn(costheta_z)\n"; |
| 145 |
|
|
| 146 |
+ |
for (int i = 0; i < nTimeBins_; ++i) { |
| 147 |
+ |
ofs << time_[i] << "\t" << |
| 148 |
+ |
histogram_[i](0) << "\t" << |
| 149 |
+ |
histogram_[i](1) << "\t" << |
| 150 |
+ |
histogram_[i](2) << "\t" << "\n"; |
| 151 |
+ |
} |
| 152 |
+ |
|
| 153 |
+ |
} else { |
| 154 |
+ |
sprintf(painCave.errMsg, |
| 155 |
+ |
"LegendreCorrFunc::writeCorrelate Error: fail to open %s\n", getOutputFileName().c_str()); |
| 156 |
+ |
painCave.isFatal = 1; |
| 157 |
+ |
simError(); |
| 158 |
+ |
} |
| 159 |
+ |
ofs.close(); |
| 160 |
+ |
} |
| 161 |
+ |
} |
