/* Copyright 2024 INAF - Osservatorio Astronomico di Cagliari

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
 */

/*! \file clffft.cpp
 *
 * \brief C++ implementation of LFFFT functions.
 */
#include <cstdio>
#include <exception>
#include <fstream>
#include <regex>
#include <string>

#ifndef INCLUDE_TYPES_H_
#include "../include/types.h"
#endif

#ifndef INCLUDE_PARSERS_H_
#include "../include/Parsers.h"
#endif

#ifndef INCLUDE_LOGGING_H_
#include "../include/logging.h"
#endif

#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif

#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.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 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 dcomplex uim = 0.0 + 1.0 * I;
  const double sq2i = 1.0 / sqrt(2.0);
  const dcomplex sq2iti = sq2i * uim;

  fstream tlfff, tlfft;
  double ****zpv = NULL;
  dcomplex *ac = NULL, *ws = NULL, *wsl = NULL;
  dcomplex **am0m = NULL;
  dcomplex **amd = NULL;
  int **indam = NULL;
  dcomplex *tmsm = NULL, *tmse = NULL, **tms = NULL;
  dcomplex **_wsum = 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 dcomplex*[le];
      for (int ti = 0; ti < le; ti++) tms[ti] = new dcomplex[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] = vreal + vimag * I;
	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
	tms[i][1] = vreal + vimag * I;
	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
	tms[i][2] = vreal + vimag * I;
      } // i loop
    } else if (is >= 1111) { // label 125
      tmsm = new dcomplex[le]();
      tmse = new dcomplex[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] = vreal + vimag * I;
	ttms.read(reinterpret_cast<char *>(&vreal), sizeof(double));
	ttms.read(reinterpret_cast<char *>(&vimag), sizeof(double));
	tmse[i] = vreal + vimag * I;
      } // i loop
    } else if (is >= 0) { // label 135
      am0m = new dcomplex*[cil->nlemt];
      for (int ai = 0; ai < cil->nlemt; ai++) am0m[ai] = new dcomplex[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] = vreal + vimag * I;
	} // j loop
      } // i loop
    } else if (is < 0) {
      nvam = le * le + (le * (le + 1) * (le * 2 + 1)) / 3;
      amd = new dcomplex*[nvam];
      for (int ai = 0; ai < nvam; ai++) amd[ai] = new dcomplex[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] = vreal + vimag * I;
	} // 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();
    TFRFME *tfrfme = NULL;
    string binary_name;
    if (jss != 1) binary_name = output_path + "/c_TFRFME.hd5";
    else binary_name = output_path + "/c_TWS";
    tfrfme = TFRFME::from_binary(binary_name, "HDF5");
    if (tfrfme != NULL) {
      int lmode, lm, nkv, nxv, nyv, nzv;
      double vk, exri, an, ff, tra;
      double spd, frsh, exril;
      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");
      _wsum = tfrfme->get_matrix();
      double *_xv = tfrfme->get_x();
      double *_yv = tfrfme->get_y();
      double *_zv = tfrfme->get_z();
      if (lm >= le) {
	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");
	if (vk == vks && exri == exris) {
	  spd = tfrfme->get_param("spd");
	  frsh = tfrfme->get_param("frsh");
	  exril = tfrfme->get_param("exril");
	  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++) {
		  double x = _xv[i];
		  tlfff.write(reinterpret_cast<char *>(&x), sizeof(double));
		}
		for (int i = 0; i < nyv; i++) {
		  double y = _yv[i];
		  tlfff.write(reinterpret_cast<char *>(&y), sizeof(double));
		}
		for (int i = 0; i < nzv; i++) {
		  double z = _zv[i];
		  tlfff.write(reinterpret_cast<char *>(&z), 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) {
	      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++) {
		  double x = _xv[i];
		  tlfft.write(reinterpret_cast<char *>(&x), sizeof(double));
		}
		for (int i = 0; i < nyv; i++) {
		  double y = _yv[i];
		  tlfft.write(reinterpret_cast<char *>(&y), sizeof(double));
		}
		for (int i = 0; i < nzv; i++) {
		  double z = _zv[i];
		  tlfft.write(reinterpret_cast<char *>(&z), 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 dcomplex[nlmmt]();
	  if (lm > le) wsl = new dcomplex[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));
		  int row = i;
		  int col = (nyv * nxv * iz475) + (nxv * iy475) + ix475;
		  dcomplex value = _wsum[row][col];
		  if (lm <= le) {
		    ws[i] = value;
		  } else { // label 170
		    wsl[i] = value;
		    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
		  }
		} // i loop
		// label 180
		if (is != 2222) {
		  if (is != 1111) {
		    if (is > 0) { // Goes to 305
		      ac = new dcomplex[cil->nlemt]();
		      camp(ac, am0m, ws, cil);
		      // Goes to 445
		    } else if (is < 0) { // Goes to 405
		      ac = new dcomplex[cil->nlemt]();
		      czamp(ac, amd, indam, ws, cil);
		      // Goes to 445
		    }
		  } else {
		    ac = new dcomplex[cil->nlemt]();
		    samp(ac, tmsm, tmse, ws, cil);
		    // Goes to 445
		  }
		} else {
		  ac = new dcomplex[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;
		}
	      } // ix475 loop
	    } // iy475 loop
	  } // iz475 loop
	  if (jss != 1) {
	    if (jft <= 0) tlfff.close();
	    if (jft >= 0) tlfft.close();
	  }
	  fclose(output66);
	  fclose(output67);
	}
      }
      delete tfrfme;
    } 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 (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("LFFT: Done.\n");
}
