/*! \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");
}