/*! \file cfrfme.cpp
*
* C++ implementation of FRFME functions.
*/
#include <complex>
#include <cstdio>
#include <exception>
#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_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif
#ifndef INCLUDE_SPH_SUBS_H_
#include "../include/sph_subs.h"
#endif
#ifndef INCLUDE_TFRFME_H_
#include "../include/tfrfme.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.hd5";
TFRFME *tfrfme = NULL;
Swap1 *tt1 = NULL;
Swap2 *tt2 = NULL;
char namef[7];
char more;
complex<double> **w = NULL;
complex<double> *wk = 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 = 0, lm = 0, nks = 0, nkv = 0;
double vk = 0.0, exri = 0.0, an = 0.0, ff = 0.0, tra = 0.0;
double exdc = 0.0, wp = 0.0, xip = 0.0, xi = 0.0;
int idfc = 0, nxi = 0;
double apfafa = 0.0, pmf = 0.0, spd = 0.0, rir = 0.0, ftcn = 0.0, fshmx = 0.0;
double vxyzmx = 0.0, delxyz = 0.0, vknmx = 0.0, delk = 0.0, delks = 0.0;
double frsh = 0.0, exril = 0.0;
double **vkzm = NULL;
double *vkv = NULL;
int nlmmt = 0, nrvc = 0;
// Vector size variables
int wsum_size;
// End of vector size variables
if (jlmf != 1) {
int nxv, nyv, nzv;
if (tfrfme == NULL) tfrfme = TFRFME::from_binary(tfrfme_name, "HDF5");
if (tfrfme != NULL) {
lmode = (int)tfrfme->get_param("lmode");
lm = (int)tfrfme->get_param("lm");
nkv = (int)tfrfme->get_param("nkv");
nxv = (int)tfrfme->get_param("nxv");
nyv = (int)tfrfme->get_param("nyv");
nzv = (int)tfrfme->get_param("nzv");
vk = tfrfme->get_param("vk");
exri = tfrfme->get_param("exri");
an = tfrfme->get_param("an");
ff = tfrfme->get_param("ff");
tra = tfrfme->get_param("tra");
spd = tfrfme->get_param("spd");
frsh = tfrfme->get_param("frsh");
exril = tfrfme->get_param("exril");
string tempname2 = output_path + "/c_TEMPTAPE2.hd5";
if (tt2 == NULL) tt2 = Swap2::from_binary(tempname2, "HDF5");
if (tt2 != NULL) {
vkv = tt2->get_vector();
vkzm = tt2->get_matrix();
apfafa = tt2->get_param("apfafa");
pmf = tt2->get_param("pmf");
spd = tt2->get_param("spd");
rir = tt2->get_param("rir");
ftcn = tt2->get_param("ftcn");
fshmx = tt2->get_param("fshmx");
vxyzmx = tt2->get_param("vxyzmx");
delxyz = tt2->get_param("delxyz");
vknmx = tt2->get_param("vknmx");
delk = tt2->get_param("delk");
delks = tt2->get_param("delks");
nlmmt = (int)tt2->get_param("nlmmt");
nrvc = (int)tt2->get_param("nrvc");
} else {
printf("ERROR: could not open TEMPTAPE2 file.\n");
}
} 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());
string tedf_name = output_path + "/" + namef + ".hd5";
ScattererConfiguration *tedf = ScattererConfiguration::from_binary(tedf_name, "HDF5");
if (tedf != NULL) {
int iduml, idum;
iduml = (int)tedf->get_param("nsph");
idum = tedf->get_iog(iduml - 1);
exdc = tedf->get_param("exdc");
wp = tedf->get_param("wp");
xip = tedf->get_param("xip");
idfc = (int)tedf->get_param("idfc");
nxi = (int)tedf->get_param("nxi");
if (idfc >= 0) {
if (ixi <= nxi) {
xi = tedf->get_scale(ixi - 1);
} else { // label 96
delete tedf;
// 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
delete tedf;
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;
// Array initialization
long swap1_size, swap2_size, tfrfme_size;
double size_mb;
printf("INFO: calculation memory requirements\n");
swap1_size = Swap1::get_memory_requirement(lm, nkv);
size_mb = 1.0 * swap1_size / 1024.0 / 1024.0;
printf("Swap 1: %.2lg MB\n", size_mb);
swap2_size = Swap2::get_memory_requirement(nkv);
size_mb = 1.0 * swap2_size / 1024.0 / 1024.0;
printf("Swap 2: %.2lg MB\n", size_mb);
tt2 = new Swap2(nkv);
vkv = tt2->get_vector();
vkzm = tt2->get_matrix();
// End of array initialization
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;
int nys = nysh * 2;
int nyv = nys + 1;
int nyshpo = nysh + 1;
int nzs = nzsh * 2;
int nzv = nzs + 1;
int nzshpo = nzsh + 1;
tfrfme_size = TFRFME::get_memory_requirement(lmode, lm, nkv, nxv, nyv, nzv);
size_mb = 1.0 * tfrfme_size / 1024.0 / 1024.0;
printf("TFRFME: %.2lg MB\n", size_mb);
size_mb = 1.0 * (swap1_size + swap2_size + tfrfme_size) / 1024.0 / 1024.0;
printf("TOTAL: %.2lg MB\n", size_mb);
tfrfme = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
double *_xv = tfrfme->get_x();
double *_yv = tfrfme->get_y();
double *_zv = tfrfme->get_z();
complex<double> **_wsum = tfrfme->get_matrix();
for (int i24 = nxshpo; i24 <= nxs; i24++) {
_xv[i24] = _xv[i24 - 1] + delxyz;
_xv[nxv - i24 - 1] = -_xv[i24];
} // i24 loop
for (int i25 = nyshpo; i25 <= nys; i25++) {
_yv[i25] = _yv[i25 - 1] + delxyz;
_yv[nyv - i25 - 1] = -_yv[i25];
} // i25 loop
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;
for (int i28 = nkshpo; i28 <= nks; i28++) {
vkv[i28] = vkv[i28 - 1] + delk;
vkv[nkv - i28 - 1] = -vkv[i28];
} // i28 loop
if (tfrfme != NULL) {
tfrfme->set_param("vk", vk);
tfrfme->set_param("exri", exri);
tfrfme->set_param("an", an);
tfrfme->set_param("ff", ff);
tfrfme->set_param("tra", tra);
tfrfme->set_param("spd", spd);
tfrfme->set_param("frsh", frsh);
tfrfme->set_param("exril", exril);
//fstream temptape1, temptape2;
string temp_name1 = output_path + "/c_TEMPTAPE1.hd5";
string temp_name2 = output_path + "/c_TEMPTAPE2.hd5";
//temptape1.open(temp_name1.c_str(), ios::out | ios::binary);
tt1 = new Swap1(lm, nkv);
if (wk != NULL) delete[] wk;
wk = frfmer(nkv, vkm, vknmx, apfafa, tra, spd, rir, ftcn, lm, lmode, pmf, tt1, tt2);
tt1->write_binary(temp_name1, "HDF5");
//delete tt1;
tt2->set_param("apfafa", apfafa);
tt2->set_param("pmf", pmf);
tt2->set_param("spd", spd);
tt2->set_param("rir", rir);
tt2->set_param("ftcn", ftcn);
tt2->set_param("fshmx", fshmx);
tt2->set_param("vxyzmx", vxyzmx);
tt2->set_param("delxyz", delxyz);
tt2->set_param("vknmx", vknmx);
tt2->set_param("delk", delk);
tt2->set_param("delks", delks);
tt2->set_param("nlmmt", 1.0 * nlmmt);
tt2->set_param("nrvc", 1.0 * nrvc);
tt2->write_binary(temp_name2, "HDF5");
for (int j80 = jlmf; j80 <= jlml; j80++) {
complex<double> *tt1_wk = tt1->get_vector();
int wk_index = 0;
// w matrix
if (w != NULL) {
for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
delete[] w;
}
w = new complex<double>*[nkv];
for (int wi = 0; wi < nkv; wi++) w[wi] = new complex<double>[nkv]();
for (int jy50 = 0; jy50 < nkv; jy50++) {
for (int jx50 = 0; jx50 < nkv; jx50++) {
for (int wi = 0; wi < nlmmt; wi++) wk[wi] = tt1_wk[wk_index++];
w[jx50][jy50] = wk[j80 - 1];
} // jx50
} // jy50 loop
int ixyz = 0;
for (int wj = 0; wj < nrvc; wj++) _wsum[j80 - 1][wj] = cc0;
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 = 2; jx55 <= nks; jx55++) {
vkx = vkv[jx55 - 1];
double vkz = vkzm[jx55 - 1][jy60];
complex<double> phas = exp(uim * (vkx * x + vky * y + vkz * z));
sumx += (w[jx55 - 1][jy60] * phas);
} // jx55 loop
if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5;
sumy += sumx;
} // jy60 loop
_wsum[j80 - 1][ixyz - 1] = sumy * delks;
} // ix65 loop
} // iy70 loop
} // iz75 loop
} // j80 loop
// label 88
tfrfme->write_binary(tfrfme_name, "HDF5");
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 != NULL) delete tfrfme;
delete[] file_lines;
if (tt2 != NULL) delete tt2;
if (w != NULL) {
for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
delete[] w;
}
if (wk != NULL) delete[] wk;
if (tt1 != NULL) delete tt1;
printf("FRFME: Done.\n");
}