From 04a47a77330b1fb43d1df83652243f4d87bba7d5 Mon Sep 17 00:00:00 2001
From: "Mulas, Giacomo" <gmulas@oa-cagliari.inaf.it>
Date: Fri, 12 Jan 2024 18:21:27 +0100
Subject: [PATCH] * Update Makefiles and make.inc to make use of wildcards for
standard compiles * rename trapping/frfme.cpp to trapping/cfrfme.cpp * rename
trapping/lffft.cpp to trapping/clffft.cpp * this makes them compliant with
make wildcards
---
src/cluster/Makefile | 32 +--
src/make.inc | 5 +
src/sphere/Makefile | 28 +--
src/trapping/Makefile | 36 ++--
src/trapping/cfrfme.cpp | 452 ++++++++++++++++++++++++++++++++++++++++
src/trapping/clffft.cpp | 391 ++++++++++++++++++++++++++++++++++
6 files changed, 896 insertions(+), 48 deletions(-)
create mode 100644 src/trapping/cfrfme.cpp
create mode 100644 src/trapping/clffft.cpp
diff --git a/src/cluster/Makefile b/src/cluster/Makefile
index d5138014..37c27125 100644
--- a/src/cluster/Makefile
+++ b/src/cluster/Makefile
@@ -14,29 +14,29 @@ edfb: edfb.o
np_cluster: $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/clu_subs.o $(BUILDDIR)/cluster.o
$(CXX) $(CXXFLAGS) -o $(BUILDDIR)/np_cluster $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/clu_subs.o $(BUILDDIR)/cluster.o $(CXXLDFLAGS)
-$(BUILDDIR)/np_cluster.o:
- $(CXX) $(CXXFLAGS) -c np_cluster.cpp -o $(BUILDDIR)/np_cluster.o
+#$(BUILDDIR)/np_cluster.o:
+# $(CXX) $(CXXFLAGS) -c np_cluster.cpp -o $(BUILDDIR)/np_cluster.o
-$(BUILDDIR)/Commons.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/Commons.cpp -o $(BUILDDIR)/Commons.o
+#$(BUILDDIR)/Commons.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/Commons.cpp -o $(BUILDDIR)/Commons.o
-$(BUILDDIR)/Configuration.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/Configuration.cpp -o $(BUILDDIR)/Configuration.o
+#$(BUILDDIR)/Configuration.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/Configuration.cpp -o $(BUILDDIR)/Configuration.o
-$(BUILDDIR)/file_io.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/file_io.cpp -o $(BUILDDIR)/file_io.o
+#$(BUILDDIR)/file_io.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/file_io.cpp -o $(BUILDDIR)/file_io.o
-$(BUILDDIR)/Parsers.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o
+#$(BUILDDIR)/Parsers.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o
-$(BUILDDIR)/cluster.o:
- $(CXX) $(CXXFLAGS) -c cluster.cpp -o $(BUILDDIR)/cluster.o
+#$(BUILDDIR)/cluster.o:
+# $(CXX) $(CXXFLAGS) -c cluster.cpp -o $(BUILDDIR)/cluster.o
-$(BUILDDIR)/clu_subs.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/clu_subs.cpp -o $(BUILDDIR)/clu_subs.o
+#$(BUILDDIR)/clu_subs.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/clu_subs.cpp -o $(BUILDDIR)/clu_subs.o
-$(BUILDDIR)/sph_subs.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/sph_subs.cpp -o $(BUILDDIR)/sph_subs.o
+#$(BUILDDIR)/sph_subs.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/sph_subs.cpp -o $(BUILDDIR)/sph_subs.o
clean:
rm -f $(BUILDDIR)/*.o
diff --git a/src/make.inc b/src/make.inc
index eafa41ac..ac95b212 100644
--- a/src/make.inc
+++ b/src/make.inc
@@ -35,3 +35,8 @@ endif
%.o : %.cpp
$(CXX) $(CXXFLAGS) -c -o $(BUILDDIR)/$@ $<
+$(BUILDDIR)/%.o : %.cpp
+ $(CXX) $(CXXFLAGS) -c -o $(BUILDDIR)/$@ $<
+
+$(BUILDDIR)/%.o : ../libnptm/%.cpp
+ $(CXX) $(CXXFLAGS) -c -o $(BUILDDIR)/$@ ../libnptm/$<
diff --git a/src/sphere/Makefile b/src/sphere/Makefile
index 79c2aae6..384a85b9 100644
--- a/src/sphere/Makefile
+++ b/src/sphere/Makefile
@@ -20,26 +20,26 @@ sph: sph.o
np_sphere: $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/sphere.o
$(CXX) $(CXXFLAGS) -o $(BUILDDIR)/np_sphere $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/sphere.o $(CXXLDFLAGS)
-$(BUILDDIR)/np_sphere.o:
- $(CXX) $(CXXFLAGS) -c np_sphere.cpp -o $(BUILDDIR)/np_sphere.o
+#$(BUILDDIR)/np_sphere.o:
+# $(CXX) $(CXXFLAGS) -c np_sphere.cpp -o $(BUILDDIR)/np_sphere.o
-$(BUILDDIR)/Commons.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/Commons.cpp -o $(BUILDDIR)/Commons.o
+#$(BUILDDIR)/Commons.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/Commons.cpp -o $(BUILDDIR)/Commons.o
-$(BUILDDIR)/Configuration.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/Configuration.cpp -o $(BUILDDIR)/Configuration.o
+#$(BUILDDIR)/Configuration.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/Configuration.cpp -o $(BUILDDIR)/Configuration.o
-$(BUILDDIR)/file_io.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/file_io.cpp -o $(BUILDDIR)/file_io.o
+#$(BUILDDIR)/file_io.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/file_io.cpp -o $(BUILDDIR)/file_io.o
-$(BUILDDIR)/Parsers.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o
+#$(BUILDDIR)/Parsers.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o
-$(BUILDDIR)/sph_subs.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/sph_subs.cpp -o $(BUILDDIR)/sph_subs.o
+#$(BUILDDIR)/sph_subs.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/sph_subs.cpp -o $(BUILDDIR)/sph_subs.o
-$(BUILDDIR)/sphere.o:
- $(CXX) $(CXXFLAGS) -c sphere.cpp -o $(BUILDDIR)/sphere.o
+#$(BUILDDIR)/sphere.o:
+# $(CXX) $(CXXFLAGS) -c sphere.cpp -o $(BUILDDIR)/sphere.o
clean:
rm -f $(BUILDDIR)/*.o
diff --git a/src/trapping/Makefile b/src/trapping/Makefile
index 9c9ff9d4..10d5bb8b 100644
--- a/src/trapping/Makefile
+++ b/src/trapping/Makefile
@@ -20,32 +20,32 @@ lffft: lffft.o
np_trapping: $(BUILDDIR)/np_trapping.o $(BUILDDIR)/cfrfme.o $(BUILDDIR)/clffft.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/file_io.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/tra_subs.o
$(CXX) $(CXXFLAGS) -o $(BUILDDIR)/np_trapping $(BUILDDIR)/np_trapping.o $(BUILDDIR)/cfrfme.o $(BUILDDIR)/clffft.o $(BUILDDIR)/file_io.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/tra_subs.o $(CXXLDFLAGS)
-$(BUILDDIR)/np_trapping.o:
- $(CXX) $(CXXFLAGS) np_trapping.cpp -c -o $(BUILDDIR)/np_trapping.o
+#$(BUILDDIR)/np_trapping.o:
+# $(CXX) $(CXXFLAGS) np_trapping.cpp -c -o $(BUILDDIR)/np_trapping.o
-$(BUILDDIR)/cfrfme.o:
- $(CXX) $(CXXFLAGS) frfme.cpp -c -o $(BUILDDIR)/cfrfme.o
+#$(BUILDDIR)/cfrfme.o:
+# $(CXX) $(CXXFLAGS) frfme.cpp -c -o $(BUILDDIR)/cfrfme.o
-$(BUILDDIR)/clffft.o:
- $(CXX) $(CXXFLAGS) lffft.cpp -c -o $(BUILDDIR)/clffft.o
+#$(BUILDDIR)/clffft.o:
+# $(CXX) $(CXXFLAGS) lffft.cpp -c -o $(BUILDDIR)/clffft.o
-$(BUILDDIR)/Commons.o:
- $(CXX) $(CXXFLAGS) ../libnptm/Commons.cpp -c -o $(BUILDDIR)/Commons.o
+#$(BUILDDIR)/Commons.o:
+# $(CXX) $(CXXFLAGS) ../libnptm/Commons.cpp -c -o $(BUILDDIR)/Commons.o
-$(BUILDDIR)/Configuration.o:
- $(CXX) $(CXXFLAGS) ../libnptm/Configuration.cpp -c -o $(BUILDDIR)/Configuration.o
+#$(BUILDDIR)/Configuration.o:
+# $(CXX) $(CXXFLAGS) ../libnptm/Configuration.cpp -c -o $(BUILDDIR)/Configuration.o
-$(BUILDDIR)/file_io.o:
- $(CXX) $(CXXFLAGS) -c ../libnptm/file_io.cpp -o $(BUILDDIR)/file_io.o
+#$(BUILDDIR)/file_io.o:
+# $(CXX) $(CXXFLAGS) -c ../libnptm/file_io.cpp -o $(BUILDDIR)/file_io.o
-$(BUILDDIR)/Parsers.o:
- $(CXX) $(CXXFLAGS) ../libnptm/Parsers.cpp -c -o $(BUILDDIR)/Parsers.o
+#$(BUILDDIR)/Parsers.o:
+# $(CXX) $(CXXFLAGS) ../libnptm/Parsers.cpp -c -o $(BUILDDIR)/Parsers.o
-$(BUILDDIR)/sph_subs.o:
- $(CXX) $(CXXFLAGS) ../libnptm/sph_subs.cpp -c -o $(BUILDDIR)/sph_subs.o
+#$(BUILDDIR)/sph_subs.o:
+# $(CXX) $(CXXFLAGS) ../libnptm/sph_subs.cpp -c -o $(BUILDDIR)/sph_subs.o
- $(BUILDDIR)/tra_subs.o:
- $(CXX) $(CXXFLAGS) ../libnptm/tra_subs.cpp -c -o $(BUILDDIR)/tra_subs.o
+#$(BUILDDIR)/tra_subs.o:
+# $(CXX) $(CXXFLAGS) ../libnptm/tra_subs.cpp -c -o $(BUILDDIR)/tra_subs.o
clean:
rm -f $(BUILDDIR)/*.o
diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp
new file mode 100644
index 00000000..251ea911
--- /dev/null
+++ b/src/trapping/cfrfme.cpp
@@ -0,0 +1,452 @@
+/*! \file frfme.cpp
+ */
+#include <complex>
+#include <cstdio>
+#include <fstream>
+#include <regex>
+#include <string>
+
+#ifndef INCLUDE_PARSERS_H_
+#include "../include/Parsers.h"
+#endif
+
+#ifndef INCLUDE_COMMONS_H_
+#include "../include/Commons.h"
+#endif
+
+#ifndef INCLUDE_SPH_SUBS_H_
+#include "../include/sph_subs.h"
+#endif
+
+#ifndef INCLUDE_TRA_SUBS_H_
+#include "../include/tra_subs.h"
+#endif
+
+using namespace std;
+
+/*! \brief C++ implementation of FRFME
+ *
+ * \param data_file: `string` Name of the input data file.
+ * \param output_path: `string` Directory to write the output files in.
+ */
+void frfme(string data_file, string output_path) {
+ string tfrfme_name = output_path + "/c_TFRFME";
+ fstream tfrfme;
+ char namef[7];
+ char more;
+ double *xv = NULL, *yv = NULL, *zv = NULL;
+ double *vkv = NULL, **vkzm = NULL;
+ complex<double> *wk = NULL, **w = NULL, **wsum = NULL;
+ const complex<double> cc0(0.0, 0.0);
+ const complex<double> uim(0.0, 1.0);
+ int line_count = 0, last_read_line = 0;
+ regex re = regex("-?[0-9]+");
+ string *file_lines = load_file(data_file, &line_count);
+ smatch m;
+ string str_target = file_lines[last_read_line++];
+ regex_search(str_target, m, re);
+ int jlmf = stoi(m.str());
+ 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;
+ // Vector size variables
+ int vkzm_size, wsum_size;
+ // End of vector size variables
+ if (jlmf != 1) {
+ int nxv, nyv, nzv;
+ tfrfme.open(tfrfme_name, ios::in | ios::binary);
+ if (tfrfme.is_open()) {
+ tfrfme.read(reinterpret_cast<char *>(&lmode), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&lm), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&nkv), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&nxv), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&nyv), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&nzv), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&vk), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&exri), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&an), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&ff), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&tra), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&spd), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&frsh), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&exril), sizeof(double));
+ xv = new double[nxv]();
+ yv = new double[nyv]();
+ zv = new double[nzv]();
+ for (int xi = 0; xi < nxv; xi++) tfrfme.read(reinterpret_cast<char *>(&(xv[xi])), sizeof(double));
+ 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";
+ 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;
+ for (int vki = 0; vki < nkv; vki++) vkzm[vki] = new double[nkv]();
+ for (int jy10 = 0; jy10 < nkv; jy10++) {
+ for (int jx10 = 0; jx10 < nkv; jx10++) {
+ temptape2.read(reinterpret_cast<char *>(&(vkzm[jx10][jy10])), sizeof(double));
+ } //jx10 loop
+ } // jy10 loop
+ temptape2.read(reinterpret_cast<char *>(&apfafa), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&pmf), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&spd), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&rir), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&ftcn), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&fshmx), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&vxyzmx), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&delxyz), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&vknmx), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&delk), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&delks), sizeof(double));
+ temptape2.read(reinterpret_cast<char *>(&nlmmt), sizeof(int));
+ temptape2.read(reinterpret_cast<char *>(&nrvc), sizeof(int));
+ temptape2.close();
+ } else {
+ printf("ERROR: could not open TEMPTAPE2 file.\n");
+ }
+ for (int ixyz12 = 0; ixyz12 < nrvc; ixyz12++) {
+ for (int j12 = 0; j12 < jlmf - 1; j12++) {
+ double vreal, vimag;
+ tfrfme.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ wsum[j12][ixyz12] = complex<double>(vreal, vimag);
+ } // j12 loop
+ } // ixyz12 loop
+ tfrfme.close();
+ } else {
+ printf("ERROR: could not open TFRFME file.\n");
+ }
+ nks = nkv - 1;
+ } else { // label 16
+ int nksh, nrsh, nxsh, nysh, nzsh;
+ str_target = file_lines[last_read_line++];
+ for (int cli = 0; cli < 7; cli++) {
+ regex_search(str_target, m, re);
+ if (cli == 0) lmode = stoi(m.str());
+ else if (cli == 1) lm = stoi(m.str());
+ else if (cli == 2) nksh = stoi(m.str());
+ else if (cli == 3) nrsh = stoi(m.str());
+ else if (cli == 4) nxsh = stoi(m.str());
+ else if (cli == 5) nysh = stoi(m.str());
+ else if (cli == 6) nzsh = stoi(m.str());
+ str_target = m.suffix().str();
+ }
+ re = regex("-?[0-9]\\.[0-9]+([dDeE][-+]?[0-9]+)?");
+ regex_search(str_target, m, re);
+ double wlenfr = stod(m.str());
+ str_target = file_lines[last_read_line++];
+ for (int cli = 0; cli < 3; cli++) {
+ regex_search(str_target, m, re);
+ if (cli == 0) an = stod(m.str());
+ else if (cli == 1) ff = stod(m.str());
+ else if (cli == 2) tra = stod(m.str());
+ str_target = m.suffix().str();
+ }
+ double spdfr, exdcl;
+ str_target = file_lines[last_read_line++];
+ for (int cli = 0; cli < 3; cli++) {
+ regex_search(str_target, m, re);
+ if (cli == 0) spd = stod(m.str());
+ else if (cli == 1) spdfr = stod(m.str());
+ else if (cli == 2) exdcl = stod(m.str());
+ str_target = m.suffix().str();
+ }
+ str_target = file_lines[last_read_line++];
+ re = regex("[eEmM]");
+ if (regex_search(str_target, m, re)) {
+ more = m.str().at(0);
+ if (more == 'm' || more == 'M') {
+ more = 'M';
+ sprintf(namef, "c_TMDF");
+ }
+ else if (more == 'e' || more == 'E') {
+ more = 'E';
+ sprintf(namef, "c_TEDF");
+ }
+ str_target = m.suffix().str();
+ re = regex("[0-9]+");
+ regex_search(str_target, m, re);
+ int ixi = stoi(m.str());
+ fstream tedf;
+ string tedf_name = output_path + "/" + namef;
+ tedf.open(tedf_name.c_str(), ios::in | ios::binary);
+ if (tedf.is_open()) {
+ int iduml, idum;
+ tedf.read(reinterpret_cast<char *>(&iduml), sizeof(int));
+ for (int i = 0; i < iduml; i++) tedf.read(reinterpret_cast<char *>(&idum), sizeof(int));
+ tedf.read(reinterpret_cast<char *>(&exdc), sizeof(double));
+ tedf.read(reinterpret_cast<char *>(&wp), sizeof(double));
+ tedf.read(reinterpret_cast<char *>(&xip), sizeof(double));
+ tedf.read(reinterpret_cast<char *>(&idfc), sizeof(int));
+ tedf.read(reinterpret_cast<char *>(&nxi), sizeof(int));
+ if (idfc >= 0) {
+ if (ixi <= nxi) {
+ for (int i = 0; i < ixi; i++) tedf.read(reinterpret_cast<char *>(&xi), sizeof(double));
+ } else { // label 96
+ tedf.close();
+ // label 98
+ string output_name = output_path + "/c_OFRFME";
+ FILE *output = fopen(output_name.c_str(), "w");
+ fprintf(output, " WRONG INPUT TAPE\n");
+ fclose(output);
+ }
+ } else { // label 18
+ xi = xip;
+ }
+ // label 20
+ tedf.close();
+ double wn = wp / 3.0e8;
+ vk = xi * wn;
+ exri = sqrt(exdc);
+ frsh = 0.0;
+ exril = 0.0;
+ fshmx = 0.0;
+ apfafa = exri / (an * ff);
+ if (lmode != 0) pmf = 2.0 * apfafa;
+ if (spd > 0.0) {
+ exril = sqrt(exdcl);
+ rir = exri / exril;
+ ftcn = 2.0 / (1.0 + rir);
+ frsh = -spd * spdfr;
+ double sthmx = an / exri;
+ double sthlmx = sthmx * rir;
+ double uy = 1.0;
+ fshmx = spd * (rir * (sqrt(uy - sthmx * sthmx) / sqrt(uy - sthlmx * sthlmx)) - uy);
+ }
+ // label 22
+ nlmmt = lm * (lm + 2) * 2;
+ nks = nksh * 2;
+ nkv = nks + 1;
+ double vkm = vk * exri;
+ vknmx = vk * an;
+ delk = vknmx / nksh;
+ delks = delk / vkm;
+ delks = delks * delks;
+ vxyzmx = acos(0.0) * 4.0 / vkm * wlenfr;
+ delxyz = vxyzmx / nrsh;
+ int nxs = nxsh * 2;
+ 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];
+ } // i24 loop
+ int nys = nysh * 2;
+ 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];
+ } // i25 loop
+ int nzs = nzsh * 2;
+ 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];
+ } // i28 loop
+ tfrfme.open(tfrfme_name.c_str(), ios::out | ios::binary);
+ if (tfrfme.is_open()) {
+ tfrfme.write(reinterpret_cast<char *>(&lmode), sizeof(int));
+ tfrfme.write(reinterpret_cast<char *>(&lm), sizeof(int));
+ tfrfme.write(reinterpret_cast<char *>(&nkv), sizeof(int));
+ tfrfme.write(reinterpret_cast<char *>(&nxv), sizeof(int));
+ tfrfme.write(reinterpret_cast<char *>(&nyv), sizeof(int));
+ tfrfme.write(reinterpret_cast<char *>(&nxv), sizeof(int));
+ tfrfme.write(reinterpret_cast<char *>(&vk), sizeof(double));
+ tfrfme.write(reinterpret_cast<char *>(&exri), sizeof(double));
+ tfrfme.write(reinterpret_cast<char *>(&an), sizeof(double));
+ tfrfme.write(reinterpret_cast<char *>(&ff), sizeof(double));
+ tfrfme.write(reinterpret_cast<char *>(&tra), sizeof(double));
+ tfrfme.write(reinterpret_cast<char *>(&spd), sizeof(double));
+ tfrfme.write(reinterpret_cast<char *>(&frsh), sizeof(double));
+ tfrfme.write(reinterpret_cast<char *>(&exril), sizeof(double));
+ for (int xi = 0; xi < nxv; xi++)
+ tfrfme.write(reinterpret_cast<char *>(&(xv[xi])), sizeof(double));
+ for (int yi = 0; yi < nyv; yi++)
+ tfrfme.write(reinterpret_cast<char *>(&(yv[yi])), sizeof(double));
+ for (int zi = 0; zi < nzv; zi++)
+ tfrfme.write(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));
+ fstream temptape1, temptape2;
+ string temp_name1 = output_path + "/c_TEMPTAPE1";
+ string temp_name2 = output_path + "/c_TEMPTAPE2";
+ temptape1.open(temp_name1.c_str(), ios::out | ios::binary);
+ temptape2.open(temp_name2.c_str(), ios::out | ios::binary);
+ for (int jx = 0; jx < nkv; jx++)
+ temptape2.write(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
+ frfmer(nkv, vkm, vkv, vknmx, apfafa, tra, spd, rir, ftcn, lm, lmode, pmf, temptape1, temptape2);
+ temptape1.close();
+ temptape2.write(reinterpret_cast<char *>(&apfafa), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&pmf), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&spd), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&rir), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&ftcn), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&fshmx), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&vxyzmx), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&delxyz), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&vknmx), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&delk), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&delks), sizeof(double));
+ temptape2.write(reinterpret_cast<char *>(&nlmmt), sizeof(int));
+ 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 jy40 = 0; jy40 < nkv; jy40++) {
+ for (int jx40 = 0; jx40 < nkv; jx40++)
+ temptape2.read(reinterpret_cast<char *>(&(vkzm[jx40][jy40])), sizeof(double));
+ } // 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]();
+ for (int j80 = jlmf - 1; j80 < jlml; j80++) {
+ temptape1.open(temp_name1.c_str(), ios::in | ios::binary);
+ for (int jy50 = 0; jy50 < nkv; jy50++) {
+ for (int jx50 = 0; jx50 < nkv; jx50++) {
+ for (int i = 0; i < nlmmt; i++) {
+ double vreal, vimag;
+ temptape1.read(reinterpret_cast<char *>(&vreal), sizeof(double));
+ temptape1.read(reinterpret_cast<char *>(&vimag), sizeof(double));
+ wk[i] = complex<double>(vreal, vimag);
+ }
+ w[jx50][jy50] = wk[j80];
+ } // jx50
+ } // jy50 loop
+ temptape1.close();
+ int ixyz = 0;
+ for (int iz75 = 0; iz75 < nzv; iz75++) {
+ double z = zv[iz75] + frsh;
+ for (int iy70 = 0; iy70 < nyv; iy70++) {
+ double y = yv[iy70];
+ for (int ix65 = 0; ix65 < nxv; ix65++) {
+ double x = xv[ix65];
+ ixyz++;
+ complex<double> sumy = cc0;
+ for (int jy60 = 0; jy60 < nkv; jy60++) {
+ double vky = vkv[jy60];
+ double vkx = vkv[nkv - 1];
+ double vkzf = vkzm[0][jy60];
+ complex<double> phasf = exp(uim * (-vkx * x + vky * y +vkzf * z));
+ 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];
+ complex<double> phas = exp(uim * (vkx * x + vky * y + vkz * z));
+ sumx += (w[jx55][jy60] * phas);
+ } // jx55 loop
+ if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5;
+ sumy += sumx;
+ } // jy60 loop
+ wsum[j80][ixyz - 1] = sumy * delks;
+ } // ix65 loop
+ } // iy70 loop
+ } // iz75 loop
+ } // j80 loop
+ if (jlmf != 1) {
+ tfrfme.open(tfrfme_name, ios::in | ios::out | ios::binary);
+ tfrfme.read(reinterpret_cast<char *>(&lmode), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&lm), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&nkv), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&nxv), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&nyv), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&nzv), sizeof(int));
+ tfrfme.read(reinterpret_cast<char *>(&vk), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&exri), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&an), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&ff), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&tra), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&spd), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&frsh), sizeof(double));
+ tfrfme.read(reinterpret_cast<char *>(&exril), sizeof(double));
+ for (int i = 0; i < nxv; i++) tfrfme.read(reinterpret_cast<char *>(&(xv[i])), sizeof(double));
+ for (int i = 0; i < nyv; i++) tfrfme.read(reinterpret_cast<char *>(&(yv[i])), sizeof(double));
+ for (int i = 0; i < nzv; i++) tfrfme.read(reinterpret_cast<char *>(&(zv[i])), sizeof(double));
+ }
+ // label 88
+ for (int ixyz = 0; ixyz < nrvc; ixyz++) {
+ for (int j = 0; j< jlml; j++) {
+ double vreal = wsum[j][ixyz].real();
+ double vimag = wsum[j][ixyz].imag();
+ tfrfme.write(reinterpret_cast<char *>(&vreal), sizeof(double));
+ tfrfme.write(reinterpret_cast<char *>(&vimag), sizeof(double));
+ } // j loop
+ } // ixyz loop
+ tfrfme.close();
+ string output_name = output_path + "/c_OFRFME";
+ FILE *output = fopen(output_name.c_str(), "w");
+ fprintf(output, " IF JLML < NLMMT, PRESERVE TEMPTAPE1, TEMPTAPE2, AND TFRFRME,\n");
+ fprintf(output, " AND RESTART LM RUN WITH JLMF = JLML+1\n");
+ if (spd > 0.0) fprintf(output, " FSHMX =%15.7lE\n", fshmx);
+ fprintf(output, " FRSH =%15.7lE\n", frsh);
+ fclose(output);
+ } else { // Should never happen.
+ printf("ERROR: could not open TFRFME file for output.\n");
+ }
+ } else {
+ printf("ERROR: could not open TEDF file.\n");
+ }
+ } else { // label 98
+ string output_name = output_path + "/c_OFRFME";
+ FILE *output = fopen(output_name.c_str(), "w");
+ fprintf(output, " WRONG INPUT TAPE\n");
+ fclose(output);
+ }
+ }
+ // label 45
+ if (tfrfme.is_open()) tfrfme.close();
+ delete[] file_lines;
+ if (xv != NULL) delete[] xv;
+ if (yv != NULL) delete[] yv;
+ if (zv != NULL) delete[] zv;
+ if (vkv != NULL) delete[] vkv;
+ if (vkzm != NULL) {
+ for (int vki = vkzm_size - 1; vki > -1; vki--) delete[] vkzm[vki];
+ delete[] vkzm;
+ }
+ 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");
+}
diff --git a/src/trapping/clffft.cpp b/src/trapping/clffft.cpp
new file mode 100644
index 00000000..d41fb77d
--- /dev/null
+++ b/src/trapping/clffft.cpp
@@ -0,0 +1,391 @@
+/*! \file lffft.cpp
+ */
+#include <complex>
+#include <cstdio>
+#include <fstream>
+#include <regex>
+#include <string>
+
+#ifndef INCLUDE_PARSERS_H_
+#include "../include/Parsers.h"
+#endif
+
+#ifndef INCLUDE_COMMONS_H_
+#include "../include/Commons.h"
+#endif
+
+#ifndef INCLUDE_SPH_SUBS_H_
+#include "../include/sph_subs.h"
+#endif
+
+#ifndef INCLUDE_TRA_SUBS_H_
+#include "../include/tra_subs.h"
+#endif
+
+using namespace std;
+
+/*! \brief C++ implementation of LFFFT
+ *
+ * \param data_file: `string` Name of the input data file.
+ * \param output_path: `string` Directory to write the output files in.
+ */
+void lffft(string data_file, string output_path) {
+ 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;
+ double ****zpv = NULL;
+ double *xv = NULL, *yv = NULL, *zv = NULL;
+ complex<double> *ac = NULL, *ws = NULL, *wsl = NULL;
+ complex<double> **am0m = NULL;
+ complex<double> **amd = NULL;
+ int **indam = NULL;
+ complex<double> *tmsm = NULL, *tmse = NULL, **tms = NULL;
+ int jft, jss, jtw;
+ int is, le, nvam = 0;
+ double vks, exris;
+ CIL *cil = new CIL();
+ CCR *ccr = new CCR();
+
+ int num_lines = 0;
+ string *file_lines = load_file(data_file, &num_lines);
+ regex re = regex("-?[0-9]+");
+ smatch m;
+ string str_target = file_lines[0];
+ for (int mi = 0; mi < 3; mi++) {
+ regex_search(str_target, m, re);
+ if (mi == 0) jft = stoi(m.str());
+ else if (mi == 1) jss = stoi(m.str());
+ else if (mi == 2) jtw = stoi(m.str());
+ str_target = m.suffix().str();
+ } // mi loop
+ string ttms_name = output_path + "/c_TTMS";
+ fstream ttms;
+ ttms.open(ttms_name, ios::in | ios::binary);
+ if (ttms.is_open()) {
+ ttms.read(reinterpret_cast<char *>(&is), sizeof(int));
+ ttms.read(reinterpret_cast<char *>(&le), sizeof(int));
+ ttms.read(reinterpret_cast<char *>(&vks), sizeof(double));
+ ttms.read(reinterpret_cast<char *>(&exris), sizeof(double));
+ 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]();
+ // 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>*[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));
+ 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));
+ cil->mxmpo = vint;
+ 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;
+ ccr->cof = 1.0 / sqk;
+ ccr->cimu = ccr->cof / sqrt(2.0);
+ if (jss != 1) {
+ string tlfff_name = output_path + "/c_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) {
+ 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
+ }
+ // label 180
+ if (is != 2222) {
+ if (is != 1111) {
+ if (is > 0) { // Goes to 305
+ ac = new complex<double>[cil->nlemt]();
+ camp(ac, am0m, ws, cil);
+ // Goes to 445
+ } else if (is < 0) { // Goes to 405
+ ac = new complex<double>[cil->nlemt]();
+ czamp(ac, amd, indam, ws, cil);
+ // Goes to 445
+ }
+ } else {
+ ac = new complex<double>[cil->nlemt]();
+ samp(ac, tmsm, tmse, ws, cil);
+ // Goes to 445
+ }
+ } else {
+ ac = new complex<double>[cil->nlemt]();
+ sampoa(ac, tms, ws, cil);
+ // Goes to 445
+ }
+ bool goto475 = false;
+ // label 445
+ if (jft <= 0) {
+ double *fffe = new double[3]();
+ double *fffs = new double[3]();
+ ffrf(zpv, ac, ws, fffe, fffs, cil, ccr);
+ 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]();
+ ffrt(ac, ws, ffte, ffts, cil);
+ 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;
+ }
+ } // i loop
+ } // ix475 loop
+ } // iy475 loop
+ } // iz475 loop
+ if (jss != 1) {
+ if (jft <= 0) tlfff.close();
+ if (jft >= 0) tlfft.close();
+ }
+ }
+ }
+ 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;
+ if (xv != NULL) delete[] xv;
+ if (yv != NULL) delete[] yv;
+ if (zv != NULL) delete[] zv;
+ if (wsl != NULL) delete[] wsl;
+ if (tmsm != NULL) delete[] tmsm;
+ if (tmse != NULL) delete[] tmse;
+ if (tms != NULL) {
+ for (int ti = le - 1; ti > -1; ti--) delete[] tms[ti];
+ delete[] tms;
+ }
+ if (am0m != NULL) {
+ for (int ai = 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 cil;
+ delete ccr;
+ delete[] file_lines;
+ printf("Done.\n");
+}
--
GitLab