/*! \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 = 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;
  int nlmmt = 0, nrvc = 0;
  // Vector size variables
  int 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));
      vkv = new double[nkv]();
      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()) {
	for (int jx = 0; jx < nkv; jx++) temptape2.read(reinterpret_cast<char *>(&(vkv[jx])), sizeof(double));
	vkzm = new double*[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;
	// Array initialization
	vkv = new double[nkv]();
	vkzm = new double*[nkv];
	for (int vi = 0; vi < nkv; vi++) vkzm[vi] = new double[nkv];
	// 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;
	xv = new double[nxv]();
	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]();
	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]();
	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
	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);
	  for (int jx = 0; jx < nkv; jx++) {
	    double value = 0.0;
	    temptape2.read(reinterpret_cast<char *>(&value), sizeof(double));
	    vkv[jx] = value;
	  }
	  for (int jy40 = 0; jy40 < nkv; jy40++) {
	    for (int jx40 = 0; jx40 < nkv; jx40++) {
	      double value = 0.0;
	      temptape2.read(reinterpret_cast<char *>(&value), sizeof(double));
	      vkzm[jx40][jy40] = value;
	    }
	  } // jy40 loop
	  temptape2.close();
	  if (wsum != NULL) {
	    for (int wsi = wsum_size - 1; wsi > -1; wsi--) delete[] wsum[wsi];
	    delete[] wsum;
	  }
	  wsum = new complex<double>*[nlmmt];
	  for (int j80 = jlmf - 1; j80 < jlml; j80++) {
	    // 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]();
	    if (wk != NULL) delete[] wk;
	    wk = new complex<double>[nlmmt]();
	    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;
	    wsum[j80] = new complex<double>[nrvc]();
	    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][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 = nkv - 1; vki > -1; vki--) delete[] vkzm[vki];
    delete[] vkzm;
  }
  if (w != NULL) {
    for (int wi = nkv - 1; wi > -1; wi--) delete[] w[wi];
    delete[] w;
  }
  if (wsum != NULL) {
    for (int wsi = wsum_size - 1; wsi > -1; wsi--) delete[] wsum[wsi];
    delete[] wsum;
  }
  if (wk != NULL) delete[] wk;
  printf("Done.\n");
}
