/* Distributed under the terms of GPLv3 or later. See COPYING for details. */

/*! \file tfrfme.cpp
 *
 * \brief Implementation of the trapping calculation objects.
 */

#include <exception>
#include <fstream>
#include <hdf5.h>
#include <string>

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

#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif

#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
#endif

#ifndef INCLUDE_TFRFME_H_
#include "../include/tfrfme.h"
#endif

#ifndef INCLUDE_FILE_IO_H_
#include "../include/file_io.h"
#endif

using namespace std;

// >>> START OF Swap1 CLASS IMPLEMENTATION <<<
Swap1::Swap1(int lm, int _nkv) {
  nkv = _nkv;
  nlmmt = 2 * lm * (lm + 2);
  const int size = nkv * nkv * nlmmt;
  wk = new dcomplex[size]();
  last_index = 0;
}

Swap1* Swap1::from_binary(string file_name, string mode) {
  Swap1 *instance = NULL;
  if (mode.compare("LEGACY") == 0) {
    instance = from_legacy(file_name);
  } else if (mode.compare("HDF5") == 0) {
    instance = from_hdf5(file_name);
  } else {
    string message = "Unknown format mode: \"" + mode + "\"";
    throw UnrecognizedFormatException(message);
  }
  return instance;
}

Swap1* Swap1::from_hdf5(string file_name) {
  Swap1 *instance = NULL;
  unsigned int flags = H5F_ACC_RDONLY;
  HDFFile *hdf_file = new HDFFile(file_name, flags);
  herr_t status = hdf_file->get_status();
  double *elements;
  string str_type;
  int _nlmmt, _nkv, lm, num_elements, index;
  dcomplex value;
  dcomplex *_wk = NULL;
  if (status == 0) {
    status = hdf_file->read("NLMMT", "INT32", &_nlmmt);
    status = hdf_file->read("NKV", "INT32", &_nkv);
    lm = (int)((-2.0 + sqrt(4.0 + 2.0 * _nlmmt)) / 2.0);
    num_elements = 2 * _nlmmt * _nkv * _nkv;
    instance = new Swap1(lm, _nkv);
    _wk = instance->get_vector();
    elements = new double[num_elements]();
    str_type = "FLOAT64_(" + to_string(num_elements) + ")";
    status = hdf_file->read("WK", str_type, elements);
    for (int wi = 0; wi < num_elements / 2; wi++) {
      index = 2 * wi;
      value = elements[index] + elements[index + 1] * I;
      _wk[wi] = value;
    } // wi loop
    delete[] elements;
    status = hdf_file->close();
    delete hdf_file;
  }
  return instance;
}

Swap1* Swap1::from_legacy(string file_name) {
  fstream input;
  Swap1 *instance = NULL;
  dcomplex *_wk = NULL;
  int _nlmmt, _nkv, lm;
  double rval, ival;
  input.open(file_name.c_str(), ios::in | ios::binary);
  if (input.is_open()) {
    input.read(reinterpret_cast<char *>(&_nlmmt), sizeof(int));
    lm = (int)((-2.0 + sqrt(4.0 + 2.0 * _nlmmt)) / 2.0);
    input.read(reinterpret_cast<char *>(&_nkv), sizeof(int));
    instance = new Swap1(lm, _nkv);
    _wk = instance->get_vector();
    int num_elements = _nlmmt * _nkv * _nkv;
    for (int j = 0; j < num_elements; j++) {
      input.read(reinterpret_cast<char *>(&rval), sizeof(double));
      input.read(reinterpret_cast<char *>(&ival), sizeof(double));
      _wk[j] = rval + ival * I;
    }
    input.close();
  } else {
    printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
  }
  return instance;
}

long Swap1::get_memory_requirement(int lm, int _nkv) {
  long size = (long)(3 * sizeof(int));
  size += (long)(sizeof(dcomplex) * 2 * lm * (lm + 2) * _nkv * _nkv);
  return size;
}

void Swap1::write_binary(string file_name, string mode) {
  if (mode.compare("LEGACY") == 0) {
    write_legacy(file_name);
  } else if (mode.compare("HDF5") == 0) {
    write_hdf5(file_name);
  } else {
    string message = "Unknown format mode: \"" + mode + "\"";
    throw UnrecognizedFormatException(message);
  }
}

void Swap1::write_hdf5(string file_name) {
  List<string> rec_name_list(1);
  List<string> rec_type_list(1);
  List<void *> rec_ptr_list(1);
  herr_t status;
  string str_type;
  int num_elements = 2 * nlmmt * nkv * nkv;
  rec_name_list.set(0, "NLMMT");
  rec_type_list.set(0, "INT32_(1)");
  rec_ptr_list.set(0, &nlmmt);
  rec_name_list.append("NKV");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&nkv);
  rec_name_list.append("WK");
  str_type = "FLOAT64_(" + to_string(num_elements) + ")";
  rec_type_list.append(str_type);
  double *ptr_elements = new double[num_elements]();
  for (int wi = 0; wi < num_elements / 2; wi++) {
    ptr_elements[2 * wi] = real(wk[wi]);
    ptr_elements[2 * wi + 1] = imag(wk[wi]);
  }
  rec_ptr_list.append(ptr_elements);

  string *rec_names = rec_name_list.to_array();
  string *rec_types = rec_type_list.to_array();
  void **rec_pointers = rec_ptr_list.to_array();
  const int rec_num = rec_name_list.length();
  FileSchema schema(rec_num, rec_types, rec_names);
  HDFFile *hdf_file = HDFFile::from_schema(schema, file_name, H5F_ACC_TRUNC);
  for (int ri = 0; ri < rec_num; ri++)
    status = hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
  status = hdf_file->close();

  delete[] ptr_elements;
  delete[] rec_names;
  delete[] rec_types;
  delete[] rec_pointers;
  delete hdf_file;
}

void Swap1::write_legacy(string file_name) {
  fstream output;
  double rval, ival;
  output.open(file_name.c_str(), ios::out | ios::binary);
  if (output.is_open()) {
    int num_elements = nlmmt * nkv * nkv;
    output.write(reinterpret_cast<char *>(&nlmmt), sizeof(int));
    output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
    for (int j = 0; j < num_elements; j++) {
      rval = real(wk[j]);
      ival = imag(wk[j]);
      output.write(reinterpret_cast<char *>(&rval), sizeof(double));
      output.write(reinterpret_cast<char *>(&ival), sizeof(double));
    }
    output.close();
  } else { // Should never occur
    printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
  }
}

bool Swap1::operator ==(Swap1 &other) {
  if (nlmmt != other.nlmmt) {
    return false;
  }
  if (nkv != other.nkv) {
    return false;
  }
  const int num_elements = nlmmt * nkv * nkv;
  for (int i = 0; i < num_elements; i++) {
    if (wk[i] != other.wk[i]) {
      return false;
    }
  }
  return true;
}
// >>> END OF Swap1 CLASS IMPLEMENTATION <<<

// >>> START OF Swap2 CLASS IMPLEMENTATION <<<
Swap2::Swap2(int _nkv) {
  nkv = _nkv;
  vkv = new double[nkv]();
  vkzm = new double*[nkv];
  for (int vi = 0; vi < nkv; vi++) vkzm[vi] = new double[nkv]();
  last_vector = 0;
  last_matrix = 0;
}

Swap2::~Swap2() {
  delete[] vkv;
  for (int vi = nkv - 1; vi > -1; vi--) delete[] vkzm[vi];
  delete[] vkzm;
}

Swap2* Swap2::from_binary(string file_name, string mode) {
  Swap2 *instance = NULL;
  if (mode.compare("LEGACY") == 0) {
    instance = from_legacy(file_name);
  } else if (mode.compare("HDF5") == 0) {
    instance = from_hdf5(file_name);
  } else {
    string message = "Unknown format mode: \"" + mode + "\"";
    throw UnrecognizedFormatException(message);
  }
  return instance;
}

Swap2* Swap2::from_hdf5(string file_name) {
  Swap2 *instance = NULL;
  unsigned int flags = H5F_ACC_RDONLY;
  HDFFile *hdf_file = new HDFFile(file_name, flags);
  herr_t status = hdf_file->get_status();
  string str_type;
  int _nkv, _nlmmt, _nrvc;
  double value;
  if (status == 0) {
    status = hdf_file->read("NKV", "INT32", &_nkv);
    instance = new Swap2(_nkv);
    str_type = "FLOAT64_(" + to_string(_nkv) + ")";
    status = hdf_file->read("VKV", str_type, instance->vkv);
    str_type = "FLOAT64_(" + to_string(_nkv) + "," + to_string(_nkv) + ")";
    status = hdf_file->read("VKZM", str_type, instance->vkzm);
    status = hdf_file->read("APFAFA", "FLOAT64", &value);
    instance->set_param("apfafa", value);
    status = hdf_file->read("PMF", "FLOAT64", &value);
    instance->set_param("pmf", value);
    status = hdf_file->read("SPD", "FLOAT64", &value);
    instance->set_param("spd", value);
    status = hdf_file->read("RIR", "FLOAT64", &value);
    instance->set_param("rir", value);
    status = hdf_file->read("FTCN", "FLOAT64", &value);
    instance->set_param("ftcn", value);
    status = hdf_file->read("FSHMX", "FLOAT64", &value);
    instance->set_param("fshmx", value);
    status = hdf_file->read("VXYZMX", "FLOAT64", &value);
    instance->set_param("vxyzmx", value);
    status = hdf_file->read("DELXYZ", "FLOAT64", &value);
    instance->set_param("delxyz", value);
    status = hdf_file->read("VKNMX", "FLOAT64", &value);
    instance->set_param("vknmx", value);
    status = hdf_file->read("DELK", "FLOAT64", &value);
    instance->set_param("delk", value);
    status = hdf_file->read("DELKS", "FLOAT64", &value);
    instance->set_param("delks", value);
    status = hdf_file->read("NLMMT", "INT32", &_nlmmt);
    instance->set_param("nlmmt", 1.0 * _nlmmt);
    status = hdf_file->read("NRVC", "INT32", &_nrvc);
    instance->set_param("nrvc", 1.0 * _nrvc);
    status = hdf_file->close();
    delete hdf_file;
  }
  return instance;
}

Swap2* Swap2::from_legacy(string file_name) {
  fstream input;
  Swap2 *instance = NULL;
  int _nkv, _nlmmt, _nrvc;
  double **_vkzm = NULL;
  double *_vkv = NULL;
  double value;
  input.open(file_name.c_str(), ios::in | ios::binary);
  if (input.is_open()) {
    input.read(reinterpret_cast<char *>(&_nkv), sizeof(int));
    instance = new Swap2(_nkv);
    _vkzm = instance->get_matrix();
    _vkv = instance->get_vector();
    for (int vj = 0; vj < _nkv; vj++) {
      input.read(reinterpret_cast<char *>(&value), sizeof(double));
      _vkv[vj] = value;
    }
    for (int mi = 0; mi < _nkv; mi++) {
      for (int mj = 0; mj < _nkv; mj++) {
	input.read(reinterpret_cast<char *>(&value), sizeof(double));
	_vkzm[mi][mj] = value;
      }
    }
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("apfafa", value);
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("pmf", value);
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("spd", value);
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("rir", value);
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("ftcn", value);
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("fshmx", value);
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("vxyzmx", value);
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("delxyz", value);
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("vknmx", value);
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("delk", value);
    input.read(reinterpret_cast<char *>(&value), sizeof(double));
    instance->set_param("delks", value);
    input.read(reinterpret_cast<char *>(&_nlmmt), sizeof(int));
    instance->set_param("nlmmt", 1.0 * _nlmmt);
    input.read(reinterpret_cast<char *>(&_nrvc), sizeof(int));
    instance->set_param("nrvc", 1.0 * _nrvc);
    input.close();
  } else {
    printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
  }
  return instance;
}

long Swap2::get_memory_requirement(int _nkv) {
  long size = (long)(5 * sizeof(int) + 11 * sizeof(double));
  size += (long)(sizeof(double) * _nkv * (_nkv + 1));
  return size;
}

double Swap2::get_param(string param_name) {
  double value;
  if (param_name.compare("nkv") == 0) value = 1.0 * nkv;
  else if (param_name.compare("apfafa") == 0) value = apfafa;
  else if (param_name.compare("pmf") == 0) value = pmf;
  else if (param_name.compare("spd") == 0) value = spd;
  else if (param_name.compare("rir") == 0) value = rir;
  else if (param_name.compare("ftcn") == 0) value = ftcn;
  else if (param_name.compare("fshmx") == 0) value = fshmx;
  else if (param_name.compare("vxyzmx") == 0) value = vxyzmx;
  else if (param_name.compare("delxyz") == 0) value = delxyz;
  else if (param_name.compare("vknmx") == 0) value = vknmx;
  else if (param_name.compare("delk") == 0) value = delk;
  else if (param_name.compare("delks") == 0) value = delks;
  else if (param_name.compare("nlmmt") == 0) value = 1.0 * nlmmt;
  else if (param_name.compare("nrvc") == 0) value = 1.0 * nrvc;
  else {
    string message = "Unrecognized parameter name \"" + param_name + "\"";
    throw UnrecognizedParameterException(message);
  }
  return value;
}

void Swap2::push_matrix(double value) {
  int col = last_matrix % (nkv - 1);
  int row = last_matrix - (nkv * row);
  vkzm[row][col] = value;
  last_matrix++;
}

void Swap2::set_param(string param_name, double value) {
  if (param_name.compare("nkv") == 0) nkv = (int)value;
  else if (param_name.compare("apfafa") == 0) apfafa = value;
  else if (param_name.compare("pmf") == 0) pmf = value;
  else if (param_name.compare("spd") == 0) spd = value;
  else if (param_name.compare("rir") == 0) rir = value;
  else if (param_name.compare("ftcn") == 0) ftcn = value;
  else if (param_name.compare("fshmx") == 0) fshmx = value;
  else if (param_name.compare("vxyzmx") == 0) vxyzmx = value;
  else if (param_name.compare("delxyz") == 0) delxyz = value;
  else if (param_name.compare("vknmx") == 0) vknmx = value;
  else if (param_name.compare("delk") == 0) delk = value;
  else if (param_name.compare("delks") == 0) delks = value;
  else if (param_name.compare("nlmmt") == 0) nlmmt = (int)value;
  else if (param_name.compare("nrvc") == 0) nrvc = (int)value;
  else {
    string message = "Unrecognized parameter name \"" + param_name + "\"";
    throw UnrecognizedParameterException(message);
  }
}

void Swap2::write_binary(string file_name, string mode) {
  if (mode.compare("LEGACY") == 0) {
    write_legacy(file_name);
  } else if (mode.compare("HDF5") == 0) {
    write_hdf5(file_name);
  } else {
    string message = "Unknown format mode: \"" + mode + "\"";
    throw UnrecognizedFormatException(message);
  }
}

void Swap2::write_hdf5(string file_name) {
  List<string> rec_name_list(1);
  List<string> rec_type_list(1);
  List<void *> rec_ptr_list(1);
  herr_t status;
  string str_type;
  rec_name_list.set(0, "NKV");
  rec_type_list.set(0, "INT32_(1)");
  rec_ptr_list.set(0, &nkv);
  rec_name_list.append("VKV");
  str_type = "FLOAT64_(" + to_string(nkv) + ")";
  rec_type_list.append(str_type);
  rec_ptr_list.append(vkv);
  rec_name_list.append("VKZM");
  str_type = "FLOAT64_(" + to_string(nkv) + "," + to_string(nkv) + ")";
  rec_type_list.append(str_type);
  rec_ptr_list.append(vkzm);
  rec_name_list.append("APFAFA");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&apfafa);
  rec_name_list.append("PMF");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&pmf);
  rec_name_list.append("SPD");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&spd);
  rec_name_list.append("RIR");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&rir);
  rec_name_list.append("FTCN");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&ftcn);
  rec_name_list.append("FSHMX");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&fshmx);
  rec_name_list.append("VXYZMX");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&vxyzmx);
  rec_name_list.append("delxyz");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&delxyz);
  rec_name_list.append("VKNMX");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&vknmx);
  rec_name_list.append("DELK");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&delk);
  rec_name_list.append("DELKS");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&delks);
  rec_name_list.append("NLMMT");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&nlmmt);
  rec_name_list.append("NRVC");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&nrvc);

  string *rec_names = rec_name_list.to_array();
  string *rec_types = rec_type_list.to_array();
  void **rec_pointers = rec_ptr_list.to_array();
  const int rec_num = rec_name_list.length();
  FileSchema schema(rec_num, rec_types, rec_names);
  HDFFile *hdf_file = HDFFile::from_schema(schema, file_name, H5F_ACC_TRUNC);
  for (int ri = 0; ri < rec_num; ri++)
    status = hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
  status = hdf_file->close();

  delete[] rec_names;
  delete[] rec_types;
  delete[] rec_pointers;
  delete hdf_file;
}

void Swap2::write_legacy(string file_name) {
  fstream output;
  double value;
  output.open(file_name.c_str(), ios::out | ios::binary);
  if (output.is_open()) {
    output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
    for (int j = 0; j < nkv; j++){
      value = vkv[j];
      output.write(reinterpret_cast<const char*>(&value), sizeof(double));
    }
    for (int mi = 0; mi < nkv; mi++) {
      for (int mj = 0; mj < nkv; mj++) {
	value = vkzm[mi][mj];
	output.write(reinterpret_cast<const char*>(&value), sizeof(double));
      }
    }
    output.write(reinterpret_cast<const char*>(&apfafa), sizeof(double));
    output.write(reinterpret_cast<const char*>(&pmf), sizeof(double));
    output.write(reinterpret_cast<const char*>(&spd), sizeof(double));
    output.write(reinterpret_cast<const char*>(&rir), sizeof(double));
    output.write(reinterpret_cast<const char*>(&ftcn), sizeof(double));
    output.write(reinterpret_cast<const char*>(&fshmx), sizeof(double));
    output.write(reinterpret_cast<const char*>(&vxyzmx), sizeof(double));
    output.write(reinterpret_cast<const char*>(&delxyz), sizeof(double));
    output.write(reinterpret_cast<const char*>(&vknmx), sizeof(double));
    output.write(reinterpret_cast<const char*>(&delk), sizeof(double));
    output.write(reinterpret_cast<const char*>(&delks), sizeof(double));
    output.write(reinterpret_cast<const char*>(&nlmmt), sizeof(int));
    output.write(reinterpret_cast<const char*>(&nrvc), sizeof(int));
    output.close();
  } else { // Should never occur
    printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
  }
}

bool Swap2::operator ==(Swap2 &other) {
  if (nlmmt != other.nlmmt) {
    return false;
  }
  if (nrvc != other.nrvc) {
    return false;
  }
  if (nkv != other.nkv) {
    return false;
  }
  if (apfafa != other.apfafa) {
    return false;
  }
  if (pmf != other.pmf) {
    return false;
  }
  if (spd != other.spd) {
    return false;
  }
  if (rir != other.rir) {
    return false;
  }
  if (ftcn != other.ftcn) {
    return false;
  }
  if (fshmx != other.fshmx) {
    return false;
  }
  if (vxyzmx != other.vxyzmx) {
    return false;
  }
  if (delxyz != other.delxyz) {
    return false;
  }
  if (vknmx != other.vknmx) {
    return false;
  }
  if (delk != other.delk) {
    return false;
  }
  if (delks != other.delks) {
    return false;
  }
  for (int vi = 0; vi < nkv; vi++) {
    if (vkv[vi] != other.vkv[vi]) {
      return false;
    }
  }
  for (int mi = 0; mi < nkv; mi++) {
    for (int mj = 0; mj < nkv; mj++) {
      if (vkzm[mi][mj] != other.vkzm[mi][mj]) {
	return false;
      }
    }
  }
  return true;
}
// >>> END OF Swap2 CLASS IMPLEMENTATION <<<

// >>> START OF TFRFME CLASS IMPLEMENTATION <<<
TFRFME::TFRFME(int _lmode, int _lm, int _nkv, int _nxv, int _nyv, int _nzv) {
  lmode = _lmode;
  lm = _lm;
  nkv = _nkv;
  nxv = _nxv;
  nyv = _nyv;
  nzv = _nzv;
  vk = 0.0;
  exri = 0.0;
  an = 0.0;
  ff = 0.0;
  tra = 0.0;
  spd = 0.0;
  frsh = 0.0;
  exril = 0.0;

  // Array initialization
  xv = new double[nxv]();
  yv = new double[nyv]();
  zv = new double[nzv]();
  nlmmt = lm * (lm + 2) * 2;
  nrvc = nxv * nyv * nzv;
  wsum = new dcomplex*[nlmmt];
  for (int wi = 0; wi < nlmmt; wi++) wsum[wi] = new dcomplex[nrvc]();
}

TFRFME::~TFRFME() {
  delete[] xv;
  delete[] yv;
  delete[] zv;
  for (int wi = nlmmt - 1; wi > -1; wi--) delete[] wsum[wi];
  delete[] wsum;
}

TFRFME* TFRFME::from_binary(string file_name, string mode) {
  TFRFME *instance = NULL;
  if (mode.compare("LEGACY") == 0) {
    instance = from_legacy(file_name);
  } else if (mode.compare("HDF5") == 0) {
    instance = from_hdf5(file_name);
  } else {
    string message = "Unknown format mode: \"" + mode + "\"";
    throw UnrecognizedFormatException(message);
  }
  return instance;
}

TFRFME* TFRFME::from_hdf5(string file_name) {
  TFRFME *instance = NULL;
  unsigned int flags = H5F_ACC_RDONLY;
  HDFFile *hdf_file = new HDFFile(file_name, flags);
  herr_t status = hdf_file->get_status();
  double *elements;
  string str_type;
  int _nlmmt, _nrvc, num_elements;
  dcomplex value;
  dcomplex **_wsum = NULL;
  if (status == 0) {
    int lmode, lm, nkv, nxv, nyv, nzv;
    double vk, exri, an, ff, tra, spd, frsh, exril;
    status = hdf_file->read("LMODE", "INT32", &lmode);
    status = hdf_file->read("LM", "INT32", &lm);
    status = hdf_file->read("NKV", "INT32", &nkv);
    status = hdf_file->read("NXV", "INT32", &nxv);
    status = hdf_file->read("NYV", "INT32", &nyv);
    status = hdf_file->read("NZV", "INT32", &nzv);
    status = hdf_file->read("VK", "FLOAT64", &vk);
    status = hdf_file->read("EXRI", "FLOAT64", &exri);
    status = hdf_file->read("AN", "FLOAT64", &an);
    status = hdf_file->read("FF", "FLOAT64", &ff);
    status = hdf_file->read("TRA", "FLOAT64", &tra);
    status = hdf_file->read("SPD", "FLOAT64", &spd);
    status = hdf_file->read("FRSH", "FLOAT64", &frsh);
    status = hdf_file->read("EXRIL", "FLOAT64", &exril);
    instance = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
    _wsum = instance->get_matrix();
    instance->set_param("vk", vk);
    instance->set_param("exri", exri);
    instance->set_param("an", an);
    instance->set_param("ff", ff);
    instance->set_param("tra", tra);
    instance->set_param("spd", spd);
    instance->set_param("frsh", frsh);
    instance->set_param("exril", exril);
    str_type = "FLOAT64_(" + to_string(nxv) + ")";
    status = hdf_file->read("XVEC", str_type, instance->xv);
    str_type = "FLOAT64_(" + to_string(nyv) + ")";
    status = hdf_file->read("YVEC", str_type, instance->yv);
    str_type = "FLOAT64_(" + to_string(nzv) + ")";
    status = hdf_file->read("ZVEC", str_type, instance->zv);
    _nlmmt = 2 * lm * (lm + 2);
    _nrvc = nxv * nyv * nzv;
    num_elements = 2 * _nlmmt * _nrvc;
    elements = new double[num_elements]();
    str_type = "FLOAT64_(" + to_string(num_elements) + ")";
    status = hdf_file->read("WSUM", str_type, elements);
    int index = 0;
    for (int wj = 0; wj < _nrvc; wj++) {
      for (int wi = 0; wi < _nlmmt; wi++) {
	value = elements[index] + elements[index + 1] * I;
	_wsum[wi][wj] = value;
	index += 2;
      } // wi loop
    } // wj loop
    delete[] elements;
    status = hdf_file->close();
    delete hdf_file;
  }
  return instance;
}

TFRFME* TFRFME::from_legacy(string file_name) {
  fstream input;
  TFRFME *instance = NULL;
  dcomplex **_wsum = NULL;
  double *coord_vec = NULL;
  input.open(file_name.c_str(), ios::in | ios::binary);
  if (input.is_open()) {
    int lmode, lm, nkv, nxv, nyv, nzv;
    double vk, exri, an, ff, tra, spd, frsh, exril;
    double dval, rval, ival;
    input.read(reinterpret_cast<char *>(&lmode), sizeof(int));
    input.read(reinterpret_cast<char *>(&lm), sizeof(int));
    input.read(reinterpret_cast<char *>(&nkv), sizeof(int));
    input.read(reinterpret_cast<char *>(&nxv), sizeof(int));
    input.read(reinterpret_cast<char *>(&nyv), sizeof(int));
    input.read(reinterpret_cast<char *>(&nzv), sizeof(int));
    input.read(reinterpret_cast<char *>(&vk), sizeof(double));
    input.read(reinterpret_cast<char *>(&exri), sizeof(double));
    input.read(reinterpret_cast<char *>(&an), sizeof(double));
    input.read(reinterpret_cast<char *>(&ff), sizeof(double));
    input.read(reinterpret_cast<char *>(&tra), sizeof(double));
    input.read(reinterpret_cast<char *>(&spd), sizeof(double));
    input.read(reinterpret_cast<char *>(&frsh), sizeof(double));
    input.read(reinterpret_cast<char *>(&exril), sizeof(double));
    instance = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
    _wsum = instance->get_matrix();
    instance->set_param("vk", vk);
    instance->set_param("exri", exri);
    instance->set_param("an", an);
    instance->set_param("ff", ff);
    instance->set_param("tra", tra);
    instance->set_param("spd", spd);
    instance->set_param("frsh", frsh);
    instance->set_param("exril", exril);
    coord_vec = instance->get_x();
    for (int xi = 0; xi < nxv; xi++) {
      input.read(reinterpret_cast<char *>(&dval), sizeof(double));
      coord_vec[xi] = dval;
    }
    coord_vec = instance->get_y();
    for (int yi = 0; yi < nyv; yi++) {
      input.read(reinterpret_cast<char *>(&dval), sizeof(double));
      coord_vec[yi] = dval;
    }
    coord_vec = instance->get_z();
    for (int zi = 0; zi < nzv; zi++) {
      input.read(reinterpret_cast<char *>(&dval), sizeof(double));
      coord_vec[zi] = dval;
    }
    int _nlmmt = 2 * lm * (lm + 2);
    int _nrvc = nxv * nyv * nzv;
    for (int wj = 0; wj < _nrvc; wj++) {
      for (int wi = 0; wi < _nlmmt; wi++) {
	input.read(reinterpret_cast<char *>(&rval), sizeof(double));
	input.read(reinterpret_cast<char *>(&ival), sizeof(double));
	dcomplex value = rval + ival * I;
	_wsum[wi][wj] = value;
      } // wi loop
    } // wj loop
    input.close();
  } else {
    printf("ERROR: could not open input file \"%s\"\n", file_name.c_str());
  }
  return instance;
}

long TFRFME::get_memory_requirement(
				    int _lmode, int _lm, int _nkv, int _nxv,
				    int _nyv, int _nzv
) {
  int _nlmmt = 2 * _lm * (_lm + 2);
  int _nrvc = _nxv * _nyv * _nzv;
  long size = (long)(8 * sizeof(int) + 8 * sizeof(double));
  size += (long)((_nxv + _nyv + _nzv) * sizeof(double));
  size += (long)(_nlmmt * _nrvc * sizeof(dcomplex));
  return size;
}

double TFRFME::get_param(string param_name) {
  double value;
  if (param_name.compare("vk") == 0) value = vk;
  else if (param_name.compare("exri") == 0) value = exri;
  else if (param_name.compare("an") == 0) value = an;
  else if (param_name.compare("ff") == 0) value = ff;
  else if (param_name.compare("tra") == 0) value = tra;
  else if (param_name.compare("spd") == 0) value = spd;
  else if (param_name.compare("frsh") == 0) value = frsh;
  else if (param_name.compare("exril") == 0) value = exril;
  else if (param_name.compare("lmode") == 0) value = 1.0 * lmode;
  else if (param_name.compare("lm") == 0) value = 1.0 * lm;
  else if (param_name.compare("nkv") == 0) value = 1.0 * nkv;
  else if (param_name.compare("nxv") == 0) value = 1.0 * nxv;
  else if (param_name.compare("nyv") == 0) value = 1.0 * nyv;
  else if (param_name.compare("nzv") == 0) value = 1.0 * nzv;
  else {
    string message = "Unrecognized parameter name \"" + param_name + "\"";
    throw UnrecognizedParameterException(message);
  }
  return value;
}

void TFRFME::set_param(string param_name, double value) {
  if (param_name.compare("vk") == 0) vk = value;
  else if (param_name.compare("exri") == 0) exri = value;
  else if (param_name.compare("an") == 0) an = value;
  else if (param_name.compare("ff") == 0) ff = value;
  else if (param_name.compare("tra") == 0) tra = value;
  else if (param_name.compare("spd") == 0) spd = value;
  else if (param_name.compare("frsh") == 0) frsh = value;
  else if (param_name.compare("exril") == 0) exril = value;
  else {
    string message = "Unrecognized parameter name \"" + param_name + "\"";
    throw UnrecognizedParameterException(message);
  }
}

void TFRFME::write_binary(string file_name, string mode) {
  if (mode.compare("LEGACY") == 0) {
    write_legacy(file_name);
  } else if (mode.compare("HDF5") == 0) {
    write_hdf5(file_name);
  } else {
    string message = "Unknown format mode: \"" + mode + "\"";
    throw UnrecognizedFormatException(message);
  }
}

void TFRFME::write_hdf5(string file_name) {
  List<string> rec_name_list(1);
  List<string> rec_type_list(1);
  List<void *> rec_ptr_list(1);
  herr_t status;
  string str_type;
  rec_name_list.set(0, "LMODE");
  rec_type_list.set(0, "INT32_(1)");
  rec_ptr_list.set(0, &lmode);
  rec_name_list.append("LM");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&lm);
  rec_name_list.append("NKV");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&nkv);
  rec_name_list.append("NXV");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&nxv);
  rec_name_list.append("NYV");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&nyv);
  rec_name_list.append("NZV");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&nzv);
  rec_name_list.append("VK");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&vk);
  rec_name_list.append("EXRI");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&exri);
  rec_name_list.append("AN");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&an);
  rec_name_list.append("FF");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&ff);
  rec_name_list.append("TRA");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&tra);
  rec_name_list.append("SPD");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&spd);
  rec_name_list.append("FRSH");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&frsh);
  rec_name_list.append("EXRIL");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&exril);
  str_type = "FLOAT64_(" + to_string(nxv) + ")";
  rec_name_list.append("XVEC");
  rec_type_list.append(str_type);
  rec_ptr_list.append(xv);
  str_type = "FLOAT64_(" + to_string(nyv) + ")";
  rec_name_list.append("YVEC");
  rec_type_list.append(str_type);
  rec_ptr_list.append(yv);
  str_type = "FLOAT64_(" + to_string(nzv) + ")";
  rec_name_list.append("ZVEC");
  rec_type_list.append(str_type);
  rec_ptr_list.append(zv);
  rec_name_list.append("WSUM");
  int num_elements = 2 * nlmmt * nrvc;
  str_type = "FLOAT64_(" + to_string(num_elements) + ")";
  rec_type_list.append(str_type);
  // The (N x M) matrix of complex is converted to a vector of double
  // with length (2 * N * M)
  double *ptr_elements = new double[num_elements]();
  int index = 0;
  for (int wj = 0; wj < nrvc; wj++) {
    for (int wi = 0; wi < nlmmt; wi++) {
      ptr_elements[index++] = real(wsum[wi][wj]);
      ptr_elements[index++] = imag(wsum[wi][wj]);
    } // wi loop
  } // wj loop
  rec_ptr_list.append(ptr_elements);

  string *rec_names = rec_name_list.to_array();
  string *rec_types = rec_type_list.to_array();
  void **rec_pointers = rec_ptr_list.to_array();
  const int rec_num = rec_name_list.length();
  FileSchema schema(rec_num, rec_types, rec_names);
  HDFFile *hdf_file = HDFFile::from_schema(schema, file_name, H5F_ACC_TRUNC);
  for (int ri = 0; ri < rec_num; ri++)
    status = hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
  status = hdf_file->close();

  delete[] ptr_elements;
  delete[] rec_names;
  delete[] rec_types;
  delete[] rec_pointers;
  delete hdf_file;
}

void TFRFME::write_legacy(string file_name) {
  fstream output;
  output.open(file_name.c_str(), ios::out | ios::binary);
  if (output.is_open()) {
    output.write(reinterpret_cast<char *>(&lmode), sizeof(int));
    output.write(reinterpret_cast<char *>(&lm), sizeof(int));
    output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
    output.write(reinterpret_cast<char *>(&nxv), sizeof(int));
    output.write(reinterpret_cast<char *>(&nyv), sizeof(int));
    output.write(reinterpret_cast<char *>(&nzv), sizeof(int));
    output.write(reinterpret_cast<char *>(&vk), sizeof(double));
    output.write(reinterpret_cast<char *>(&exri), sizeof(double));
    output.write(reinterpret_cast<char *>(&an), sizeof(double));
    output.write(reinterpret_cast<char *>(&ff), sizeof(double));
    output.write(reinterpret_cast<char *>(&tra), sizeof(double));
    output.write(reinterpret_cast<char *>(&spd), sizeof(double));
    output.write(reinterpret_cast<char *>(&frsh), sizeof(double));
    output.write(reinterpret_cast<char *>(&exril), sizeof(double));
    for (int xi = 0; xi < nxv; xi++)
      output.write(reinterpret_cast<char *>(&(xv[xi])), sizeof(double));
    for (int yi = 0; yi < nyv; yi++)
      output.write(reinterpret_cast<char *>(&(yv[yi])), sizeof(double));
    for (int zi = 0; zi < nzv; zi++)
      output.write(reinterpret_cast<char *>(&(zv[zi])), sizeof(double));
    for (int wj = 0; wj < nrvc; wj++) {
      for (int wi = 0; wi < nlmmt; wi++) {
	double rval = real(wsum[wi][wj]);
	double ival = imag(wsum[wi][wj]);
	output.write(reinterpret_cast<char *>(&rval), sizeof(double));
	output.write(reinterpret_cast<char *>(&ival), sizeof(double));
      } // wi loop
    } // wj loop
    output.close();
  } else { // Should never occur
    printf("ERROR: could not open output file \"%s\"\n", file_name.c_str());
  }
}

bool TFRFME::operator ==(TFRFME &other) {
  if (lmode != other.lmode) {
    return false;
  }
  if (lm != other.lm) {
    return false;
  }
  if (nkv != other.nkv) {
    return false;
  }
  if (nxv != other.nxv) {
    return false;
  }
  if (nyv != other.nyv) {
    return false;
  }
  if (nzv != other.nzv) {
    return false;
  }
  if (vk != other.vk) {
    return false;
  }
  if (exri != other.exri) {
    return false;
  }
  if (an != other.an) {
    return false;
  }
  if (ff != other.ff) {
    return false;
  }
  if (tra != other.tra) {
    return false;
  }
  if (spd != other.spd) {
    return false;
  }
  if (frsh != other.frsh) {
    return false;
  }
  if (exril != other.exril) {
    return false;
  }
  for (int xi = 0; xi < nxv; xi++) {
    if (xv[xi] != other.xv[xi]) {
      return false;
    }
  }
  for (int yi = 0; yi < nyv; yi++) {
    if (yv[yi] != other.yv[yi]) {
      return false;
    }
  }
  for (int zi = 0; zi < nzv; zi++) {
    if (zv[zi] != other.zv[zi]) {
      return false;
    }
  }
  for (int wi = 0; wi < nlmmt; wi++) {
    for (int wj = 0; wj < nrvc; wj++) {
      if (wsum[wi][wj] != other.wsum[wi][wj]) {
	return false;
      }
    } // wj loop
  } // wi loop
  return true;
}
// >>> END OF TFRFME CLASS IMPLEMENTATION <<<