diff --git a/src/trapping/frfme.cpp b/src/trapping/frfme.cpp index b67788179cdabb863de650592afe47505449aea9..b34cb37ba6147deb187b2035fc155ba01cd62a2a 100644 --- a/src/trapping/frfme.cpp +++ b/src/trapping/frfme.cpp @@ -1,13 +1,14 @@ #include <cstdio> #include <fstream> +#include <regex> #include <string> #include <complex> #ifndef INCLUDE_PARSERS_H_ #include "../include/Parsers.h" #endif -#ifndef INCLUDE_SPH_SUBS_H_ -#include "../include/sph_subs.h" -#endif +//#ifndef INCLUDE_SPH_SUBS_H_ +//#include "../include/sph_subs.h" +//#endif using namespace std; @@ -19,17 +20,324 @@ using namespace std; //void frfme(string data_file, string output_path) { int main() { string data_file = "../../test_data/trapping/DFRFME"; - char namef[5]; - chare more; - double *xv, *yv, *zv; - double *vkv, **vkzm; - complex<double> *wk, **w, **wsum; + string output_path = "."; + 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); int line_count = 0, last_read_line = 0; + regex re = regex("-?[0-9]+"); string *file_lines = load_file(data_file, &line_count); - for (int fli = 0; fli < line_count; fli++) { - printf("%s\n", file_lines[fli]); + 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; + fstream tfrfme; + tfrfme.open("c_TFRFME", 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; + temptape2.open("c_TEMPTAPE2", 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"); + } + wsum = new complex<double>*[nlmmt]; + wsum_size = nlmmt; + for (int wsi = 0; wsi < jlmf - 1; wsi++) wsum[wsi] = new complex<double>[nrvc](); + 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; + tedf.open(namef, 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; + //wsum = new complex<double>*[nlmmt]; + wsum_size = nlmmt; + 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; + 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 + fstream tfrfme; + string tfrfme_name = output_path + "/c_TFRFME"; + 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)); + // Would call FRFMER(NKV,VKM,VKV,VKNMX,APFAFA,TRA,SPD,RIR,FTCN,LM,LMODE,PMF,ITT1,ITT2) + temptape1.write(reinterpret_cast<char *>(&apfafa), sizeof(double)); // Place-holder + 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(); + tfrfme.close(); + } 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 + 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; + } return 0; }