85 |
|
MinW_ = -0.25; |
86 |
|
MaxW_ = 0.25; |
87 |
|
deltaW_ = (MaxW_ - MinW_) / nbins; |
88 |
– |
} |
88 |
|
|
89 |
+ |
// Make arrays for Wigner3jm |
90 |
+ |
double* THRCOF = new double[2*lMax_+1]; |
91 |
+ |
// Variables for Wigner routine |
92 |
+ |
double lPass, m1Pass, m2m, m2M; |
93 |
+ |
int error, mSize; |
94 |
+ |
mSize = 2*lMax_+1; |
95 |
+ |
|
96 |
+ |
for (int l = 0; l <= lMax_; l++) { |
97 |
+ |
lPass = (double)l; |
98 |
+ |
for (int m1 = -l; m1 <= l; m1++) { |
99 |
+ |
m1Pass = (double)m1; |
100 |
+ |
|
101 |
+ |
std::pair<int,int> lm = std::make_pair(l, m1); |
102 |
+ |
|
103 |
+ |
// Zero work array |
104 |
+ |
for (int ii = 0; ii < 2*l + 1; ii++){ |
105 |
+ |
THRCOF[ii] = 0.0; |
106 |
+ |
} |
107 |
+ |
|
108 |
+ |
// Get Wigner coefficients |
109 |
+ |
Wigner3jm(&lPass, &lPass, &lPass, |
110 |
+ |
&m1Pass, &m2m, &m2M, |
111 |
+ |
THRCOF, &mSize, &error); |
112 |
+ |
|
113 |
+ |
m2Min[lm] = (int)floor(m2m); |
114 |
+ |
m2Max[lm] = (int)floor(m2M); |
115 |
+ |
|
116 |
+ |
for (int mmm = 0; mmm < (int)(m2M - m2m); mmm++) { |
117 |
+ |
w3j[lm].push_back(THRCOF[mmm]); |
118 |
+ |
} |
119 |
+ |
} |
120 |
+ |
} |
121 |
+ |
} |
122 |
+ |
|
123 |
|
BondOrderParameter::~BondOrderParameter() { |
124 |
|
Q_histogram_.clear(); |
125 |
|
W_histogram_.clear(); |
161 |
|
int nBonds, Nbonds; |
162 |
|
SphericalHarmonic sphericalHarmonic; |
163 |
|
int i, j; |
131 |
– |
// Make arrays for Wigner3jm |
132 |
– |
double* THRCOF = new double[2*lMax_+1]; |
133 |
– |
// Variables for Wigner routine |
134 |
– |
double lPass, m1Pass, m2Min, m2Max; |
135 |
– |
int error, m1, m2, m3, mSize; |
136 |
– |
mSize = 2*lMax_+1; |
164 |
|
|
165 |
|
DumpReader reader(info_, dumpFilename_); |
166 |
|
int nFrames = reader.getNFrames(); |
272 |
|
|
273 |
|
for (int l = 0; l <= lMax_; l++) { |
274 |
|
w[l] = 0.0; |
248 |
– |
lPass = (double)l; |
275 |
|
for (int m1 = -l; m1 <= l; m1++) { |
276 |
< |
// Zero work array |
277 |
< |
for (int ii = 0; ii < 2*l + 1; ii++){ |
278 |
< |
THRCOF[ii] = 0.0; |
279 |
< |
} |
280 |
< |
// Get Wigner coefficients |
281 |
< |
m1Pass = (double)m1; |
256 |
< |
|
257 |
< |
Wigner3jm(&lPass, &lPass, &lPass, |
258 |
< |
&m1Pass, &m2Min, &m2Max, |
259 |
< |
THRCOF, &mSize, &error); |
260 |
< |
|
261 |
< |
for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) { |
262 |
< |
m2 = (int)floor(m2Min) + mmm; |
263 |
< |
m3 = -m1-m2; |
264 |
< |
w[l] += THRCOF[mmm] * |
265 |
< |
q[std::make_pair(l,m1)] * |
266 |
< |
q[std::make_pair(l,m2)] * |
267 |
< |
q[std::make_pair(l,m3)]; |
276 |
> |
std::pair<int,int> lm = std::make_pair(l, m1); |
277 |
> |
for (int mmm = 0; mmm < (m2Max[lm] - m2Min[lm]); mmm++) { |
278 |
> |
int m2 = m2Min[lm] + mmm; |
279 |
> |
int m3 = -m1-m2; |
280 |
> |
w[l] += w3j[lm][mmm] * q[lm] * |
281 |
> |
q[std::make_pair(l,m2)] * q[std::make_pair(l,m3)]; |
282 |
|
} |
283 |
|
} |
284 |
|
|
317 |
|
|
318 |
|
for (int l = 0; l <= lMax_; l++) { |
319 |
|
W[l] = 0.0; |
306 |
– |
lPass = (double)l; |
320 |
|
for (int m1 = -l; m1 <= l; m1++) { |
321 |
< |
// Zero work array |
322 |
< |
for (int ii = 0; ii < 2*l + 1; ii++){ |
323 |
< |
THRCOF[ii] = 0.0; |
321 |
> |
std::pair<int,int> lm = std::make_pair(l, m1); |
322 |
> |
for (int mmm = 0; mmm < (m2Max[lm] - m2Min[lm]); mmm++) { |
323 |
> |
int m2 = m2Min[lm] + mmm; |
324 |
> |
int m3 = -m1-m2; |
325 |
> |
W[l] += w3j[lm][mmm] * QBar[lm] * |
326 |
> |
QBar[std::make_pair(l,m2)] * QBar[std::make_pair(l,m3)]; |
327 |
|
} |
312 |
– |
// Get Wigner coefficients |
313 |
– |
m1Pass = (double)m1; |
314 |
– |
|
315 |
– |
Wigner3jm(&lPass, &lPass, &lPass, |
316 |
– |
&m1Pass, &m2Min, &m2Max, |
317 |
– |
THRCOF, &mSize, &error); |
318 |
– |
|
319 |
– |
for (int mmm = 0; mmm < (int)(m2Max - m2Min); mmm++) { |
320 |
– |
m2 = (int)floor(m2Min) + mmm; |
321 |
– |
m3 = -m1-m2; |
322 |
– |
W[l] += THRCOF[mmm] * |
323 |
– |
QBar[std::make_pair(l,m1)] * |
324 |
– |
QBar[std::make_pair(l,m2)] * |
325 |
– |
QBar[std::make_pair(l,m3)]; |
326 |
– |
} |
328 |
|
} |
329 |
|
|
330 |
|
W_hat[l] = W[l] / pow(Q2[l], 1.5); |