/*! \file Configuration.cpp
 */

#include <cmath>
#include <complex>
#include <cstdio>
#include <exception>
#include <fstream>
#include <regex>
#include <string>
#include <hdf5.h>

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

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

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

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

using namespace std;

GeometryConfiguration::GeometryConfiguration(
					     int _nsph, int _lm, int _in_pol, int _npnt, int _npntts, int _isam,
					     int _li, int _le, int _mxndm, int _iavm,
					     double *x, double *y, double *z,
					     double in_th_start, double in_th_step, double in_th_end,
					     double sc_th_start, double sc_th_step, double sc_th_end,
					     double in_ph_start, double in_ph_step, double in_ph_end,
					     double sc_ph_start, double sc_ph_step, double sc_ph_end,
					     int _jwtm
					     ) {
  number_of_spheres = _nsph;
  l_max = _lm;
  in_pol = _in_pol;
  npnt = _npnt;
  npntts = _npntts;
  meridional_type = _isam;
  li = _li;
  le = _le;
  mxndm = _mxndm;
  iavm = _iavm;
  in_theta_start = in_th_start;
  in_theta_step = in_th_step;
  in_theta_end = in_th_end;
  in_phi_start = in_ph_start;
  in_phi_step = in_ph_step;
  in_phi_end = in_ph_end;
  sc_theta_start = sc_th_start;
  sc_theta_step = sc_th_step;
  sc_theta_end = sc_th_end;
  sc_phi_start = sc_ph_start;
  sc_phi_step = sc_ph_step;
  sc_phi_end = sc_ph_end;
  jwtm = _jwtm;
  sph_x = x;
  sph_y = y;
  sph_z = z;
}

GeometryConfiguration::~GeometryConfiguration() {
  delete[] sph_x;
  delete[] sph_y;
  delete[] sph_z;
}

GeometryConfiguration* GeometryConfiguration::from_legacy(string file_name) {
  int num_lines = 0;
  int last_read_line = 0;
  string *file_lines;
  string str_target, str_num;
  smatch m;
  try {
    file_lines = load_file(file_name, &num_lines);
  } catch (...) {
    OpenConfigurationFileException ex(file_name);
    throw ex;
  }
  int _nsph = 0, _lm = 0, _in_pol = 0, _npnt = 0, _npntts = 0, _isam = 0;
  int _li = 0, _le = 0, _mxndm = 0, _iavm = 0;
  regex re = regex("-?[0-9]+");
  str_target = file_lines[last_read_line++];
  regex_search(str_target, m, re);
  _nsph = stoi(m.str());
  if (_nsph == 1) {
    for (int ri = 0; ri < 5; ri++) {
      str_target = m.suffix().str();
      regex_search(str_target, m, re);
      if (ri == 0) _lm = stoi(m.str());
      if (ri == 1) _in_pol = stoi(m.str());
      if (ri == 2) _npnt = stoi(m.str());
      if (ri == 3) _npntts = stoi(m.str());
      if (ri == 4) _isam = stoi(m.str());
    }
  } else {
    for (int ri = 0; ri < 8; ri++) {
      str_target = m.suffix().str();
      regex_search(str_target, m, re);
      if (ri == 0) _li = stoi(m.str());
      if (ri == 1) _le = stoi(m.str());
      if (ri == 2) _mxndm = stoi(m.str());
      if (ri == 3) _in_pol = stoi(m.str());
      if (ri == 4) _npnt = stoi(m.str());
      if (ri == 5) _npntts = stoi(m.str());
      if (ri == 6) _iavm = stoi(m.str());
      if (ri == 7) _isam = stoi(m.str());
    }
  }
  double *x, *y, *z;
  x = new double[_nsph];
  y = new double[_nsph];
  z = new double[_nsph];
  if (_nsph == 1) {
    x[0] = 0.0;
    y[0] = 0.0;
    z[0] = 0.0;
  } else {
    for (int i = 0; i < _nsph; i++) {
      str_target = file_lines[last_read_line++];
      re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
      for (int ri = 0; ri < 3; ri++) {
	regex_search(str_target, m, re);
	str_num = regex_replace(m.str(), regex("D"), "e");
	str_num = regex_replace(str_num, regex("d"), "e");
	if (ri == 0) x[i] = stod(str_num);
	if (ri == 1) y[i] = stod(str_num);
	if (ri == 2) z[i] = stod(str_num);
	str_target = m.suffix().str();
      }
    }
  }
  double in_th_start, in_th_end, in_th_step, sc_th_start, sc_th_end, sc_th_step;
  re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
  str_target = file_lines[last_read_line++];
  for (int ri = 0; ri < 6; ri++) {
    regex_search(str_target, m, re);
    str_num = regex_replace(m.str(), regex("D"), "e");
    str_num = regex_replace(str_num, regex("d"), "e");
    if (ri == 0) in_th_start = stod(str_num);
    if (ri == 1) in_th_step = stod(str_num);
    if (ri == 2) in_th_end = stod(str_num);
    if (ri == 3) sc_th_start = stod(str_num);
    if (ri == 4) sc_th_step = stod(str_num);
    if (ri == 5) sc_th_end = stod(str_num);
    str_target = m.suffix().str();
  }

  double in_ph_start, in_ph_end, in_ph_step, sc_ph_start, sc_ph_end, sc_ph_step;
  str_target = file_lines[last_read_line++];
  for (int ri = 0; ri < 6; ri++) {
    regex_search(str_target, m, re);
    str_num = regex_replace(m.str(), regex("D"), "e");
    str_num = regex_replace(str_num, regex("d"), "e");
    if (ri == 0) in_ph_start = stod(str_num);
    if (ri == 1) in_ph_step = stod(str_num);
    if (ri == 2) in_ph_end = stod(str_num);
    if (ri == 3) sc_ph_start = stod(str_num);
    if (ri == 4) sc_ph_step = stod(str_num);
    if (ri == 5) sc_ph_end = stod(str_num);
    str_target = m.suffix().str();
  }

  int _jwtm;
  re = regex("[0-9]+");
  str_target = file_lines[last_read_line++];
  regex_search(str_target, m, re);
  _jwtm = stoi(m.str());
  GeometryConfiguration *conf = new GeometryConfiguration(
							  _nsph, _lm, _in_pol, _npnt, _npntts, _isam,
							  _li, _le, _mxndm, _iavm,
							  x, y, z,
							  in_th_start, in_th_step, in_th_end,
							  sc_th_start, sc_th_step, sc_th_end,
							  in_ph_start, in_ph_step, in_ph_end,
							  sc_ph_start, sc_ph_step, sc_ph_end,
							  _jwtm
							  );
  return conf;
}

ScattererConfiguration::ScattererConfiguration(
					       int nsph,
					       double *scale_vector,
					       int nxi,
					       string variable_name,
					       int *iog_vector,
					       double *ros_vector,
					       int *nshl_vector,
					       double **rcf_vector,
					       int dielectric_func_type,
					       complex<double> ***dc_matrix,
					       bool is_external,
					       double ex,
					       double w,
					       double x
					       ) {
  number_of_spheres = nsph;
  number_of_scales = nxi;
  reference_variable_name = variable_name;
  iog_vec = iog_vector;
  radii_of_spheres = ros_vector;
  nshl_vec = nshl_vector;
  rcf = rcf_vector;
  idfc = dielectric_func_type;
  dc0_matrix = dc_matrix;
  use_external_sphere = is_external;
  exdc = ex;
  wp = w;
  xip = x;
  if (variable_name == "XIV") scale_vec = scale_vector;
  else {
    scale_vec = new double[number_of_scales]();
    const double pi2 = 2.0 * acos(-1.0);
    const double evc = 6.5821188e-16;
    for (int si = 0; si < number_of_scales; si++) {
      double value = scale_vector[si];
      if (variable_name.compare("WNS") == 0) value *= (3.0e8 / wp);
      else if (variable_name.compare("WLS") == 0) value = pi2 / value * 3.0e8 / wp;
      else if (variable_name.compare("PUS") == 0) value /= wp;
      else if (variable_name.compare("EVS") == 0) value /= (evc * wp);
      scale_vec[si] = value;
    }
  }
}

ScattererConfiguration::~ScattererConfiguration() {
  int max_ici = 0;
  for (int i = 1; i <= number_of_spheres; i++) {
    if (iog_vec[i - 1] < i) continue;
    int nsh = nshl_vec[i - 1];
    if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
  }
  for (int i = 0; i < max_ici; i++) {
    for (int j = 0; j < number_of_spheres; j++) {
      delete[] dc0_matrix[i][j];
    }
  }
  delete[] dc0_matrix;
  for (int i = 0; i < number_of_spheres; i++) {
    delete[] rcf[i];
  }
  delete[] nshl_vec;
  delete[] radii_of_spheres;
  delete[] iog_vec;
  delete[] scale_vec;
}

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

ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
  ScattererConfiguration *conf = NULL;
  return conf;
}

ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
  int nsph;
  int *iog;
  double _exdc, _wp, _xip;
  int _idfc, nxi;
  int *nshl_vector;
  double *xi_vec;
  double *ros_vector;
  double **rcf_vector;
  complex<double> ***dc0m;
  
  int max_ici = 0;
  fstream input;
  input.open(file_name.c_str(), ios::in | ios::binary);
  if (input.is_open()) {
    input.read(reinterpret_cast<char *>(&nsph), sizeof(int));
    iog = new int[nsph]();
    for (int i = 0; i < nsph; i++) {
      input.read(reinterpret_cast<char *>(&(iog[i])), sizeof(int));
    }
    input.read(reinterpret_cast<char *>(&_exdc), sizeof(double));
    input.read(reinterpret_cast<char *>(&_wp), sizeof(double));
    input.read(reinterpret_cast<char *>(&_xip), sizeof(double));
    input.read(reinterpret_cast<char *>(&_idfc), sizeof(int));
    input.read(reinterpret_cast<char *>(&nxi), sizeof(int));
    try {
      xi_vec = new double[nxi]();
    } catch (const bad_alloc &ex) {
      throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of scales " + to_string(nxi));
    }
    for (int i = 0; i < nxi; i++) {
      input.read(reinterpret_cast<char *>(&(xi_vec[i])), sizeof(double));
    }
    nshl_vector = new int[nsph]();
    ros_vector = new double[nsph]();
    rcf_vector = new double*[nsph];
    for (int i115 = 1; i115 <= nsph; i115++) {
      if (iog[i115 - 1] < i115) {
	rcf_vector[i115 - 1] = new double[1]();
	continue;
      }
      input.read(reinterpret_cast<char *>(&(nshl_vector[i115 - 1])), sizeof(int));
      input.read(reinterpret_cast<char *>(&(ros_vector[i115 - 1])), sizeof(double));
      int nsh = nshl_vector[i115 - 1];
      if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
      try {
	rcf_vector[i115 - 1] = new double[nsh]();
      } catch (const bad_alloc &ex) {
	throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of layers " + to_string(nsh));
      }
      for (int nsi = 0; nsi < nsh; nsi++) {
	input.read(reinterpret_cast<char *>(&(rcf_vector[i115 - 1][nsi])), sizeof(double));
      }
    }
    dc0m = new complex<double>**[max_ici];
    int dim3 = (_idfc == 0) ? nxi : 1;
    for (int dim1 = 0; dim1 < max_ici; dim1++) {
      dc0m[dim1] = new complex<double>*[nsph];
      for (int dim2 = 0; dim2 < nsph; dim2++) {
	dc0m[dim1][dim2] = new complex<double>[dim3]();
      }
    }
    for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
      if (_idfc != 0 && jxi468 > 1) continue;
      for (int i162 = 1; i162 <= nsph; i162++) {
	if (iog[i162 - 1] < i162) continue;
	int nsh = nshl_vector[i162 - 1];
	int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
	for (int i157 = 0; i157 < ici; i157++) {
	  double dc0_real, dc0_img;
	  input.read(reinterpret_cast<char *>(&dc0_real), sizeof(double));
	  input.read(reinterpret_cast<char *>(&dc0_img), sizeof(double));
	  dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;
	}
      }
    }
    input.close();
  } else { // Opening of the input file did not succeed.
    OpenConfigurationFileException ex(file_name);
    throw ex;
  }
    
  ScattererConfiguration *conf = new ScattererConfiguration(
							    nsph,
							    xi_vec,
							    nxi,
							    "XIV",
							    iog,
							    ros_vector,
							    nshl_vector,
							    rcf_vector,
							    _idfc,
							    dc0m,
							    false,
							    _exdc,
							    _wp,
							    _xip
							    );
  return conf;
}

ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_name) {
  int num_lines = 0;
  int last_read_line = 0;
  string *file_lines;
  regex re;
  smatch m;
  try {
    file_lines = load_file(dedfb_file_name, &num_lines);
  } catch (...) {
    OpenConfigurationFileException ex(dedfb_file_name);
    throw ex;
  }
  int nsph, ies;
  int max_ici = 0;
  re = regex("[0-9]+");
  string str_target = file_lines[last_read_line];
  for (int ri = 0; ri < 2; ri++) {
    regex_search(str_target, m, re);
    if (ri == 0) nsph = stoi(m.str());
    if (ri == 1) ies = stoi(m.str());
    str_target = m.suffix().str();
  }
  if (ies != 0) ies = 1;
  double _exdc, _wp, _xip;
  str_target = file_lines[++last_read_line];
  re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
  for (int ri = 0; ri < 3; ri++) {
    regex_search(str_target, m, re);
    string str_number = m.str();
    str_number = regex_replace(str_number, regex("D"), "e");
    str_number = regex_replace(str_number, regex("d"), "e");
    if (ri == 0) _exdc = stod(str_number);
    if (ri == 1) _wp = stod(str_number);
    if (ri == 2) _xip = stod(str_number);
    str_target = m.suffix().str();
  }
  int _idfc, nxi, instpc, insn;
  re = regex("-?[0-9]+");
  for (int ri = 0; ri < 4; ri++) {
    regex_search(str_target, m, re);
    if (ri == 0) _idfc = stoi(m.str());
    if (ri == 1) nxi = stoi(m.str());
    if (ri == 2) instpc = stoi(m.str());
    if (ri == 3) insn = stoi(m.str());
    str_target = m.suffix().str();
  }

  double *variable_vector;
  string variable_name;

  if (_idfc < 0) { // Diel. functions at XIP value and XI is scale factor
    variable_name = "XIV";
    if (instpc < 1) { // The variable vector is explicitly defined.
      double xi;
      List<double> xi_vector;
      str_target = file_lines[++last_read_line];
      re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
      regex_search(str_target, m, re);
      string str_number = m.str();
      str_number = regex_replace(str_number, regex("D"), "e");
      str_number = regex_replace(str_number, regex("d"), "e");
      xi = stod(str_number);
      xi_vector.set(0, xi);
      for (int jxi310 = 1; jxi310 < nxi; jxi310++) {
	str_target = file_lines[++last_read_line];
	regex_search(str_target, m, re);
	str_number = m.str();
	str_number = regex_replace(str_number, regex("D"), "e");
	str_number = regex_replace(str_number, regex("d"), "e");
	xi = stod(str_number);
	xi_vector.append(xi);
      }
      variable_vector = xi_vector.to_array();
    } else { // instpc >= 1: the variable vector is defined in steps
      double xi, xi_step;
      str_target = file_lines[++last_read_line];
      re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
      regex_search(str_target, m, re);
      for (int ri = 0; ri < 2; ri++) {
	regex_search(str_target, m, re);
	string str_number = m.str();
	str_number = regex_replace(str_number, regex("D"), "e");
	str_number = regex_replace(str_number, regex("d"), "e");
	if (ri == 0) xi = stod(str_number);
	if (ri == 1) xi_step = stod(str_number);
	str_target = m.suffix().str();
      }
      variable_vector = new double[nxi]();
      for (int jxi320 = 0; jxi320 < nxi; jxi320++) {
	variable_vector[jxi320] = xi;
	xi += xi_step;
      }
    }
  } else { // idfc >= 0
    variable_vector = new double[nxi]();
    if (instpc == 0) { // The variable vector is explicitly defined
      double vs;
      for (int jxi_r = 0; jxi_r < nxi; jxi_r++) {
	str_target = file_lines[++last_read_line];
	re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
	regex_search(str_target, m, re);
	string str_number = m.str();
	str_number = regex_replace(str_number, regex("D"), "e");
	str_number = regex_replace(str_number, regex("d"), "e");
	vs = stod(str_number);
	variable_vector[jxi_r] = vs;
      }
      switch (insn) {
      case 1: //xi vector definition
	variable_name = "XIV";
	break;
      case 2: //wave number vector definition
	variable_name = "WNS";
	break;
      case 3: //wavelength vector definition
	variable_name = "WLS";
	break;
      case 4: //pu vector definition
	variable_name = "PUS";
	break;
      case 5: //eV vector definition
	variable_name = "EVS";
	break;
      }
    } else { // The variable vector needs to be computed in steps
      double vs, vs_step;
      str_target = file_lines[++last_read_line];
      re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
      regex_search(str_target, m, re);
      for (int ri = 0; ri < 2; ri++) {
	regex_search(str_target, m, re);
	string str_number = m.str();
	str_number = regex_replace(str_number, regex("D"), "e");
	str_number = regex_replace(str_number, regex("d"), "e");
	if (ri == 0) vs = stod(str_number);
	if (ri == 1) vs_step = stod(str_number);
	str_target = m.suffix().str();
      }
      for (int jxi110w = 0; jxi110w < nxi; jxi110w++) {
	variable_vector[jxi110w] = vs;
	vs += vs_step;
      }
      switch (insn) {
      case 1: //xi vector definition
	variable_name = "XIV";
	break;
      case 2: //wave number vector definition
	variable_name = "WNS";
	break;
      case 3: //wavelength vector definition
	variable_name = "WLS";
	break;
      case 4: //pu vector definition
	variable_name = "PUS";
	break;
      case 5: //eV vector definition
	variable_name = "EVS";
	break;
      }
    }
  }
  int *iog_vector = new int[nsph]();
  double *ros_vector = new double[nsph]();
  double **rcf_vector = new double*[nsph];
  int *nshl_vector = new int[nsph]();
  int filled_iogs = 0;
  re = regex("[0-9]+");
  while (filled_iogs < nsph) {
    string str_target = file_lines[++last_read_line];
    while(regex_search(str_target, m, re)) {
      iog_vector[filled_iogs++] = stoi(m.str());
      str_target = m.suffix().str();
      if (filled_iogs == nsph) break;
    }
  }
  for (int i113 = 1; i113 <= nsph; i113++) {
    int i_val, nsh;
    double ros_val;
    if (iog_vector[i113 - 1] < i113) {
      rcf_vector[i113 - 1] = new double[1]();
      continue;
    }
    re = regex("[0-9]+");
    str_target = file_lines[++last_read_line];
    regex_search(str_target, m, re);
    i_val = stoi(m.str());
    str_target = m.suffix();
    re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
    regex_search(str_target, m, re);
    string str_number = m.str();
    str_number = regex_replace(str_number, regex("D"), "e");
    str_number = regex_replace(str_number, regex("d"), "e");
    ros_val = stod(str_number);
    nshl_vector[i113 - 1] = i_val;
    if (max_ici < (i_val + 1) / 2) max_ici = (i_val + 1) / 2;
    ros_vector[i113 - 1] = ros_val;
    nsh = nshl_vector[i113 - 1];
    if (i113 == 1) nsh += ies;
    rcf_vector[i113 - 1] = new double[nsh]();
    for (int ns = 0; ns < nsh; ns++) {
      double ns_rcf;
      str_target = file_lines[++last_read_line];
      regex_search(str_target, m, re);
      str_number = m.str();
      str_number = regex_replace(str_number, regex("D"), "e");
      str_number = regex_replace(str_number, regex("d"), "e");
      ns_rcf = stod(str_number);
      rcf_vector[i113 - 1][ns] = ns_rcf;
    }
  }
  complex<double> ***dc0m = new complex<double>**[max_ici];
  int dim3 = (_idfc == 0) ? nxi : 1;
  for (int dim1 = 0; dim1 < max_ici; dim1++) {
    dc0m[dim1] = new complex<double>*[nsph];
    for (int dim2 = 0; dim2 < nsph; dim2++) {
      dc0m[dim1][dim2] = new complex<double>[dim3]();
    }
  }
  re = regex("-?[0-9]+\\.[0-9]+([eEdD][-+]?[0-9]+)?");
  for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
    if (_idfc != 0 && jxi468 > 1) continue;
    for (int i162 = 1; i162 <= nsph; i162++) {
      if (iog_vector[i162 - 1] < i162) continue;
      int nsh = nshl_vector[i162 - 1];
      int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
      if (i162 == 1) ici = ici + ies;
      for (int i157 = 0; i157 < ici; i157++) {
	double dc0_real, dc0_img;
	str_target = file_lines[++last_read_line];
	regex_search(str_target, m, re);
	string str_number = m.str();
	str_number = regex_replace(str_number, regex("D"), "e");
	str_number = regex_replace(str_number, regex("d"), "e");
	dc0_real = stod(str_number);
	str_target = m.suffix().str();
	regex_search(str_target, m, re);
	str_number = m.str();
	str_number = regex_replace(str_number, regex("D"), "e");
	str_number = regex_replace(str_number, regex("d"), "e");
	dc0_img = stod(str_number);
	dc0m[i157][i162 - 1][jxi468 - 1] = complex<double>(dc0_real, dc0_img);
      }
    }
  }

  ScattererConfiguration *config = new ScattererConfiguration(
							      nsph,
							      variable_vector,
							      nxi,
							      variable_name,
							      iog_vector,
							      ros_vector,
							      nshl_vector,
							      rcf_vector,
							      _idfc,
							      dc0m,
							      (ies > 0),
							      _exdc,
							      _wp,
							      _xip
							      );
  delete[] file_lines;
  return config;
}

void ScattererConfiguration::print() {
  int ies = (use_external_sphere)? 1 : 0;
  int max_ici = 0;
  printf("### CONFIGURATION DATA ###\n");
  printf("NSPH  = %d\n", number_of_spheres);
  printf("ROS   = [");
  for (int i = 0; i < number_of_spheres; i++) printf("\t%lg", radii_of_spheres[i]);
  printf("\t]\n");
  printf("IOG   = [");
  for (int i = 0; i < number_of_spheres; i++) printf("\t%d", iog_vec[i]);
  printf("\t]\n");
  printf("NSHL  = [");
  for (int i = 0; i < number_of_spheres; i++) printf("\t%d", nshl_vec[i]);
  printf("\t]\n");
  printf("RCF   = [\n");
  for (int i = 1; i <= number_of_spheres; i++) {
    int nsh = nshl_vec[i - 1];
    if (i == 1) nsh += ies;
    if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
    printf("         [");
    for (int ns = 0; ns < nsh; ns++) {
      printf("\t%lg", rcf[i - 1][ns]);
    }
    printf("\t]\n");
  }
  printf("        ]\n");
  printf("SCALE = %s\n", reference_variable_name.c_str());
  printf("NXI   = %d\n", number_of_scales);
  printf("VEC   = [");
  for (int i = 0; i < number_of_scales; i++) printf("\t%lg", scale_vec[i]);
  printf("\t]\n");
  printf("DC0M  = [\n");
  for (int i = 0; i < max_ici; i++) {
    printf("         [\n");
    for (int j = 0; j < number_of_spheres; j++) {
      printf("          [");
      for (int k = 0; k < number_of_scales; k++) {
	if (idfc != 0 and k > 0) continue;
	printf("\t%lg + i(%lg)", dc0_matrix[i][j][k].real(), dc0_matrix[i][j][k].imag());
      }
      printf("\t]\n");
    }
    printf("         ]\n");
  }
  printf("        ]\n");
}

void ScattererConfiguration::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 UnrecognizedConfigurationException(message);
  }
}

void ScattererConfiguration::write_hdf5(string file_name) {
  const double two_pi = acos(0.0) * 4.0;
  const double evc = 6.5821188e-16;
  int ies = (use_external_sphere)? 1 : 0;
  int max_ici = 0;
  List<string> rec_name_list(1);
  List<string> rec_type_list(1);
  List<void *> rec_ptr_list(1);
  string str_type, str_name;
  rec_name_list.set(0, "NSPH");
  rec_type_list.set(0, "INT32_(1)");
  rec_ptr_list.set(0, &number_of_spheres);
  rec_name_list.append("IOGVEC");
  str_type = "INT32_(" + to_string(number_of_spheres) + ")";
  rec_type_list.append(str_type);
  rec_ptr_list.append(iog_vec);
  rec_name_list.append("EXDC");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&exdc);
  rec_name_list.append("WP");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&wp);
  rec_name_list.append("XIP");
  rec_type_list.append("FLOAT64_(1)");
  rec_ptr_list.append(&xip);
  rec_name_list.append("IDFC");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&idfc);
  rec_name_list.append("NXI");
  rec_type_list.append("INT32_(1)");
  rec_ptr_list.append(&number_of_scales);
  rec_name_list.append("XIVEC");
  str_type = "FLOAT64_(" + to_string(number_of_scales) + ")";
  rec_type_list.append(str_type);
  rec_ptr_list.append(scale_vec);
  for (int i115 = 1; i115 <= number_of_spheres; i115++) {
    if (iog_vec[i115 - 1] < i115) continue;
    str_name = "NSHL_" + to_string(i115);
    rec_name_list.append(str_name);
    rec_type_list.append("INT32_(1)");
    rec_ptr_list.append(&(nshl_vec[i115 - 1]));
    str_name = "ROS_" + to_string(i115);
    rec_name_list.append(str_name);
    rec_type_list.append("FLOAT64_(1)");
    rec_ptr_list.append(&(radii_of_spheres[i115 - 1]));
    int nsh = nshl_vec[i115 - 1];
    if (i115 == 1) nsh += ies;
    if (max_ici < (nsh + 1) / 2) max_ici = nsh + 1 / 2;
    str_name = "RCF_" + to_string(i115);
    str_type = "FLOAT64_(" + to_string(nsh) + ")";
    rec_name_list.append(str_name);
    rec_type_list.append(str_type);
    rec_ptr_list.append(&(rcf[i115 - 1][0]));
  }

  int dim3 = (idfc == 0) ? number_of_scales : 1;
  int dc0m_size = 2 * dim3 * number_of_spheres * max_ici;
  double *dc0m = new double[dc0m_size];
  int dc0_index = 0;
  for (int jxi468 = 1; jxi468 <= number_of_scales; jxi468++) {
    if (idfc != 0 && jxi468 > 1) continue;
    for (int i162 = 1; i162 <= number_of_spheres; i162++) {
      if (iog_vec[i162 - 1] < i162) continue;
      int nsh = nshl_vec[i162 - 1];
      int ici = (nsh + 1) / 2;
      if (i162 == 1) ici = ici + ies;
      for (int i157 = 0; i157 < ici; i157++) {
	double dc0_real, dc0_imag;
	dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
	dc0_imag = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
	dc0m[dc0_index++] = dc0_real;
	dc0m[dc0_index++] = dc0_imag;
      }
    }
  }
  str_type = "FLOAT64_(" + to_string(dc0m_size) + ")";
  rec_name_list.append("DC0M");
  rec_type_list.append(str_type);
  rec_ptr_list.append(dc0m);

  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, "c_TEDF.hd5", H5F_ACC_TRUNC);
  for (int ri = 0; ri < rec_num; ri++)
    hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
  hdf_file->close();
  
  // Clean memory
  delete[] dc0m;
  delete[] rec_names;
  delete[] rec_types;
  delete[] rec_pointers;
  delete hdf_file;
}

void ScattererConfiguration::write_legacy(string file_name) {
  const double two_pi = acos(0.0) * 4.0;
  const double evc = 6.5821188e-16;
  fstream output;
  int ies = (use_external_sphere)? 1 : 0;
  output.open(file_name.c_str(), ios::out | ios::binary);
  output.write(reinterpret_cast<char *>(&number_of_spheres), sizeof(int));
  for (int i = 0; i < number_of_spheres; i++)
    output.write(reinterpret_cast<char *>(&(iog_vec[i])), sizeof(int));
  output.write(reinterpret_cast<char *>(&exdc), sizeof(double));
  output.write(reinterpret_cast<char *>(&wp), sizeof(double));
  output.write(reinterpret_cast<char *>(&xip), sizeof(double));
  output.write(reinterpret_cast<char *>(&idfc), sizeof(int));
  output.write(reinterpret_cast<char *>(&number_of_scales), sizeof(int));
  for (int i = 0; i < number_of_scales; i++)
    output.write(reinterpret_cast<char *>(&(scale_vec[i])), sizeof(double));
  for (int i115 = 1; i115 <= number_of_spheres; i115++) {
    if (iog_vec[i115 - 1] < i115) continue;
    output.write(reinterpret_cast<char *>(&(nshl_vec[i115 - 1])), sizeof(int));
    output.write(reinterpret_cast<char *>(&(radii_of_spheres[i115 - 1])), sizeof(double));
    int nsh = nshl_vec[i115 - 1];
    if (i115 == 1) nsh += ies;
    for (int nsi = 0; nsi < nsh; nsi++)
      output.write(reinterpret_cast<char *>(&(rcf[i115 - 1][nsi])), sizeof(double));
  }
  for (int jxi468 = 1; jxi468 <= number_of_scales; jxi468++) {
    if (idfc != 0 && jxi468 > 1) continue;
    for (int i162 = 1; i162 <= number_of_spheres; i162++) {
      if (iog_vec[i162 - 1] < i162) continue;
      int nsh = nshl_vec[i162 - 1];
      int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
      if (i162 == 1) ici = ici + ies;
      for (int i157 = 0; i157 < ici; i157++) {
	double dc0_real, dc0_img;
	dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
	dc0_img = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
	// The FORTRAN code writes the complex numbers as a 16-byte long binary stream.
	// Here we assume that the 16 bytes are equally split in 8 bytes to represent the
	// real part and 8 bytes to represent the imaginary one.
	output.write(reinterpret_cast<char *>(&dc0_real), sizeof(double));
	output.write(reinterpret_cast<char *>(&dc0_img), sizeof(double));
      }
    }
  }
  output.close();
}

void ScattererConfiguration::write_formatted(string file_name) {
  const double evc = 6.5821188e-16;
  const double two_pi = acos(0.0) * 4.0;
  double *xi_vec;
  int ici;
  int ies = (use_external_sphere)? 1: 0;
  FILE *output = fopen(file_name.c_str(), "w");
  int scale_type = -1;
  if (reference_variable_name.compare("XIV") == 0) scale_type = 0;
  else if (reference_variable_name.compare("WNS") == 0) scale_type = 1;
  else if (reference_variable_name.compare("WLS") == 0) scale_type = 2;
  else if (reference_variable_name.compare("PUS") == 0) scale_type = 3;
  else if (reference_variable_name.compare("EVS") == 0) scale_type = 4;
  if (idfc >= 0) { // Dielectric functions are constant or they depend on XI
    double  *pu_vec, *ev_vec, *wn_vec, *wl_vec;
    xi_vec = new double[number_of_scales];
    pu_vec = new double[number_of_scales];
    ev_vec = new double[number_of_scales];
    wn_vec = new double[number_of_scales];
    wl_vec = new double[number_of_scales];
    switch (scale_type) {
    case 0:
      fprintf(output, "  JXI     XIV          WNS          WLS          PUS          EVS\n");
      for (int i = 0; i < number_of_scales; i++) {
	xi_vec[i] = scale_vec[i];
	pu_vec[i] = xi_vec[i] * wp;
	ev_vec[i] = pu_vec[i] * evc;
	wn_vec[i] = pu_vec[i] / 3.0e8;
	wl_vec[i] = two_pi / wn_vec[i];
	fprintf(
		output,
		"%5d%13.4lE%13.4lE%13.4lE%13.4lE%13.4lE\n",
		(i + 1),
		xi_vec[i],
		wn_vec[i],
		wl_vec[i],
		pu_vec[i],
		ev_vec[i]
		);
      }
      break;
    case 1:
      fprintf(output, "  JXI     WNS          WLS          PUS          EVS          XIV\n");
      for (int i = 0; i < number_of_scales; i++) {
	xi_vec[i] = scale_vec[i];
	pu_vec[i] = xi_vec[i] * wp;
	wn_vec[i] = pu_vec[i] / 3.0e8;
	wl_vec[i] = two_pi / wn_vec[i];
	ev_vec[i] = pu_vec[i] * evc;
	fprintf(
		output,
		"%5d%13.4lE%13.4lE%13.4lE%13.4lE%13.4lE\n",
		(i + 1),
		wn_vec[i],
		wl_vec[i],
		pu_vec[i],
		ev_vec[i],
		xi_vec[i]
		);
      }
      break;
    case 2:
      fprintf(output, "  JXI     WLS          WNS          PUS          EVS          XIV\n");
      for (int i = 0; i < number_of_scales; i++) {
	xi_vec[i] = scale_vec[i];
	pu_vec[i] = xi_vec[i] * wp;
	wn_vec[i] = pu_vec[i] / 3.0e8;
	wl_vec[i] = two_pi / wn_vec[i];
	ev_vec[i] = pu_vec[i] * evc;
	fprintf(
		output,
		"%5d%13.4lE%13.4lE%13.4lE%13.4lE%13.4lE\n",
		(i + 1),
		wl_vec[i],
		wn_vec[i],
		pu_vec[i],
		ev_vec[i],
		xi_vec[i]
		);
      }
      break;
    case 3:
      fprintf(output, "  JXI     PUS          WNS          WLS          EVS          XIV\n");
      for (int i = 0; i < number_of_scales; i++) {
	xi_vec[i] = scale_vec[i];
	pu_vec[i] = xi_vec[i] * wp;
	wn_vec[i] = pu_vec[i] / 3.0e8;
	wl_vec[i] = two_pi / wn_vec[i];
	ev_vec[i] = pu_vec[i] * evc;
	fprintf(
		output,
		"%5d%13.4lE%13.4lE%13.4lE%13.4lE%13.4lE\n",
		(i + 1),
		pu_vec[i],
		wn_vec[i],
		wl_vec[i],
		ev_vec[i],
		xi_vec[i]
		);
      }
      break;
    case 4:
      fprintf(output, "  JXI     EVS          WNS          WLS          PUS          XIV\n");
      for (int i = 0; i < number_of_scales; i++) {
	xi_vec[i] = scale_vec[i];
	pu_vec[i] = xi_vec[i] * wp;
	wn_vec[i] = pu_vec[i] / 3.0e8;
	wl_vec[i] = two_pi / wn_vec[i];
	ev_vec[i] = pu_vec[i] * evc;
	fprintf(
		output,
		"%5d%13.4lE%13.4lE%13.4lE%13.4lE%13.4lE\n",
		(i + 1),
		ev_vec[i],
		wn_vec[i],
		wl_vec[i],
		pu_vec[i],
		xi_vec[i]
		);
      }
      break;
    default:
      throw UnrecognizedConfigurationException(
					       "Wonrg parameter set: unrecognized scale variable type "
					       + reference_variable_name
					       );
      break;
    }
  } else { // idfc < 0, Dielectric functions are at XIP and XI is scale for dimensions
    double pu, wn;
    xi_vec = scale_vec;
    pu = xip * wp;
    wn = pu / 3.0e8;
    fprintf(output, "          XIP          WN           WL           PU           EV\n");
    fprintf(output, "     %13.4lE", xip);
    fprintf(output, "%13.4lE", wn);
    fprintf(output, "%13.4lE", two_pi / wn);
    fprintf(output, "%13.4lE", pu);
    fprintf(output, "%13.4lE\n", pu * evc);
    fprintf(output, "  SCALE FACTORS XI\n");
    for (int i = 0; i < number_of_scales; i++)
      fprintf(output, "%5d%13.4lE\n", (i + 1), xi_vec[i]);
  }
  if (idfc != 0) {
    fprintf(output, "  DIELECTRIC CONSTANTS\n");
    for (int i473 = 1; i473 <= number_of_spheres; i473++) {
      if (iog_vec[i473 - 1] != i473) continue;
      ici = (nshl_vec[i473 - 1] + 1) / 2;
      if (i473 == 1) ici += ies;
      fprintf(output, " SPHERE N.%4d\n", i473);
      for (int ic472 = 0; ic472 < ici; ic472++) {
	double dc0_real = dc0_matrix[ic472][i473 - 1][0].real(), dc0_img = dc0_matrix[ic472][i473 - 1][0].imag();
	fprintf(output, "%5d %12.4lE%12.4lE\n", (ic472 + 1), dc0_real, dc0_img);
      }
    }
  } else {
    fprintf(output, "  DIELECTRIC FUNCTIONS\n");
    for (int i478 = 1; i478 <= number_of_spheres; i478++) {
      if (iog_vec[i478 - 1] != i478) continue;
      ici = (nshl_vec[i478 - 1] + 1) / 2;
      if (i478 == 1) ici += ies;
      fprintf(output, " SPHERE N.%4d\n", i478);
      for (int ic477 = 1; ic477 <= ici; ic477++) {
	fprintf(output, " NONTRANSITION LAYER N.%2d, SCALE = %3s\n", ic477, reference_variable_name.c_str());
	for (int jxi476 = 0; jxi476 < number_of_scales; jxi476++) {
	  double dc0_real = dc0_matrix[ic477 - 1][i478 - 1][jxi476].real();
	  double dc0_img = dc0_matrix[ic477 - 1][i478 - 1][jxi476].imag();
	  fprintf(output, "%5d %12.4lE%12.4lE\n", (jxi476 + 1), dc0_real, dc0_img);
	}
      }
    }
  }
  fclose(output);
}
