diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index 780c078c3dfede172593736f0e5e5d0b14d06baa..bf3b108a9d2ada949446800aa4d09c6721c2cbfe 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -7,11 +7,11 @@
#endif
//Structures for TRAPPING (that will probably move to Commons)
-struct cil {
+struct CIL {
int le, nlem, nlemt, mxmpo, mxim;
};
-struct ccr {
+struct CCR {
double cof, cimu;
};
//End of TRAPPING structures
@@ -286,3 +286,40 @@ void frfmer(
} // jy90 loop
delete[] wk;
}
+
+/*! C++ porting of SAMPOA
+ *
+ * \param ac: Vector of complex. QUESTION: definition?
+ * \param tms: Matrix of complex. QUESTION: definition?
+ * \param ws: Vector of complex. QUESTION: definition?
+ * \param cil: `CIL *` Pointer to a CIL structure.
+ *
+ * This function builds the AC vector using TMS and WS.
+ */
+void sampoa(
+ std::complex<double> *ac, std::complex<double> **tms, std::complex<double> *ws,
+ CIL *cil
+) {
+ std::complex<double> **tm = new std::complex<double>*[2];
+ tm[0] = new std::complex<double>[2]();
+ tm[1] = new std::complex<double>[2]();
+ int i = 0;
+ ac = new std::complex<double>[cil->nlemt]();
+ for (int l20 = 0; l20 < cil->le; l20++) {
+ tm[0][0] = tms[l20][0];
+ tm[0][1] = tms[l20][1];
+ tm[1][1] = tms[l20][2];
+ tm[1][0] = tm[0][1];
+ int l = l20 + 1;
+ int ltpo = l + l + 1;
+ for (int im20 = 0; im20 < ltpo; im20++) {
+ int ie = i + cil->nlem;
+ ac[i] = tm[0][0] * ws[i] + tm[0][1] * ws[ie];
+ ac[ie] = tm[1][0] * ws[i] + tm[1][1] * ws[ie];
+ i++;
+ } // im20 loop
+ } // l20 loop
+ delete[] tm[1];
+ delete[] tm[0];
+ delete[] tm;
+}
diff --git a/src/trapping/lffft.cpp b/src/trapping/lffft.cpp
index 500f4230139f5e82846af1b1d2b76c8e60f12272..33d516a5ae5929d5154f4d0821072751986f31de 100644
--- a/src/trapping/lffft.cpp
+++ b/src/trapping/lffft.cpp
@@ -22,14 +22,15 @@ using namespace std;
*/
int main() {
const complex<double> uim(0.0, 1.0);
+ const double sq2i = 1.0 / sqrt(2.0);
+ const complex<double> sq2iti = sq2i * uim;
+
+ fstream tlfff, tlfft;
string data_file = "../../test_data/trapping/DLFFFT";
string output_path = ".";
- double ****zpv = NULL, *ac = NULL;
- complex<double> *ws = NULL;
- double *fffe = NULL, *fffs = NULL;
- double *ffte = NULL, *ffts = NULL;
+ double ****zpv = NULL;
double *xv = NULL, *yv = NULL, *zv = NULL;
- complex<double> *wsl = NULL;
+ complex<double> *ac = NULL, *ws = NULL, *wsl = NULL;
complex<double> **am0m = NULL;
complex<double> **amd = NULL;
int **indam = NULL;
@@ -37,8 +38,8 @@ int main() {
int jft, jss, jtw;
int is, le, nvam = 0;
double vks, exris;
- cil *p_cil = new cil();
- ccr *p_ccr = new ccr();
+ CIL *cil = new CIL();
+ CCR *ccr = new CCR();
int num_lines = 0;
string *file_lines = load_file(data_file, &num_lines);
@@ -60,9 +61,9 @@ int main() {
ttms.read(reinterpret_cast<char *>(&le), sizeof(int));
ttms.read(reinterpret_cast<char *>(&vks), sizeof(double));
ttms.read(reinterpret_cast<char *>(&exris), sizeof(double));
- p_cil->le = le;
- p_cil->nlem = le * (le + 2);
- p_cil->nlemt = p_cil->nlem + p_cil->nlem;
+ cil->le = le;
+ cil->nlem = le * (le + 2);
+ cil->nlemt = cil->nlem + cil->nlem;
if (is > 2222) { // label 120
tms = new complex<double>*[le];
for (int ti = 0; ti < le; ti++) tms[ti] = new complex<double>[3]();
@@ -93,10 +94,10 @@ int main() {
tmse[i] = complex<double>(vreal, vimag);
} // i loop
} else if (is >= 0) { // label 135
- am0m = new complex<double>*[p_cil->nlemt];
- for (int ai = 0; ai < p_cil->nlemt; ai++) am0m[ai] = new complex<double>[p_cil->nlemt]();
- for (int i = 0; i < p_cil->nlemt; i++) {
- for (int j = 0; j < p_cil->nlemt; j++) {
+ am0m = new complex<double>*[cil->nlemt];
+ for (int ai = 0; ai < cil->nlemt; ai++) am0m[ai] = new complex<double>[cil->nlemt]();
+ for (int i = 0; i < cil->nlemt; i++) {
+ for (int j = 0; j < cil->nlemt; j++) {
double vreal, vimag;
ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
@@ -125,8 +126,8 @@ int main() {
} // j loop
} // i loop
ttms.read(reinterpret_cast<char *>(&vint), sizeof(int));
- p_cil->mxmpo = vint;
- p_cil->mxim = vint * 2 - 1;
+ cil->mxmpo = vint;
+ cil->mxim = vint * 2 - 1;
}
// label 150
ttms.close();
@@ -174,11 +175,10 @@ int main() {
thdps(le, zpv);
double exdc = exri * exri;
double sqk = vk * vk * exdc;
- p_ccr->cof = 1.0 / sqk;
- p_ccr->cimu = p_ccr->cof / sqrt(2.0);
+ ccr->cof = 1.0 / sqk;
+ ccr->cimu = ccr->cof / sqrt(2.0);
if (jss != 1) {
string tlfff_name = output_path + "/c_TLFFF";
- fstream tlfff;
tlfff.open(tlfff_name.c_str(), ios::out | ios::binary);
if (tlfff.is_open()) {
tlfff.write(reinterpret_cast<char *>(&lmode), sizeof(int));
@@ -209,13 +209,121 @@ int main() {
}
// label 155
if (!goto160) {
- const double sq2i = 1.0 / sqrt(2.0);
- const complex<double> sq2iti = sq2i * uim;
if (jss != 1) {
// Would open the ITT file.
+ string tlfft_name = output_path + "/c_TLFFT";
+ tlfft.open(tlfft_name.c_str(), ios::out | ios::binary);
+ if (tlfft.is_open()) {
+ tlfft.write(reinterpret_cast<char *>(&lmode), sizeof(int));
+ tlfft.write(reinterpret_cast<char *>(&le), sizeof(int));
+ tlfft.write(reinterpret_cast<char *>(&nkv), sizeof(int));
+ tlfft.write(reinterpret_cast<char *>(&nxv), sizeof(int));
+ tlfft.write(reinterpret_cast<char *>(&nyv), sizeof(int));
+ tlfft.write(reinterpret_cast<char *>(&nzv), sizeof(int));
+ tlfft.write(reinterpret_cast<char *>(&vk), sizeof(double));
+ tlfft.write(reinterpret_cast<char *>(&exri), sizeof(double));
+ tlfft.write(reinterpret_cast<char *>(&an), sizeof(double));
+ tlfft.write(reinterpret_cast<char *>(&ff), sizeof(double));
+ tlfft.write(reinterpret_cast<char *>(&tra), sizeof(double));
+ tlfft.write(reinterpret_cast<char *>(&spd), sizeof(double));
+ tlfft.write(reinterpret_cast<char *>(&frsh), sizeof(double));
+ tlfft.write(reinterpret_cast<char *>(&exril), sizeof(double));
+ for (int i = 0; i < nxv; i++)
+ tlfft.write(reinterpret_cast<char *>(&(xv[i])), sizeof(double));
+ for (int i = 0; i < nyv; i++)
+ tlfft.write(reinterpret_cast<char *>(&(yv[i])), sizeof(double));
+ for (int i = 0; i < nzv; i++)
+ tlfft.write(reinterpret_cast<char *>(&(zv[i])), sizeof(double));
+ } else { // Should never happen.
+ printf("ERROR: could not open TLFFT file.\n");
+ }
}
}
// label 160
+ const int nlmm = lm * (lm + 2);
+ const int nlmmt = nlmm + nlmm;
+ ws = new complex<double>[nlmmt]();
+ if (lm > le) wsl = new complex<double>[nlmmt]();
+ for (int iz475 = 0; iz475 < nzv; iz475++) {
+ for (int iy475 = 0; iy475 < nyv; iy475++) {
+ for (int ix475 = 0; ix475 < nxv; ix475++) {
+ for (int i = 0; i < nlmmt; i++) {
+ double vreal, vimag;
+ binary_input.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ binary_input.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ if (lm <= le) {
+ ws[i] = complex<double>(vreal, vimag);
+ } else { // label 170
+ wsl[i] = complex<double>(vreal, vimag);
+ for (int i175 = 0; i175 < cil->nlem; i175++) {
+ int ie = i175 + cil->nlem;
+ int iel = i175 + nlmm;
+ ws[i175] = wsl[i175];
+ ws[ie] = wsl[iel];
+ } // i175 loop
+ }
+ } // i loop
+ // label 180
+ bool goto475 = false;
+ if (is != 2222) {
+ if (is != 1111) {
+ if (is > 0) { // Goes to 305
+ // Would call CAMP(AC,AM0M,WS)
+ } else if (is < 0) { // Goes to 405
+ // Would call CZAMP(AC,AMD,INDAM,WS)
+ }
+ } else {
+ // Would call SAMP(AC,TMSM,TMSE,WS)
+ }
+ } else {
+ sampoa(ac, tms, ws, cil);
+ }
+ // label 445
+ if (jft <= 0) {
+ double *fffe = new double[3]();
+ double *fffs = new double[3]();
+ // Would call FFRF(ZPV,AC,WS,FFFE,FFFS)
+ if (jss == 1) {
+ // Writes to 66
+ } else { // label 450
+ for (int i = 0; i < 3; i++) {
+ double value = fffe[i] - fffs[i];
+ tlfff.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ if (jtw == 1) {
+ // Writes to 66
+ }
+ }
+ if (jft < 0) goto475 = true;
+ delete[] fffe;
+ delete[] fffs;
+ }
+ // label 460
+ if (!goto475) {
+ double *ffte = new double[3]();
+ double *ffts = new double[3]();
+ // Would call FFRT(AC,WS,FFTE,FFTS)
+ if (jss == 1) {
+ // Writes to 67
+ } else { // label 470
+ for (int i = 0; i < 3; i++) {
+ double value = ffte[i] - ffts[i];
+ tlfft.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ if (jtw == 1) {
+ // Writes to 67
+ }
+ }
+ delete[] ffte;
+ delete[] ffts;
+ }
+ } // ix475 loop
+ } // iy475 loop
+ } // iz475 loop
+ if (jss != 1) {
+ if (jft <= 0) tlfff.close();
+ if (jft >= 0) tlfft.close();
+ }
}
}
binary_input.close();
@@ -229,10 +337,6 @@ int main() {
// Clean up memory
if (ac != NULL) delete[] ac;
if (ws != NULL) delete[] ws;
- if (fffe != NULL) delete[] fffe;
- if (ffte != NULL) delete[] ffte;
- if (fffs != NULL) delete[] fffs;
- if (ffts != NULL) delete[] ffts;
if (xv != NULL) delete[] xv;
if (yv != NULL) delete[] yv;
if (zv != NULL) delete[] zv;
@@ -244,7 +348,7 @@ int main() {
delete[] tms;
}
if (am0m != NULL) {
- for (int ai = p_cil->nlemt - 1; ai > -1; ai--) delete[] am0m[ai];
+ for (int ai = cil->nlemt - 1; ai > -1; ai--) delete[] am0m[ai];
delete[] am0m;
}
if (amd != NULL) {
@@ -265,8 +369,8 @@ int main() {
} // zi loop
delete[] zpv;
}
- delete p_cil;
- delete p_ccr;
+ delete cil;
+ delete ccr;
delete[] file_lines;
printf("Done.\n");
return 0;