From 9a2b8adb2ee9eabba69ecb2759f7789dafe7fc10 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Wed, 21 Feb 2024 17:36:07 +0100
Subject: [PATCH] Divert TFRFME I/O to RAM

---
 src/trapping/cfrfme.cpp | 212 ++++++++++++++++------------------------
 src/trapping/clffft.cpp | 104 +++++++++++---------
 2 files changed, 141 insertions(+), 175 deletions(-)

diff --git a/src/trapping/cfrfme.cpp b/src/trapping/cfrfme.cpp
index 22c46110..29d21e8f 100644
--- a/src/trapping/cfrfme.cpp
+++ b/src/trapping/cfrfme.cpp
@@ -4,6 +4,7 @@
  */
 #include <complex>
 #include <cstdio>
+#include <exception>
 #include <fstream>
 #include <regex>
 #include <string>
@@ -16,10 +17,18 @@
 #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
@@ -32,13 +41,13 @@ using namespace std;
  *  \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;
+  string tfrfme_name = output_path + "/c_TFRFME.hd5";
+  //fstream tfrfme;
+  TFRFME *tfrfme = NULL;
   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;
+  complex<double> *wk = NULL, **w = 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;
@@ -64,29 +73,23 @@ void frfme(string data_file, string output_path) {
   // 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));
+    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("nkv");
+      nyv = (int)tfrfme->get_param("nkv");
+      nzv = (int)tfrfme->get_param("nkv");
+      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");
       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);
@@ -116,15 +119,6 @@ void frfme(string data_file, string output_path) {
       } 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");
     }
@@ -179,23 +173,22 @@ void frfme(string data_file, string output_path) {
       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()) {
+      string tedf_name = output_path + "/" + namef + ".hd5";
+      ScattererConfiguration *tedf = ScattererConfiguration::from_binary(tedf_name, "HDF5");
+      if (tedf != NULL) {
 	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));
+	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) {
-	    for (int i = 0; i < ixi; i++) tedf.read(reinterpret_cast<char *>(&xi), sizeof(double));
+	    xi = tedf->get_scale(ixi - 1);
 	  } else { // label 96
-	    tedf.close();
+	    delete tedf;
 	    // label 98
 	    string output_name = output_path + "/c_OFRFME";
 	    FILE *output = fopen(output_name.c_str(), "w");
@@ -206,7 +199,7 @@ void frfme(string data_file, string output_path) {
 	  xi = xip;
 	}
 	// label 20
-	tedf.close();
+	delete tedf;
 	double wn = wp / 3.0e8;
 	vk = xi * wn;
 	exri = sqrt(exdc);
@@ -244,26 +237,24 @@ void frfme(string data_file, string output_path) {
 	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]();
+	tfrfme = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
+	for (int i24 = nxshpo; i24 <= nxs; i24++) {
+	  tfrfme->set_x(i24, tfrfme->get_x(i24 - 1) + delxyz);
+	  tfrfme->set_x(nxv - i24 - 1, -tfrfme->get_x(i24));
+	} // i24 loop
+	for (int i25 = nyshpo; i25 <= nys; i25++) {
+	  tfrfme->set_y(i25, tfrfme->get_y(i25 - 1) + delxyz);
+	  tfrfme->set_y(nyv - i25 - 1, -tfrfme->get_y(i25));
+	} // i25 loop
 	for (int i27 = nzshpo; i27 <= nzs; i27++) {
-	  zv[i27] = zv[i27 - 1] + delxyz;
-	  zv[nzv - i27 - 1] = -zv[i27];
+	  tfrfme->set_z(i27, tfrfme->get_z(i27 - 1) + delxyz);
+	  tfrfme->set_z(nzv - i27 - 1, -tfrfme->get_z(i27));
 	} // i27 loop
 	int nrvc = nxv * nyv * nzv;
 	int nkshpo = nksh + 1;
@@ -271,28 +262,15 @@ void frfme(string data_file, string output_path) {
 	  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 *>(&nzv), 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));
+	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";
 	  string temp_name2 = output_path + "/c_TEMPTAPE2";
@@ -330,11 +308,6 @@ void frfme(string data_file, string output_path) {
 	    }
 	  } // 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) {
@@ -359,13 +332,13 @@ void frfme(string data_file, string output_path) {
 	    } // jy50 loop
 	    temptape1.close();
 	    int ixyz = 0;
-	    wsum[j80] = new complex<double>[nrvc]();
+	    for (int wj = 0; wj < nrvc; wj++) tfrfme->set_matrix_element(j80, wj, cc0);
 	    for (int iz75 = 0; iz75 < nzv; iz75++) {
-	      double z = zv[iz75]  + frsh;
+	      double z = tfrfme->get_z(iz75)  + frsh;
 	      for (int iy70 = 0; iy70 < nyv; iy70++) {
-		double y = yv[iy70];
+		double y = tfrfme->get_y(iy70);
 		for (int ix65 = 0; ix65 < nxv; ix65++) {
-		  double x = xv[ix65];
+		  double x = tfrfme->get_x(ix65);
 		  ixyz++;
 		  complex<double> sumy = cc0;
 		  for (int jy60 = 0; jy60 < nkv; jy60++) {
@@ -385,41 +358,29 @@ void frfme(string data_file, string output_path) {
 		    if (jy60 == 0 || jy60 == nkv - 1) sumx *= 0.5;
 		    sumy += sumx;
 		  } // jy60 loop
-		  wsum[j80][ixyz - 1] = sumy * delks;
+		  tfrfme->set_matrix_element(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));
+	    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");
 	  }
 	  // 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();
+	  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");
@@ -441,11 +402,8 @@ void frfme(string data_file, string output_path) {
     }
   }
   // label 45
-  if (tfrfme.is_open()) tfrfme.close();
+  if (tfrfme != NULL) delete tfrfme;
   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];
@@ -455,10 +413,6 @@ void frfme(string data_file, string output_path) {
     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");
 }
diff --git a/src/trapping/clffft.cpp b/src/trapping/clffft.cpp
index 9cd9b50b..380e7e9e 100644
--- a/src/trapping/clffft.cpp
+++ b/src/trapping/clffft.cpp
@@ -4,6 +4,7 @@
  */
 #include <complex>
 #include <cstdio>
+#include <exception>
 #include <fstream>
 #include <regex>
 #include <string>
@@ -20,6 +21,10 @@
 #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
@@ -38,7 +43,6 @@ void lffft(string data_file, string output_path) {
 
   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;
@@ -140,37 +144,31 @@ void lffft(string data_file, string output_path) {
     }
     // label 150
     ttms.close();
-    fstream binary_input;
+    TFRFME *tfrfme = NULL;
     string binary_name;
-    if (jss != 1) binary_name = output_path + "/c_TFRFME";
+    if (jss != 1) binary_name = output_path + "/c_TFRFME.hd5";
     else binary_name = output_path + "/c_TWS";
-    binary_input.open(binary_name, ios::in | ios::binary);
-    if (binary_input.is_open()) {
+    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;
-      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));
+      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");
       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));
+	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) {
-	  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));
+	  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];
@@ -204,12 +202,18 @@ void lffft(string data_file, string output_path) {
 		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));
+		for (int i = 0; i < nxv; i++) {
+		  double x = tfrfme->get_x(i);
+		  tlfff.write(reinterpret_cast<char *>(&x), sizeof(double));
+		}
+		for (int i = 0; i < nyv; i++) {
+		  double y = tfrfme->get_y(i);
+		  tlfff.write(reinterpret_cast<char *>(&y), sizeof(double));
+		}
+		for (int i = 0; i < nzv; i++) {
+		  double z = tfrfme->get_z(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");
@@ -236,12 +240,18 @@ void lffft(string data_file, string output_path) {
 		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));
+		for (int i = 0; i < nxv; i++) {
+		  double x = tfrfme->get_x(i);
+		  tlfft.write(reinterpret_cast<char *>(&x), sizeof(double));
+		}
+		for (int i = 0; i < nyv; i++) {
+		  double y = tfrfme->get_y(i);
+		  tlfft.write(reinterpret_cast<char *>(&y), sizeof(double));
+		}
+		for (int i = 0; i < nzv; i++) {
+		  double z = tfrfme->get_z(i);
+		  tlfft.write(reinterpret_cast<char *>(&z), sizeof(double));
+		}
 	      } else { // Should never happen.
 		printf("ERROR: could not open TLFFT file.\n");
 	      }
@@ -262,13 +272,18 @@ void lffft(string data_file, string output_path) {
 	    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));
+		  //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 = (nzv * iz475) + (nyv * iy475) + ix475;
+		  complex<double> value = tfrfme->get_matrix_element(row, col);
 		  if (lm <= le) {
-		    ws[i] = complex<double>(vreal, vimag);
+		    //ws[i] = complex<double>(vreal, vimag);
+		    ws[i] = value;
 		  } else { // label 170
-		    wsl[i] = complex<double>(vreal, vimag);
+		    //wsl[i] = complex<double>(vreal, vimag);
+		    wsl[i] = value;
 		    for (int i175 = 0; i175 < cil->nlem; i175++) {
 		      int ie = i175 + cil->nlem;
 		      int iel = i175 + nlmm;
@@ -384,7 +399,7 @@ void lffft(string data_file, string output_path) {
 	  fclose(output67);
 	}
       }
-      binary_input.close();
+      delete tfrfme;
     } else {
       printf("ERROR: could not open binary input file %s.\n", binary_name.c_str());
     }
@@ -395,9 +410,6 @@ void lffft(string data_file, string output_path) {
   // 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;
-- 
GitLab