Something went wrong on our end
-
Giovanni La Mura authoredGiovanni La Mura authored
clffft.cpp 15.13 KiB
/*! \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]();
// FORTRAN writes two output formatted files without opening them
// explicitly. It is assumed thay can be opened here.
string out66_name = output_path + "/c_out66.txt";
string out67_name = output_path + "/c_out67.txt";
FILE *output66 = fopen(out66_name.c_str(), "w");
FILE *output67 = fopen(out67_name.c_str(), "w");
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
fprintf(
output66, " %18.16lE%18.16lE%18.16lE\n",
fffe[0], fffs[0], fffe[0] - fffs[0]
);
fprintf(
output66, " %18.16lE%18.16lE%18.16lE\n",
fffe[1], fffs[1], fffe[1] - fffs[1]
);
fprintf(
output66, " %18.16lE%18.16lE%18.16lE\n",
fffe[2], fffs[2], fffe[2] - fffs[2]
);
} 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
fprintf(
output66, " %5d%4d%4d%15.4lE%15.4lE%15.4lE\n",
ix475 + 1, iy475 + 1, iz475 + 1,
fffe[0] - fffs[0], fffe[1] - fffs[1], fffe[2] - fffs[2]
);
}
}
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
fprintf(
output67, " %18.16lE%18.16lE%18.16lE\n",
ffte[0], ffts[0], ffte[0] - ffts[0]
);
fprintf(
output67, " %18.16lE%18.16lE%18.16lE\n", ffte[1], ffts[1], ffte[1] - ffts[1]
);
fprintf(
output67, " %18.16lE%18.16lE%18.16lE\n",
ffte[2], ffts[2], ffte[2] - ffts[2]
);
} 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
fprintf(
output67, " %5d%4d%4d%15.4lE%15.4lE%15.4lE\n",
ix475 + 1, iy475 + 1, iz475 + 1,
ffte[0] - ffts[0], ffte[1] - ffts[1], ffte[2] - ffts[2]
);
}
}
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();
}
fclose(output66);
fclose(output67);
}
}
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");
}