diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp
index 13a5aa9ab809515a982e0b7dd03a43a9ae328f43..e9122e7dbd277bf18f3a876b1c833a648855a11e 100644
--- a/src/libnptm/tra_subs.cpp
+++ b/src/libnptm/tra_subs.cpp
@@ -297,7 +297,7 @@ void pwmalp(complex<double> **w, double *up, double *un, complex<double> *ylm, i
int lf = l20 + 1;
int lftl = lf * l20;
double x = 1.0 * lftl;
- complex<double> cl = std::pow((four_pi / sqrt(x)) * uim, 1.0 * l20);
+ complex<double> cl = (four_pi / sqrt(x)) * std::pow(uim, 1.0 * l20);
int mv = l20 + lf;
int m = -lf;
for (int mf20 = 1; mf20 <= mv; mf20++) {
@@ -387,7 +387,7 @@ void wamff(
if (!(onz && onx && ony)) {
rho = sqrt(x * x + y * y);
cph = x / rho;
- sph = y = rho;
+ sph = y / rho;
} else {
if (onz) { // label 10
cph = 1.0;
@@ -429,6 +429,7 @@ void wamff(
}
}
if (lmode == 0 || sth != 0.0) { // label 25
+ bool skip62 = false;
ylm[nlmp - 1] = cc0;
sphar(cth, sth, cph, sph, lm, ylm);
up[0] = cph * cth;
@@ -461,16 +462,19 @@ void wamff(
cfam *= (sth * pmf * ftc1);
for (int i52 = 0; i52 < nlmmt; i52++) wk[i52] = w[i52][0] * cfam;
// returns
+ skip62 = true;
} else if (lmode == 3) { // label 53
ftc2 = 2.0 * cthl / (cthl + cth * rir);
cfam *= (sth * pmf * ftc2);
for (int i55 = 0; i55 < nlmmt; i55++) wk[i55] = w[i55][1] * cfam;
// returns
+ skip62 = true;
}
if (lmode == 0 || lmode == 1) { //label 48
cf1 = cph * ftc1 * cfam;
cf2 = -sph * ftc2 * cfam;
// goes to 62
+ skip62 = false;
}
} else { // label 57
double fam = tra * std::exp(-apfa * apfa) / sqrt(cth);
@@ -479,23 +483,29 @@ void wamff(
double f2 = -sph * fam;
for (int i58 = 0; i58 < nlmmt; i58++) wk[i58] = w[i58][0] * f1 + w[i58][1] * f2;
// returns
+ skip62 = true;
} else if (lmode == 1) { // label 60
cfam = (pmf * sth * fam) * (cph * uim * sph);
cf1 = cph * cfam;
cf2 = -sph * cfam;
// follows on to 62
+ skip62 = false;
} else if (lmode == 2) { // label 65
fam *= (pmf * sth);
for (int i67 = 0; i67 < nlmmt; i67++) wk[i67] = w[i67][0] * fam;
// returns
+ skip62 = true;
} else if (lmode == 3) { // label 68
fam *= (pmf * sth);
for (int i70 = 0; i70 < nlmmt; i70++) wk[i70] = w[i70][1] * fam;
// returns
+ skip62 = true;
}
}
- if (lmode == 0 || lmode == 1) { // label 62
- for (int i63 = 0; i63 < nlmmt; i63++) wk[i63] = w[i63][0] * cf1 + w[i63][1] * cf2;
+ if (!skip62) {
+ if (lmode == 0 || lmode == 1) { // label 62
+ for (int i63 = 0; i63 < nlmmt; i63++) wk[i63] = w[i63][0] * cf1 + w[i63][1] * cf2;
+ }
}
}
// Clean up memory
diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp
index 251ea9117b5100bfc70b6bfccde7e49ef787a1e2..55f6f0e170834c22e3513e063739c0054e9d30b5 100644
--- a/src/trapping/cfrfme.cpp
+++ b/src/trapping/cfrfme.cpp
@@ -49,14 +49,14 @@ void frfme(string data_file, string output_path) {
str_target = m.suffix().str();
regex_search(str_target, m, re);
int jlml = stoi(m.str());
- int lmode, lm, nks, nkv;
- double vk, exri, an, ff, tra;
- double exdc, wp, xip, xi;
- int idfc, nxi;
- double apfafa, pmf, spd, rir, ftcn, fshmx;
- double vxyzmx, delxyz, vknmx, delk, delks;
- double frsh, exril;
- int nlmmt, nrvc;
+ int lmode = 0, lm = 0, nks = 0, nkv = 0;
+ double vk = 0.0, exri = 0.0, an = 0.0, ff = 0.0, tra = 0.0;
+ double exdc = 0.0, wp = 0.0, xip = 0.0, xi = 0.0;
+ int idfc = 0, nxi = 0;
+ double apfafa = 0.0, pmf = 0.0, spd = 0.0, rir = 0.0, ftcn = 0.0, fshmx = 0.0;
+ double vxyzmx = 0.0, delxyz = 0.0, vknmx = 0.0, delk = 0.0, delks = 0.0;
+ double frsh = 0.0, exril = 0.0;
+ int nlmmt = 0, nrvc = 0;
// Vector size variables
int vkzm_size, wsum_size;
// End of vector size variables
@@ -78,6 +78,7 @@ void frfme(string data_file, string output_path) {
tfrfme.read(reinterpret_cast<char *>(&spd), sizeof(double));
tfrfme.read(reinterpret_cast<char *>(&frsh), sizeof(double));
tfrfme.read(reinterpret_cast<char *>(&exril), sizeof(double));
+ vkv = new double[nkv]();
xv = new double[nxv]();
yv = new double[nyv]();
zv = new double[nzv]();
@@ -85,10 +86,9 @@ void frfme(string data_file, string output_path) {
for (int yi = 0; yi < nxv; yi++) tfrfme.read(reinterpret_cast<char *>(&(yv[yi])), sizeof(double));
for (int zi = 0; zi < nxv; zi++) tfrfme.read(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));
fstream temptape2;
- string tempname2 = output_path + "c_TEMPTAPE2";
+ string tempname2 = output_path + "/c_TEMPTAPE2";
temptape2.open(tempname2.c_str(), ios::in | ios::binary);
if (temptape2.is_open()) {
- //vkv = new double[nkv]();
for (int jx = 0; jx < nkv; jx++) temptape2.read(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
vkzm = new double*[nkv];
vkzm_size = nkv;
@@ -228,6 +228,11 @@ void frfme(string data_file, string output_path) {
nlmmt = lm * (lm + 2) * 2;
nks = nksh * 2;
nkv = nks + 1;
+ // Array initialization
+ vkv = new double[nkv]();
+ vkzm = new double*[nkv];
+ for (int vi = 0; vi < nkv; vi++) vkzm[vi] = new double[nkv];
+ // End of array initialization
double vkm = vk * exri;
vknmx = vk * an;
delk = vknmx / nksh;
@@ -239,7 +244,6 @@ void frfme(string data_file, string output_path) {
int nxv = nxs + 1;
int nxshpo = nxsh + 1;
xv = new double[nxv]();
- //xv[nxsh] = 0.0;
for (int i24 = nxshpo; i24 <= nxs; i24++) {
xv[i24] = xv[i24 - 1] + delxyz;
xv[nxv - i24 - 1] = -xv[i24];
@@ -248,7 +252,6 @@ void frfme(string data_file, string output_path) {
int nyv = nys + 1;
int nyshpo = nysh + 1;
yv = new double[nyv]();
- //yv[nysh] = 0.0;
for (int i25 = nyshpo; i25 <= nys; i25++) {
yv[i25] = yv[i25 - 1] + delxyz;
yv[nyv - i25 - 1] = -yv[i25];
@@ -257,18 +260,12 @@ void frfme(string data_file, string output_path) {
int nzv = nzs + 1;
int nzshpo = nzsh + 1;
zv = new double[nzv]();
- //zv[nysh] = 0.0;
for (int i27 = nzshpo; i27 <= nzs; i27++) {
zv[i27] = zv[i27 - 1] + delxyz;
zv[nzv - i27 - 1] = -zv[i27];
} // i27 loop
int nrvc = nxv * nyv * nzv;
int nkshpo = nksh + 1;
- wsum = new complex<double>*[nlmmt];
- wsum_size = nlmmt;
- for (int wsi = 0; wsi < nlmmt; wsi++) wsum[wsi] = new complex<double>[nrvc]();
- vkv = new double[nkv]();
- // vkv[nksh] = 0.0;
for (int i28 = nkshpo; i28 <= nks; i28++) {
vkv[i28] = vkv[i28 - 1] + delk;
vkv[nkv - i28 - 1] = -vkv[i28];
@@ -319,21 +316,34 @@ void frfme(string data_file, string output_path) {
temptape2.write(reinterpret_cast<char *>(&nrvc), sizeof(int));
temptape2.close();
temptape2.open("c_TEMPTAPE2", ios::in | ios::binary);
- vkv = new double[nkv]();
- for (int jx = 0; jx < nkv; jx++)
- temptape2.read(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
- vkzm = new double*[nkv];
- vkzm_size = nkv;
- for (int vki = 0; vki < nkv; vki++) vkzm[vki] = new double[nkv]();
+ for (int jx = 0; jx < nkv; jx++) {
+ double value = 0.0;
+ temptape2.read(reinterpret_cast<char *>(&value), sizeof(double));
+ vkv[jx] = value;
+ }
for (int jy40 = 0; jy40 < nkv; jy40++) {
- for (int jx40 = 0; jx40 < nkv; jx40++)
- temptape2.read(reinterpret_cast<char *>(&(vkzm[jx40][jy40])), sizeof(double));
+ for (int jx40 = 0; jx40 < nkv; jx40++) {
+ double value = 0.0;
+ temptape2.read(reinterpret_cast<char *>(&value), sizeof(double));
+ vkzm[jx40][jy40] = value;
+ }
} // jy40 loop
temptape2.close();
- wk = new complex<double>[nlmmt];
- w = new complex<double>*[nkv];
- for (int wi = 0; wi < nkv; wi++) w[wi] = new complex<double>[nkv]();
+ if (wsum != NULL) {
+ for (int wsi = wsum_size - 1; wsi > -1; wsi--) delete[] wsum[wsi];
+ delete[] wsum;
+ }
+ wsum = new complex<double>*[nlmmt];
for (int j80 = jlmf - 1; j80 < jlml; j80++) {
+ // w matrix
+ if (w != NULL) {
+ for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
+ delete[] w;
+ }
+ w = new complex<double>*[nkv];
+ for (int wi = 0; wi < nkv; wi++) w[wi] = new complex<double>[nkv]();
+ if (wk != NULL) delete[] wk;
+ wk = new complex<double>[nlmmt]();
temptape1.open(temp_name1.c_str(), ios::in | ios::binary);
for (int jy50 = 0; jy50 < nkv; jy50++) {
for (int jx50 = 0; jx50 < nkv; jx50++) {
@@ -348,6 +358,7 @@ void frfme(string data_file, string output_path) {
} // jy50 loop
temptape1.close();
int ixyz = 0;
+ wsum[j80] = new complex<double>[nrvc]();
for (int iz75 = 0; iz75 < nzv; iz75++) {
double z = zv[iz75] + frsh;
for (int iy70 = 0; iy70 < nyv; iy70++) {
@@ -364,11 +375,11 @@ void frfme(string data_file, string output_path) {
double vkzl = vkzm[nkv - 1][jy60];
complex<double> phasl = exp(uim * (vkx * x + vky * y + vkzl * z));
complex<double> sumx = 0.5 * (w[0][jy60] * phasf + w[nkv - 1][jy60] * phasl);
- for (int jx55 = 1; jx55 < nks; jx55++) {
- vkx = vkv[jx55];
- double vkz = vkzm[jx55][jy60];
+ for (int jx55 = 2; jx55 <= nks; jx55++) {
+ vkx = vkv[jx55 - 1];
+ double vkz = vkzm[jx55 - 1][jy60];
complex<double> phas = exp(uim * (vkx * x + vky * y + vkz * z));
- sumx += (w[jx55][jy60] * phas);
+ sumx += (w[jx55 - 1][jy60] * phas);
} // jx55 loop
if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5;
sumy += sumx;
@@ -439,14 +450,14 @@ void frfme(string data_file, string output_path) {
for (int vki = vkzm_size - 1; vki > -1; vki--) delete[] vkzm[vki];
delete[] vkzm;
}
+ if (w != NULL) {
+ for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
+ delete[] w;
+ }
if (wsum != NULL) {
for (int wsi = wsum_size - 1; wsi > -1; wsi--) delete[] wsum[wsi];
delete[] wsum;
}
if (wk != NULL) delete[] wk;
- if (w != NULL) {
- for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
- delete[] w;
- }
printf("Done.\n");
}