Merge branch 'setup_doxygen' into 'master'

Setup doxygen

See merge request giacomo.mulas/np_tmcode!3

.gitignore

 build/cluster/*
 build/sphere/*
 build/trapping/*
-
+doc/build/*

doc/src/README.md

+# Folder instructions
+
+This directory contains the material to build the project documentation with `doxygen`.
+
+## Instructions
+
+The project documentation is managed by `doxygen`, a documentation generator that is able to extract documents directly from properly formatted comment sections of the source code. To build a local instance of project documents, make sure that you have `doxygen` installed, then `cd` into the document source folder (the folder containing the `conf.dox` file, specifically `np_tmcode/doc/src`) and finally run:
+```
+doxygen conf.dox
+```
+
+`doxygen` will automatically build the HTML structure to cover all the documented source code and it will additionally provide the fundamental structure to prepare a LaTeX formatted version of the documents. These two outputs will be placed, respectively, under the folders `np_tmcode/doc/build/html` and `np_tmcode/doc/build/latex`.
\ No newline at end of file

src/sphere/List.h

+/*! \file List.h
+ */
+
+#ifndef LIST_OUT_OF_BOUNDS_EXCEPTION
+#define LIST_OUT_OF_BOUNDS_EXCEPTION 1
+#endif
+
+/**
+ * \brief A class to represent dynamic lists.
+ *
+ * This class helps in the creation and management of dynamic lists of
+ * objects, whose size is not known in advance. List offers the advantage
+ * of saving memory, since only the necessary space will be allocated,
+ * but it has the disadvantage of creating an object which is not contiguous
+ * in memory and, therefore, is very inefficient for subsequent manipulation.
+ *
+ * For this reason, the best use of List objects is to collect all the
+ * desired members and then, once the element number is known, to convert
+ * the List to C array, by calling List.to_array(). This function returns
+ * a contiguous array of type T[SIZE] that can be used for indexed access.
+ */
+template<class T> class List {
+ protected:
+  int size; //!< Size of the List.
+  struct element {
+    T value; //!< Value of the list element.
+    element* p_prev; //!< Pointer to the previous element in the list.
+  }; //!< List element connector.
+  element *current, //!< Pointer to element affected by last operation.
+    *first, //!< Pointer to the first element in the List.
+    *last; //!< Pointer to the last element in the List.
+
+ public:
+  /*! \brief List constructor.
+   *
+   * Use the constructor List<T>([int length]) to create a new list with a given
+   * size. If the required size is not known in advance, it is recommended
+   * to create a List with SIZE=1 (this is the default behavior) and then
+   * to append the elements dynamically, using List.append(ELEMENT) (where
+   * ELEMENT needs to be a value of type T, corresponding to the class
+   * template specialization). Note that, due to the default behavior, the
+   * following calls are equivalent and they both produce an integer List
+   * with size equal to 1:
+   *
+   * a = List<int>(1);
+   *
+   * b = List<int>();
+   *
+   * \param length: `int` The size of the list to be constructed [OPTIONAL, default=1].
+   */
+  List(int length = 1) {
+    size = length;
+    first = new element;
+    first->p_prev = NULL;
+    element *current = first;
+    element *p_prev = first;
+    for (int i = 1; i < size; i++) {
+      current = new element;
+      current->p_prev = p_prev;
+      p_prev = current;
+    }
+    last = current;
+  }
+
+  ~List() {
+    current = last;
+    element *old;
+    while (current->p_prev) {
+      old = current;
+      current = old->p_prev;
+      delete old;
+    }
+  }
+  
+  /*! \brief Append an element at the end of the list.
+   *
+   * To dynamically create a list whose size is not known in advance,
+   * elements can be appended in an iterative way. Note that element
+   * manipulation is much more effective in a C array than in a List
+   * object. For this reason, after the List has been created, it is
+   * strongly advised to convert it to a C array by calling the function
+   * List.to_array().
+   *
+   * \param value: `T` The value of the element to be appended.
+   */
+  void append(T value) {
+    element *p_prev = last;
+    current = new element;
+    current->value = value;
+    current->p_prev = p_prev;
+    last = current;
+    size++;
+  }
+
+  /*! \brief Get the element at given index.
+   *
+   * Get the element specified by the index argument. The first element
+   * has index 0 and the last one has index [size - 1].
+   *
+   * \param index: `int` The index of the element to be retrieved. 0 for first.
+   * \return value `T` The value of the element at the requested position.
+   * \throws LIST_OUT_OF_BOUNDS_EXCEPTION: Raised if the index is out of bounds.
+   */
+  T get(int index) {
+    if (index < 0 || index > size - 1) {
+      throw LIST_OUT_OF_BOUNDS_EXCEPTION;
+    }
+    current = last;
+    for (int i = size - 1; i > index; i--) current = current->p_prev;
+    return current->value;
+  }
+
+  /*! \brief Get the number of elements in the list.
+   *
+   * Get the number of elements currently stored in the list.
+   *
+   * \return size `int` The size of the list.
+   */
+  int length() {
+    return size;
+  }
+
+  /*! \brief Set an element by index and value.
+   *
+   * Set the element at the position specified by the index to the value
+   * specified by the value argument.
+   *
+   * \param index: `int` The index of the element to be set. 0 for first.
+   * \param value: `int` The value to store in the pointed element.
+   * \throws LIST_OUT_OF_BOUNDS_EXCEPTION: Raised if the index is out of bounds.
+   */
+  void set(int index, T value) {
+    if (index < 0 || index > size - 1) {
+      throw LIST_OUT_OF_BOUNDS_EXCEPTION;
+    }
+    current = last;
+    for (int i = size - 1; i > index; i--) current = current->p_prev;
+    current->value = value;
+  }
+
+  /*! \brief Convert the list to a C array.
+   *
+   * The List object is useful to dynamically manage a set of objects
+   * when the number of elements is not known in advance. However, the
+   * resulting object is not contiguosly stored in memory. As a result,
+   * access to specific elements in the middle of the list is not very
+   * effective, because the list needs to be walked every time up to
+   * the desired position. In order to avoid this, List.to_array() makes
+   * a conversion from List to C array, returning a contiguous object,
+   * where indexed access can be used.
+   *
+   * \return array `T*` A C array of type T and size equal to the List size.
+   */
+  T* to_array() {
+    T *array = new T[size];
+    current = last;
+    for (int i = size; i > 0; i--) {
+      array[i - 1] = current->value;
+      current = current->p_prev;
+    }
+    return array;
+  }
+};

src/sphere/edfb.cpp

+/*! \file edfb.cpp
+ */
+
+#include <cstdio>
+#include <cmath>
+#include <complex>
+#include <cstring>
+#include <iostream>
+#include <fstream>
+#include "List.h"
+
+using namespace std;
+
+/*! \brief Load a text file as a sequence of strings in memory.
+ *
+ * The configuration of the field expansion code in FORTRAN uses
+ * shared memory access and file I/O operations managed by different
+ * functions. Although this approach could be theoretically replicated,
+ * it is more convenient to handle input and output to distinct files
+ * using specific functions. load_file() helps in the task of handling
+ * input such as configuration files or text data structures that need
+ * to be loaded entirely. The function performs a line-by line scan of
+ * the input file and returns an array of strings that can be later
+ * parsed and ingested by the concerned code blocks. An optional pointer
+ * to integer allows the function to keep track of the number of file
+ * lines that were read, if needed.
+ *
+ * \param file_name: `string` The path of the file to be read.
+ * \param count: `int*` Pointer to an integer recording the number of
+ * read lines [OPTIONAL, default=NULL].
+ * \return array_lines `string*` An array of strings, one for each input
+ * file line.
+ */
+string *load_file(string file_name, int *count);
+
+/*! \brief C++ implementation of EDFB
+ *
+ *  This code aims at replicating the original work-flow in C++.
+ */
+int main(int argc, char **argv) {
+  // Common variables set
+  complex<double> *dc0, ***dc0m;
+  double *ros, **rcf;
+  int *iog, *nshl;
+  double *xiv, *wns, *wls, *pus, *evs, *vss;
+  string vns[5];
+
+  int max_nsh = 0; // A helper variable to set the size of dc0m
+  int ici;
+
+  // Input file reading section
+  int num_lines = 0;
+  int last_read_line = 0; // Keep track of where the input stream was left
+  string *file_lines = load_file("../../test_data/sphere/DEDFB", &num_lines);
+
+  // Configuration code
+  int nsph, ies;
+  sscanf(file_lines[last_read_line].c_str(), " %d %d", &nsph, &ies);
+  if (ies != 0) ies = 1;
+  double exdc, wp, xip;
+  int exdc_exp, wp_exp, xip_exp;
+  int idfc, nxi, instpc, insn;
+  int nsh;
+  sscanf(
+    file_lines[++last_read_line].c_str(),
+    " %9lf D%d %9lf D%d %8lf D%d %d %d %d %d",
+    &exdc, &exdc_exp,
+    &wp, &wp_exp,
+    &xip, &xip_exp,
+    &idfc, &nxi, &instpc, &insn
+  );
+  exdc *= pow(10.0, exdc_exp);
+  wp *= pow(10.0, wp_exp);
+  xip *= pow(10.0, xip_exp);
+
+  FILE *output = fopen("c_OEDFB", "w");
+  // FORTRAN starts subroutine INXI at this point
+  const double pigt = acos(0.0) * 4.0;
+  const double evc = 6.5821188e-16;
+  if (idfc >= 0) {
+    // Not walked by default input data
+    // This part of the code in not tested
+    vss = new double[nxi];
+    xiv = new double[nxi];
+    pus = new double[nxi];
+    evs = new double[nxi];
+    wns = new double[nxi];
+    wls = new double[nxi];
+    if (instpc == 0) { // The variable vector is explicitly defined
+      double vs;
+      int vs_exp;
+      for (int jxi_r = 0; jxi_r < nxi; jxi_r++) {
+	sscanf(file_lines[++last_read_line].c_str(), " %lf D%d", &vs, &vs_exp);
+	vs *= pow(10.0, vs_exp);
+	vss[jxi_r] = vs;
+      }
+      switch (insn) {
+      case 1: //xi vector definition
+	vns[insn - 1] = "XIV";
+	fprintf(output, "  JXI     XIV          WNS          WLS          PUS          EVS\n");
+	for (int jxi210w = 0; jxi210w < nxi; jxi210w++) {
+	  xiv[jxi210w] = vss[jxi210w];
+	  pus[jxi210w] = xiv[jxi210w] * wp;
+	  evs[jxi210w] = pus[jxi210w] * evc;
+	  wns[jxi210w] = pus[jxi210w] / 3.0e8;
+	  wls[jxi210w] = pigt / wns[jxi210w];
+	  fprintf(
+	    output,
+	    "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+	    (jxi210w + 1),
+	    xiv[jxi210w],
+	    wns[jxi210w],
+	    wls[jxi210w],
+	    pus[jxi210w],
+	    evs[jxi210w]
+	  );
+	}
+	break;
+      case 2: //wave number vector definition
+	vns[insn - 1] = "WNS";
+	fprintf(output, "  JXI     WNS          WLS          PUS          EVS          XIV\n");
+	for (int jxi230w = 0; jxi230w < nxi; jxi230w++) {
+	  wns[jxi230w] = vss[jxi230w];
+	  wls[jxi230w] = pigt / wns[jxi230w];
+	  xiv[jxi230w] = 3.0e8 * wns[jxi230w] / wp;
+	  pus[jxi230w] = xiv[jxi230w] * wp;
+	  evs[jxi230w] = pus[jxi230w] * evc;
+	  fprintf(
+	    output,
+	    "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+	    (jxi230w + 1),
+	    wns[jxi230w],
+	    wls[jxi230w],
+	    pus[jxi230w],
+	    evs[jxi230w],
+	    xiv[jxi230w]
+	  );
+	}
+	break;
+      case 3: //wavelength vector definition
+	vns[insn - 1] = "WLS";
+	fprintf(output, "  JXI     WLS          WNS          PUS          EVS          XIV\n");
+	for (int jxi250w = 0; jxi250w < nxi; jxi250w++) {
+	  wls[jxi250w] = vss[jxi250w];
+	  wns[jxi250w] = pigt / wls[jxi250w];
+	  xiv[jxi250w] = 3.0e8 * wns[jxi250w] / wp;
+	  pus[jxi250w] = xiv[jxi250w] * wp;
+	  evs[jxi250w] = pus[jxi250w] * evc;
+	  fprintf(
+	    output,
+	    "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+	    (jxi250w + 1),
+	    wls[jxi250w],
+	    wns[jxi250w],
+	    pus[jxi250w],
+	    evs[jxi250w],
+	    xiv[jxi250w]
+	  );
+	}
+	break;
+      case 4: //pu vector definition
+	vns[insn - 1] = "PUS";
+	fprintf(output, "  JXI     PUS          WNS          WLS          EVS          XIV\n");
+	for (int jxi270w = 0; jxi270w < nxi; jxi270w++) {
+	  pus[jxi270w] = vss[jxi270w];
+	  xiv[jxi270w] = pus[jxi270w] / wp;
+	  wns[jxi270w] = pus[jxi270w] / 3.0e8;
+	  wls[jxi270w] = pigt / wns[jxi270w];
+	  evs[jxi270w] = pus[jxi270w] * evc;
+	  fprintf(
+	    output,
+	    "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+	    (jxi270w + 1),
+	    pus[jxi270w],
+	    wns[jxi270w],
+	    wls[jxi270w],
+	    evs[jxi270w],
+	    xiv[jxi270w]
+	  );
+	}
+	break;
+      case 5: //eV vector definition
+	vns[insn - 1] = "EVS";
+	fprintf(output, "  JXI     EVS          WNS          WLS          PUS          XIV\n");
+	for (int jxi290w = 0; jxi290w < nxi; jxi290w++) {
+	  evs[jxi290w] = vss[jxi290w];
+	  pus[jxi290w] = evs[jxi290w] / evc;
+	  xiv[jxi290w] = pus[jxi290w] / wp;
+	  wns[jxi290w] = pus[jxi290w] / 3.0e8;
+	  wls[jxi290w] = pigt / wns[jxi290w];
+	  fprintf(
+	    output,
+	    "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+	    (jxi290w + 1),
+	    evs[jxi290w],
+	    wns[jxi290w],
+	    wls[jxi290w],
+	    pus[jxi290w],
+	    xiv[jxi290w]
+	  );
+	}
+	break;
+      }
+    } else { // The variable vector needs to be computed in steps
+      double vs, vs_step;
+      int vs_exp, vs_step_exp;
+      sscanf(file_lines[++last_read_line].c_str(), " %lf D%d %lf D%d", &vs, &vs_exp, &vs_step, &vs_step_exp);
+      vs *= pow(10.0, vs_exp);
+      vs_step *= pow(10.0, vs_step_exp);
+      switch (insn) {
+      case 1: //xi vector definition
+	vns[insn - 1] = "XIV";
+	fprintf(output, "  JXI     XIV          WNS          WLS          PUS          EVS\n");
+	for (int jxi110w = 0; jxi110w < nxi; jxi110w++) {
+	  vss[jxi110w] = vs;
+	  xiv[jxi110w] = vss[jxi110w];
+	  pus[jxi110w] = xiv[jxi110w] * wp;
+	  wns[jxi110w] = pus[jxi110w] / 3.0e8;
+	  evs[jxi110w] = pus[jxi110w] * evc;
+	  wls[jxi110w] = pigt / wns[jxi110w];
+	  fprintf(
+	    output,
+	    "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+	    (jxi110w + 1),
+	    xiv[jxi110w],
+	    wns[jxi110w],
+	    wls[jxi110w],
+	    pus[jxi110w],
+	    evs[jxi110w]
+	  );
+	  vs += vs_step;
+	}
+	break;
+      case 2: //wave number vector definition
+	vns[insn - 1] = "WNS";
+	fprintf(output, "  JXI     WNS          WLS          PUS          EVS          XIV\n");
+	for (int jxi130w = 0; jxi130w < nxi; jxi130w++) {
+	  vss[jxi130w] = vs;
+	  wns[jxi130w] = vss[jxi130w];
+	  xiv[jxi130w] = 3.0e8 * wns[jxi130w] / wp;
+	  pus[jxi130w] = xiv[jxi130w] * wp;
+	  wls[jxi130w] = pigt / wns[jxi130w];
+	  evs[jxi130w] = pus[jxi130w] * evc;
+	  fprintf(
+	    output,
+	    "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+	    (jxi130w + 1),
+	    wns[jxi130w],
+	    wls[jxi130w],
+	    pus[jxi130w],
+	    evs[jxi130w],
+	    xiv[jxi130w]
+	  );
+	  vs += vs_step;
+	}
+	break;
+      case 3: //wavelength vector definition
+	vns[insn - 1] = "WLS";
+	fprintf(output, "  JXI     WLS          WNS          PUS          EVS          XIV\n");
+	for (int jxi150w = 0; jxi150w < nxi; jxi150w++) {
+	  vss[jxi150w] = vs;
+	  wls[jxi150w] = vss[jxi150w];
+	  wns[jxi150w] = pigt / wls[jxi150w];
+	  xiv[jxi150w] = 3.0e8 * wns[jxi150w] / wp;
+	  pus[jxi150w] = xiv[jxi150w] * wp;
+	  evs[jxi150w] = pus[jxi150w] * evc;
+	  fprintf(
+	    output,
+	    "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+	    (jxi150w + 1),
+	    wls[jxi150w],
+	    wns[jxi150w],
+	    pus[jxi150w],
+	    evs[jxi150w],
+	    xiv[jxi150w]
+	  );
+	  vs += vs_step;
+	}
+	break;
+      case 4: //pu vector definition
+	vns[insn - 1] = "PUS";
+	fprintf(output, "  JXI     PUS          WNS          WLS          EVS          XIV\n");
+	for (int jxi170w = 0; jxi170w < nxi; jxi170w++) {
+	  vss[jxi170w] = vs;
+	  pus[jxi170w] = vss[jxi170w];
+	  xiv[jxi170w] = pus[jxi170w] / wp;
+	  wns[jxi170w] = pus[jxi170w] / 3.0e8;
+	  wls[jxi170w] = pigt / wns[jxi170w];
+	  evs[jxi170w] = pus[jxi170w] * evc;
+	  fprintf(
+	    output,
+	    "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+	    (jxi170w + 1),
+	    pus[jxi170w],
+	    wns[jxi170w],
+	    wls[jxi170w],
+	    evs[jxi170w],
+	    xiv[jxi170w]
+	  );
+	  vs += vs_step;
+	}
+	break;
+      case 5: //eV vector definition
+	vns[insn - 1] = "EVS";
+	fprintf(output, "  JXI     EVS          WNS          WLS          PUS          XIV\n");
+	for (int jxi190w = 0; jxi190w < nxi; jxi190w++) {
+	  vss[jxi190w] = vs;
+	  evs[jxi190w] = vss[jxi190w];
+	  pus[jxi190w] = evs[jxi190w] / evc;
+	  xiv[jxi190w] = pus[jxi190w] / wp;
+	  wns[jxi190w] = pus[jxi190w] / 3.0e8;
+	  wls[jxi190w] = pigt / wns[jxi190w];
+	  fprintf(
+	    output,
+	    "%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+	    (jxi190w + 1),
+	    evs[jxi190w],
+	    wns[jxi190w],
+	    wls[jxi190w],
+	    pus[jxi190w],
+	    xiv[jxi190w]
+	  );
+	  vs += vs_step;
+	}
+	break;
+      }
+    }
+    // End of the untested code section.
+  } else {
+    if (instpc < 1) {
+      // In this case the XI vector is explicitly defined.
+      // Test input comes this way.
+      double xi, pu, wn;
+      int xi_exp;
+      vns[insn - 1] = "XIV";
+      List<double> xi_vector;
+      sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
+      xi *= pow(10.0, xi_exp);
+      xi_vector.set(0, xi);
+      for (int jxi310 = 1; jxi310 < nxi; jxi310++) {
+	sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
+	xi *= pow(10.0, xi_exp);
+	xi_vector.append(xi);
+      }
+      vss = xi_vector.to_array();
+      xiv = xi_vector.to_array();
+      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", pigt / wn);
+      fprintf(output, "%13.4lE", pu);
+      fprintf(output, "%13.4lE\n", pu * evc);
+      fprintf(output, "  SCALE FACTORS XI\n", pu * evc);
+      for (int jxi6612 = 1; jxi6612 <= nxi; jxi6612++)
+	fprintf(output, "%5d%13.4lE\n", jxi6612, xiv[jxi6612 - 1]);
+      //INXI branch ends here.
+    }
+  }
+  last_read_line++;
+  iog = new int[nsph];
+  for (int i = 0; i < nsph; i++) {
+    sscanf(file_lines[last_read_line].c_str(), " %d", (iog + i));
+  }
+  nshl = new int[nsph];
+  ros = new double[nsph];
+  rcf = new double*[nsph];
+  for (int i113 = 1; i113 <= nsph; i113++) {
+    int i_val;
+    double ros_val;
+    int ros_val_exp;
+    if (iog[i113 - 1] < i113) continue;
+    sscanf(file_lines[++last_read_line].c_str(), " %d %9lf D%d", &i_val, &ros_val, &ros_val_exp);
+    nshl[i113 - 1] = i_val;
+    ros[i113 - 1] = ros_val * pow(10.0, ros_val_exp);
+    nsh = nshl[i113 -1];
+    if (i113 == 1) nsh += ies;
+    if ((nsh + 1) / 2 + ies > max_nsh) max_nsh = (nsh + 1) / 2 + ies;
+    rcf[i113 - 1] = new double[nsh];
+    for (int ns = 0; ns < nsh; ns++) {
+      double ns_rcf;
+      int ns_rcf_exp;
+      sscanf(file_lines[++last_read_line].c_str(), " %8lf D%d", &ns_rcf, &ns_rcf_exp);
+      rcf[i113 -1][ns] = ns_rcf * pow(10.0, ns_rcf_exp);
+    }
+  }
+  // The FORTRAN code writes an auxiliary file in binary format. This should
+  // be avoided or possibly replaced with the use of standard file formats for
+  // scientific use (e.g. FITS).
+  ofstream tedf_file;
+  tedf_file.open("c_TEDF", ofstream::binary);
+  tedf_file.write(reinterpret_cast<char *>(&nsph), sizeof(nsph));
+  for (int iogi = 0; iogi < nsph; iogi++)
+    tedf_file.write(reinterpret_cast<char *>(iog + iogi), sizeof(iog[iogi]));
+  tedf_file.write(reinterpret_cast<char *>(&exdc), sizeof(exdc));
+  tedf_file.write(reinterpret_cast<char *>(&wp), sizeof(wp));
+  tedf_file.write(reinterpret_cast<char *>(&xip), sizeof(xip));
+  tedf_file.write(reinterpret_cast<char *>(&idfc), sizeof(idfc));
+  tedf_file.write(reinterpret_cast<char *>(&nxi), sizeof(nxi));
+  for (int i115 = 1; i115 <= nsph; i115++) {
+    if (iog[i115 - 1] < i115) continue;
+    tedf_file.write(reinterpret_cast<char *>(nshl + i115 - 1), sizeof(nshl[i115 - 1]));
+    tedf_file.write(reinterpret_cast<char *>(ros + i115 - 1), sizeof(ros[i115 - 1]));
+    nsh = nshl[i115 - 1];
+    if (i115 == 1) nsh += ies;
+    for (int ins = 0; ins < nsh; ins++)
+      tedf_file.write(reinterpret_cast<char *>(rcf[i115 - 1] + ins), sizeof(rcf[i115 - 1][ins]));
+  }
+  // Remake the dc0m matrix.
+  dc0m = new complex<double>**[max_nsh];
+  for (int dim1 = 0; dim1 < max_nsh; dim1++) {
+    dc0m[dim1] = new complex<double>*[nsph];
+    for (int dim2 = 0; dim2 < nxi; dim2++) {
+      dc0m[dim1][dim2] = new complex<double>[nxi];
+    }
+  }
+  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;
+      nsh = nshl[i162 - 1];
+      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;
+	int dc0_real_exp, dc0_img_exp;
+	sscanf(file_lines[++last_read_line].c_str(), " (%8lf D%d, %8lf D%d)", &dc0_real, &dc0_real_exp, &dc0_img, &dc0_img_exp);
+	dc0_real *= pow(10.0, dc0_real_exp);
+	dc0_img *= pow(10.0, dc0_img_exp);
+	dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;
+	// 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.
+	tedf_file.write(reinterpret_cast<char *>(&dc0_real), sizeof(dc0_real));
+	tedf_file.write(reinterpret_cast<char *>(&dc0_img), sizeof(dc0_img));
+      }
+    }
+  }
+  tedf_file.close();
+  if (idfc != 0) {
+    fprintf(output, "  DIELECTRIC CONSTANTS\n");
+    for (int i473 = 1; i473 <= nsph; i473++) {
+      if (iog[i473 - 1] != i473) continue;
+      ici = (nshl[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 = dc0m[ic472][i473 - 1][0].real(), dc0_img = dc0m[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 <= nsph; i478++) {
+      if (iog[i478 - 1] != i478) continue;
+      ici = (nshl[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 =  %3c\n", ic477, vns[insn - 1].c_str());
+	for (int jxi476 = 0; jxi476 < nxi; jxi476++) {
+	  double dc0_real = dc0m[ic477 - 1][i478 - 1][jxi476].real();
+	  double dc0_img = dc0m[ic477 - 1][i478 - 1][jxi476].imag();
+	  fprintf(output, "%5d (%12.4lE,%12.4lE)\n", (jxi476 + 1), dc0_real, dc0_img);
+	}
+      }
+    }
+  }
+  fclose(output);
+  return 0;
+}
+
+string *load_file(string file_name, int *count = 0) {
+  fstream input_file(file_name.c_str(), ios::in);
+  List<string> file_lines = List<string>();
+  string line;
+  if (input_file.is_open()) {
+    getline(input_file, line);
+    file_lines.set(0, line);
+    while (getline(input_file, line)) {
+      file_lines.append(line);
+    }
+    input_file.close();
+  }
+  string *array_lines = file_lines.to_array();
+  if (count != 0) *count = file_lines.length();
+  return array_lines;
+}