From b966ce93353b733ef86ccb8f2b418d03a668b9eb Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Tue, 2 Jan 2024 19:03:16 +0100
Subject: [PATCH] Extend migration of LFFFT to C++
---
src/include/tra_subs.h | 1 +
src/trapping/lffft.cpp | 206 +++++++++++++++++++++++++++++++++++++++--
2 files changed, 198 insertions(+), 9 deletions(-)
diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index b839b405..780c078c 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -22,6 +22,7 @@ extern void sphar(
double cosrth, double sinrth, double cosrph, double sinrph,
int ll, std::complex<double> *ylm
);
+extern void thdps(int lm, double ****zpv);
//End of externally defined subroutines
/*! C++ porting of PWMALP
diff --git a/src/trapping/lffft.cpp b/src/trapping/lffft.cpp
index 49606385..500f4230 100644
--- a/src/trapping/lffft.cpp
+++ b/src/trapping/lffft.cpp
@@ -21,6 +21,7 @@ using namespace std;
* \param output_path: `string` Directory to write the output files in.
*/
int main() {
+ const complex<double> uim(0.0, 1.0);
string data_file = "../../test_data/trapping/DLFFFT";
string output_path = ".";
double ****zpv = NULL, *ac = NULL;
@@ -32,10 +33,12 @@ int main() {
complex<double> **am0m = NULL;
complex<double> **amd = NULL;
int **indam = NULL;
- complex<double> *tmsm = NULL, *tmse = NULL, *tms = NULL;
+ complex<double> *tmsm = NULL, *tmse = NULL, **tms = NULL;
int jft, jss, jtw;
- int is, le;
+ int is, le, nvam = 0;
double vks, exris;
+ cil *p_cil = new cil();
+ ccr *p_ccr = new ccr();
int num_lines = 0;
string *file_lines = load_file(data_file, &num_lines);
@@ -57,15 +60,172 @@ 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));
- const int nlem = le * (le + 2);
- const int nlemt = nlem + nlem;
- printf("DEBUG: is = %d\n", is);
- printf("DEBUG: le = %d\n", le);
- printf("DEBUG: vks = %lE\n", vks);
- printf("DEBUG: exris = %lE\n", exris);
+ p_cil->le = le;
+ p_cil->nlem = le * (le + 2);
+ p_cil->nlemt = p_cil->nlem + p_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]();
+ // QUESTION|WARNING: original code uses LM without defining it. Where does it come from?
+ int lm = le;
+ for (int i = 0; i < lm; i++) {
+ double vreal, vimag;
+ ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ tms[i][0] = complex<double>(vreal, vimag);
+ ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ tms[i][1] = complex<double>(vreal, vimag);
+ ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ tms[i][2] = complex<double>(vreal, vimag);
+ } // i loop
+ } else if (is > 1111) { // label 125
+ tmsm = new complex<double>[le]();
+ tmse = new complex<double>[le]();
+ for (int i = 0; i < le; i++) {
+ double vreal, vimag;
+ ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ tmsm[i] = complex<double>(vreal, vimag);
+ ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ 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++) {
+ double vreal, vimag;
+ ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ am0m[i][j] = complex<double>(vreal, vimag);
+ } // j loop
+ } // i loop
+ } else if (is < 0) {
+ nvam = le * le + (le * (le + 1) * (le * 2 + 1)) / 3;
+ amd = new complex<double>*[nvam];
+ for (int ai = 0; ai < nvam; ai++) amd[ai] = new complex<double>[4]();
+ for (int i = 0; i < nvam; i++) {
+ for (int j = 0; j < 4; j++) {
+ double vreal, vimag;
+ ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ amd[i][j] = complex<double>(vreal, vimag);
+ } // j loop
+ } // i loop
+ indam = new int*[le];
+ int vint;
+ for (int ii = 0; ii < le; ii++) indam[ii] = new int[le]();
+ for (int i = 0; i < le; i++) {
+ for (int j = 0; j < le; j++) {
+ ttms.read(reinterpret_cast<char *>(&vint), sizeof(int));
+ indam[i][j] = vint;
+ } // j loop
+ } // i loop
+ ttms.read(reinterpret_cast<char *>(&vint), sizeof(int));
+ p_cil->mxmpo = vint;
+ p_cil->mxim = vint * 2 - 1;
+ }
+ // label 150
+ ttms.close();
+ fstream binary_input;
+ string binary_name;
+ if (jss != 1) binary_name = output_path + "/c_TFRFME";
+ else binary_name = output_path + "c_TWS";
+ binary_input.open(binary_name, ios::in | ios::binary);
+ if (binary_input.is_open()) {
+ int lmode, lm, nkv, nxv, nyv, nzv;
+ double vk, exri, an, ff, tra;
+ double spd, frsh, exril;
+ binary_input.read(reinterpret_cast<char *>(&lmode), sizeof(int));
+ binary_input.read(reinterpret_cast<char *>(&lm), sizeof(int));
+ binary_input.read(reinterpret_cast<char *>(&nkv), sizeof(int));
+ binary_input.read(reinterpret_cast<char *>(&nxv), sizeof(int));
+ binary_input.read(reinterpret_cast<char *>(&nyv), sizeof(int));
+ binary_input.read(reinterpret_cast<char *>(&nzv), sizeof(int));
+ if (lm >= le) {
+ binary_input.read(reinterpret_cast<char *>(&vk), sizeof(double));
+ binary_input.read(reinterpret_cast<char *>(&exri), sizeof(double));
+ binary_input.read(reinterpret_cast<char *>(&an), sizeof(double));
+ binary_input.read(reinterpret_cast<char *>(&ff), sizeof(double));
+ binary_input.read(reinterpret_cast<char *>(&tra), sizeof(double));
+ if (vk == vks && exri == exris) {
+ binary_input.read(reinterpret_cast<char *>(&spd), sizeof(double));
+ binary_input.read(reinterpret_cast<char *>(&frsh), sizeof(double));
+ binary_input.read(reinterpret_cast<char *>(&exril), sizeof(double));
+ xv = new double[nxv];
+ for (int i = 0; i < nxv; i++) binary_input.read(reinterpret_cast<char *>(&(xv[i])), sizeof(double));
+ yv = new double[nyv];
+ for (int i = 0; i < nyv; i++) binary_input.read(reinterpret_cast<char *>(&(yv[i])), sizeof(double));
+ zv = new double[nzv];
+ for (int i = 0; i < nzv; i++) binary_input.read(reinterpret_cast<char *>(&(zv[i])), sizeof(double));
+ bool goto160 = false;
+ if (jft <= 0) {
+ zpv = new double***[le];
+ for (int zi = 0; zi < le; zi++) {
+ zpv[zi] = new double**[3];
+ for (int zj = 0; zj < 3; zj++) {
+ zpv[zi][zj] = new double*[2];
+ for (int zk = 0; zk < 2; zk++) zpv[zi][zj][zk] = new double[2]();
+ } // zj loop
+ } // zi loop
+ 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);
+ 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));
+ tlfff.write(reinterpret_cast<char *>(&le), sizeof(int));
+ tlfff.write(reinterpret_cast<char *>(&nkv), sizeof(int));
+ tlfff.write(reinterpret_cast<char *>(&nxv), sizeof(int));
+ tlfff.write(reinterpret_cast<char *>(&nyv), sizeof(int));
+ tlfff.write(reinterpret_cast<char *>(&nzv), sizeof(int));
+ tlfff.write(reinterpret_cast<char *>(&vk), sizeof(double));
+ tlfff.write(reinterpret_cast<char *>(&exri), sizeof(double));
+ tlfff.write(reinterpret_cast<char *>(&an), sizeof(double));
+ tlfff.write(reinterpret_cast<char *>(&ff), sizeof(double));
+ tlfff.write(reinterpret_cast<char *>(&tra), sizeof(double));
+ tlfff.write(reinterpret_cast<char *>(&spd), sizeof(double));
+ tlfff.write(reinterpret_cast<char *>(&frsh), sizeof(double));
+ tlfff.write(reinterpret_cast<char *>(&exril), sizeof(double));
+ for (int i = 0; i < nxv; i++)
+ tlfff.write(reinterpret_cast<char *>(&(xv[i])), sizeof(double));
+ for (int i = 0; i < nyv; i++)
+ tlfff.write(reinterpret_cast<char *>(&(yv[i])), sizeof(double));
+ for (int i = 0; i < nzv; i++)
+ tlfff.write(reinterpret_cast<char *>(&(zv[i])), sizeof(double));
+ if (jft < 0) goto160 = true;
+ } else { // Should never happen.
+ printf("ERROR: could not open TLFFF file.\n");
+ }
+ }
+ }
+ // 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.
+ }
+ }
+ // label 160
+ }
+ }
+ binary_input.close();
+ } else {
+ printf("ERROR: could not open binary input file %s.\n", binary_name.c_str());
+ }
} else {
printf("ERROR: could not open TTMS file.\n");
}
+
// Clean up memory
if (ac != NULL) delete[] ac;
if (ws != NULL) delete[] ws;
@@ -79,7 +239,35 @@ int main() {
if (wsl != NULL) delete[] wsl;
if (tmsm != NULL) delete[] tmsm;
if (tmse != NULL) delete[] tmse;
- if (tms != NULL) delete[] tms;
+ if (tms != NULL) {
+ for (int ti = le - 1; ti > -1; ti--) delete[] tms[ti];
+ delete[] tms;
+ }
+ if (am0m != NULL) {
+ for (int ai = p_cil->nlemt - 1; ai > -1; ai--) delete[] am0m[ai];
+ delete[] am0m;
+ }
+ if (amd != NULL) {
+ for (int ai = nvam - 1; ai > -1; ai--) delete[] amd[ai];
+ delete[] amd;
+ }
+ if (indam != NULL) {
+ for (int ii = le - 1; ii > -1; ii--) delete[] indam[ii];
+ delete[] indam;
+ }
+ if (zpv != NULL) {
+ for (int zi = le - 1; zi > -1; zi--) {
+ for (int zj = 2; zj > -1; zj--) {
+ for (int zk = 1; zk > -1; zk--) delete[] zpv[zi][zj][zk];
+ delete[] zpv[zi][zj];
+ } // zj loop
+ delete[] zpv[zi];
+ } // zi loop
+ delete[] zpv;
+ }
+ delete p_cil;
+ delete p_ccr;
+ delete[] file_lines;
printf("Done.\n");
return 0;
}
--
GitLab