From 2159e85c942a0ca1ce5151cd354d278df28db413 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Sun, 10 Dec 2023 20:21:48 +0100
Subject: [PATCH] Reconfigure doxygen to cover python scripts

---
 doc/src/config.dox                      |   3 +-
 src/scripts/{pycompare => pycompare.py} |  68 +++-
 src/sphere/edfb.cpp                     | 495 ------------------------
 3 files changed, 49 insertions(+), 517 deletions(-)
 rename src/scripts/{pycompare => pycompare.py} (82%)
 delete mode 100644 src/sphere/edfb.cpp

diff --git a/doc/src/config.dox b/doc/src/config.dox
index ee0f8baf..036f8e1b 100644
--- a/doc/src/config.dox
+++ b/doc/src/config.dox
@@ -986,7 +986,8 @@ INPUT_FILE_ENCODING    =
 FILE_PATTERNS          = *.cpp \
 		         *.f \
 		         *.h \
-			 *.md
+			 *.md \
+			 *.py
 
 # The RECURSIVE tag can be used to specify whether or not subdirectories should
 # be searched for input files as well.
diff --git a/src/scripts/pycompare b/src/scripts/pycompare.py
similarity index 82%
rename from src/scripts/pycompare
rename to src/scripts/pycompare.py
index 4adeec21..04224c08 100755
--- a/src/scripts/pycompare
+++ b/src/scripts/pycompare.py
@@ -1,12 +1,32 @@
 #!/bin/python
 
+## @package pycompare
+#  Script to perform output consistency tests
+#
+#  Comparing the numeric output can be rendered hard by the amount of information
+#  contained in a typical output file and the necessity to determine whether a
+#  difference is actually significant or just caused by numeric noise hitting
+#  negligible values. The task of this script is to compare two output files, in
+#  the assumption that they were written by the FORTRAN and the C++ versions of
+#  the code and to flag all the possible inconsistencies according to various
+#  severity levels (namely: NOISE, WARNING, and ERROR).
+
 import re
 
 from math import log10
 from sys import argv
 
+## \cond
 number_reg = re.compile(r'-?[0-9]\.[0-9]+E[-+][0-9]{2,2}')
+## \endcond
 
+## \brief Main execution code
+#
+# `main()` is the function that handles the creation of the script configuration
+# and the execution of the comparison. It returns an integer value corresponding
+# to the number of detected error-level inconsistencies.
+#
+# \returns errors: `int` Number of detected error-level inconsistencies.
 def main():
     config = parse_arguments()
     errors, warnings, noisy = (0, 0, 0)
@@ -32,6 +52,7 @@ def main():
             ))
     return errors
 
+## \brief Perform the comparison of two files.
 def compare_files(config):
     mismatch_count = {
         'errors': 0,
@@ -78,8 +99,14 @@ def compare_files(config):
             l_file.write("  </body>\n")
             l_file.write("</html>\n")
             l_file.close()
+    else:
+        mismatch_count['errors'] = len(c_lines)
+        print("ERROR: {0:s} and {1:s} have different numbers of lines!".format(
+            config['fortran_file_name'], config['c_file_name']
+        ))
     return mismatch_count
 
+## \brief Perform the comparison of two file lines.
 def compare_lines(f_line, c_line, config, line_num=0, num_len=1, log_file=None):
     errors = 0
     warnings = 0
@@ -156,26 +183,21 @@ def compare_lines(f_line, c_line, config, line_num=0, num_len=1, log_file=None):
             log_file.write(log_line + "</code></pre></div>\n")
     return (errors, warnings, noisy)
 
+## \brief Determine the severity of a numerical mismatch.
+#
+#  The severity scale is currently designed with the following integer codes:
+#  0 - the values are equal
+#  1 - the values are subject to suspect numerical noise (green fonts)
+#  2 - the values are different but below error threshold (blue fonts)
+#  3 - the values differ more than error threshold (red fonts)
+#
+#  \param str_f_values: `array(string)` The strings representing the numeric
+#     values read from the FORTRAN output file.
+#  \param str_c_values: `array(string)` The strings representing the numeric
+#     values read from the C++ output file.
+#  \param config: `dict` A dictionary containing the configuration options from
+#     which to read the warning and the error threshold.
 def mismatch_severities(str_f_values, str_c_values, config):
-    """Determine the severity of a numerical mismatch.
-
-       The severiti scale is currently designed with the following integer codes:
-       0 - the values are equal
-       1 - the values are subject to suspect numerical noise (green fonts)
-       2 - the values are different but below error threshold (blue fonts)
-       3 - the values differ more than error threshold (red fonts)
-
-    -----------
-    Parameters:
-    str_f_values: `array(string)`
-        The strings representing the numeric values read from the FORTRAN output
-        file.
-    str_c_values: `array(string)`
-        The strings representing the numeric values read from the C++ output file.
-    config: `dict`
-        A dictionary containing the configuration options from which to read the
-        warning and the error threshold.
-    """
     result = [0 for ri in range(len(str_f_values))]
     for i in range(len(str_f_values)):
         if (str_f_values[i] != str_c_values[i]):
@@ -204,6 +226,7 @@ def mismatch_severities(str_f_values, str_c_values, config):
                 else: result[i] = 3
     return result
     
+## \brief Parse the command line arguments.
 def parse_arguments():
     config = {
         'fortran_file_name': '',
@@ -234,13 +257,14 @@ def parse_arguments():
             raise Exception("Unrecognized argument \'{0:s}\'".format(arg))
     return config
 
+## \brief Print a command-line help summary.
 def print_help():
     print("                                            ")
     print("***              PYCOMPARE               ***")
     print("                                            ")
     print("Compare the output of C++ and FORTRAN codes.")
     print("                                            ")
-    print("Usage: \"./pycompare OPTIONS\"              ")
+    print("Usage: \"./pycompare.py OPTIONS\"           ")
     print("                                            ")
     print("Valid options are:                          ")
     print("--ffile=FORTRAN_OUTPUT   File containing the output of the FORTRAN code (mandatory).")
@@ -253,7 +277,9 @@ def print_help():
     print("                                            ")
     
 
-### PROGRAM EXECUTION ###
+# ### PROGRAM EXECUTION ###
+## \cond
 res = main()
+## \endcond
 if (res > 0): exit(1)
 exit(0)
diff --git a/src/sphere/edfb.cpp b/src/sphere/edfb.cpp
deleted file mode 100644
index 5858e70f..00000000
--- a/src/sphere/edfb.cpp
+++ /dev/null
@@ -1,495 +0,0 @@
-/*! \file edfb.cpp
- */
-
-#include <cstdio>
-#include <cmath>
-#include <complex>
-#include <cstring>
-#include <iostream>
-#include <fstream>
-#include "../include/file_io.h"
-#include "../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 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");
-      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++) {
-    string read_format = "";
-    for (int j = 0; j < i; j++) read_format += " %*d";
-    read_format += " %d";
-    sscanf(file_lines[last_read_line].c_str(), read_format.c_str(), (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;
-    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).
-  int uid = 27;
-  string bin_file_name = "c_TEDF";
-  string status = "UNKNOWN";
-  string mode = "UNFORMATTED";
-  open_file_(&uid, bin_file_name.c_str(), status.c_str(), mode.c_str());
-  write_int_(&uid, &nsph);
-  for (int iogi = 0; iogi < nsph; iogi++)
-    write_int_(&uid, (iog + iogi));
-  write_double_(&uid, &exdc);
-  write_double_(&uid, &wp);
-  write_double_(&uid, &xip);
-  write_int_(&uid, &idfc);
-  write_int_(&uid, &nxi);
-  for (int xivi = 0; xivi < nxi; xivi++)
-    write_double_(&uid, (xiv + xivi));
-  for (int i115 = 1; i115 <= nsph; i115++) {
-    if (iog[i115 - 1] < i115) continue;
-    write_int_(&uid, (nshl + i115 -1));
-    write_double_(&uid, (ros + i115 - 1));
-    nsh = nshl[i115 - 1];
-    if (i115 == 1) nsh += ies;
-    for (int ins = 0; ins < nsh; ins++)
-      write_double_(&uid, (rcf[i115 - 1] + ins));
-  }
-  // Remake the dc0m matrix.
-  dc0m = new complex<double>**[nsph];
-  for (int dim1 = 0; dim1 < nsph; 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.
-	write_complex_(&uid, &dc0_real, &dc0_img);
-      }
-    }
-  }
-  close_file_(&uid);
-  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;
-}
-- 
GitLab