diff --git a/src/trapping/frfme.cpp b/src/trapping/frfme.cpp
deleted file mode 100644
index 251ea9117b5100bfc70b6bfccde7e49ef787a1e2..0000000000000000000000000000000000000000
--- a/src/trapping/frfme.cpp
+++ /dev/null
@@ -1,452 +0,0 @@
-/*! \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/lffft.cpp b/src/trapping/lffft.cpp
deleted file mode 100644
index d41fb77dba055f42b5f2658c33baba50c612d78f..0000000000000000000000000000000000000000
--- a/src/trapping/lffft.cpp
+++ /dev/null
@@ -1,391 +0,0 @@
-/*! \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");
-}