diff --git a/src/Configuration.cpp b/src/Configuration.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4ece13ccb37fcdf9ffbae8ed3b2d77608fda5c0f
--- /dev/null
+++ b/src/Configuration.cpp
@@ -0,0 +1,801 @@
+/*! \file Configuration.cpp
+*/
+
+#include <cmath>
+#include <cstdio>
+#include <fstream>
+#include <string>
+#include "include/List.h"
+#include "include/Parsers.h"
+#include "include/Configuration.h"
+
+using namespace std;
+
+GeometryConfiguration::GeometryConfiguration(
+ int nsph, int lm, int _in_pol, int _npnt, int _npntts, int isam,
+ 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;
+ 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;
+ try {
+ file_lines = load_file(file_name, &num_lines);
+ } catch (exception &ex) {
+ throw OpenConfigurationFileException(file_name);
+ }
+ int nsph, lm, _in_pol, _npnt, _npntts, isam;
+ sscanf(
+ file_lines[last_read_line++].c_str(),
+ " %d %d %d %d %d %d",
+ &nsph, &lm, &_in_pol, &_npnt, &_npntts, &isam
+ );
+ 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++) {
+ double sph_x, sph_y, sph_z;
+ int sph_x_exp, sph_y_exp, sph_z_exp;
+ sscanf(
+ file_lines[last_read_line++].c_str(),
+ " %lf D%d %lf D%d %lf D%d",
+ &sph_x, &sph_x_exp, &sph_y, &sph_y_exp, &sph_z, &sph_z_exp
+ );
+ x[i] = sph_x * pow(10.0, 1.0 * sph_x_exp);
+ y[i] = sph_y * pow(10.0, 1.0 * sph_y_exp);
+ z[i] = sph_z * pow(10.0, 1.0 * sph_z_exp);
+ }
+ }
+ double in_th_start, in_th_end, in_th_step, sc_th_start, sc_th_end, sc_th_step;
+ int in_th_start_exp, in_th_end_exp, in_th_step_exp, sc_th_start_exp, sc_th_end_exp, sc_th_step_exp;
+ sscanf(
+ file_lines[last_read_line++].c_str(),
+ " %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d",
+ &in_th_start, &in_th_start_exp,
+ &in_th_step, &in_th_step_exp,
+ &in_th_end, &in_th_end_exp,
+ &sc_th_start, &sc_th_start_exp,
+ &sc_th_step, &sc_th_step_exp,
+ &sc_th_end, &sc_th_end_exp
+ );
+ in_th_start *= pow(10.0, 1.0 * in_th_start_exp);
+ in_th_step *= pow(10.0, 1.0 * in_th_step_exp);
+ in_th_end *= pow(10.0, 1.0 * in_th_end_exp);
+ sc_th_start *= pow(10.0, 1.0 * sc_th_start_exp);
+ sc_th_step *= pow(10.0, 1.0 * sc_th_step_exp);
+ sc_th_end *= pow(10.0, 1.0 * sc_th_end_exp);
+ double in_ph_start, in_ph_end, in_ph_step, sc_ph_start, sc_ph_end, sc_ph_step;
+ int in_ph_start_exp, in_ph_end_exp, in_ph_step_exp, sc_ph_start_exp, sc_ph_end_exp, sc_ph_step_exp;
+ sscanf(
+ file_lines[last_read_line++].c_str(),
+ " %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d",
+ &in_ph_start, &in_ph_start_exp,
+ &in_ph_step, &in_ph_step_exp,
+ &in_ph_end, &in_ph_end_exp,
+ &sc_ph_start, &sc_ph_start_exp,
+ &sc_ph_step, &sc_ph_step_exp,
+ &sc_ph_end, &sc_ph_end_exp
+ );
+ in_ph_start *= pow(10.0, 1.0 * in_ph_start_exp);
+ in_ph_step *= pow(10.0, 1.0 * in_ph_step_exp);
+ in_ph_end *= pow(10.0, 1.0 * in_ph_end_exp);
+ sc_ph_start *= pow(10.0, 1.0 * sc_ph_start_exp);
+ sc_ph_step *= pow(10.0, 1.0 * sc_ph_step_exp);
+ sc_ph_end *= pow(10.0, 1.0 * sc_ph_end_exp);
+ int _jwtm;
+ sscanf(file_lines[last_read_line++].c_str(), " %d", &_jwtm);
+ GeometryConfiguration *conf = new GeometryConfiguration(
+ nsph, lm, _in_pol, _npnt, _npntts, isam,
+ x, y, z,
+ in_th_start, in_th_end, in_th_step,
+ sc_th_start, sc_th_end, sc_th_step,
+ in_ph_start, in_ph_end, in_ph_step,
+ sc_ph_start, sc_ph_end, sc_ph_step,
+ _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;
+ scale_vec = scale_vector;
+ 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;
+}
+
+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];
+ }
+ }
+ 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) {
+ 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;
+ if (mode.compare("LEGACY") == 0) { // Legacy mode was chosen.
+ 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 (bad_alloc &ex) {
+ throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of scales " + 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) 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 (bad_alloc &ex) {
+ throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of layers " + 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];
+ 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>[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;
+ 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.
+ throw OpenConfigurationFileException(file_name);
+ }
+ } else { // A different binary format was chosen.
+ //TODO: this part is not yet implemented.
+ // Functions to write optimized file formats may be invoked here.
+ }
+ 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;
+ try {
+ file_lines = load_file(dedfb_file_name, &num_lines);
+ } catch (exception &ex) {
+ throw OpenConfigurationFileException(dedfb_file_name);
+ }
+ int nsph, ies;
+ int max_ici = 0;
+ sscanf(file_lines[last_read_line].c_str(), " %d %d", &nsph, &ies);
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ //printf("DEBUG: nsph = %d, ies = %d\n", nsph, ies);
+ if (ies != 0) ies = 1;
+ double _exdc, _wp, _xip;
+ int exdc_exp, wp_exp, xip_exp;
+ int _idfc, nxi, instpc, insn;
+ sscanf(
+ file_lines[++last_read_line].c_str(),
+ " %lf D%d %lf D%d %lf D%d %d %d %d %d",
+ &_exdc, &exdc_exp,
+ &_wp, &wp_exp,
+ &_xip, &xip_exp,
+ &_idfc, &nxi, &instpc, &insn
+ );
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ _exdc *= pow(10.0, 1.0 * 1.0 * exdc_exp);
+ _wp *= pow(10.0, 1.0 * wp_exp);
+ _xip *= pow(10.0, 1.0 * xip_exp);
+ //printf("DEBUG: exdc = %lg, wp = %lg, xip = %lg, idfc = %d, nxi = %d, instpc = %d, insn = %d\n", _exdc, _wp, _xip, _idfc, nxi, instpc, insn);
+
+ 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;
+ int xi_exp;
+ List<double> xi_vector;
+ sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ xi *= pow(10.0, 1.0 * xi_exp);
+ xi_vector.set(0, xi);
+ //printf("DEBUG: xi = %lg\n", xi);
+ for (int jxi310 = 1; jxi310 < nxi; jxi310++) {
+ sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ xi *= pow(10.0, 1.0 * xi_exp);
+ xi_vector.append(xi);
+ //printf("DEBUG: xi = %lg\n", xi);
+ }
+ variable_vector = xi_vector.to_array();
+ } else { // instpc >= 1: the variable vector is defined in steps
+ double xi, xi_step;
+ int xi_exp, xi_step_exp;
+ sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d %9lE D%d", &xi, &xi_exp, &xi_step, &xi_step_exp);
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ xi *= pow(10.0, 1.0 * xi_exp);
+ xi_step *= pow(10.0, 1.0 * xi_step_exp);
+ //printf("DEBUG: xi = %lg, xi_step = %lg\n", xi, xi_step);
+ 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;
+ 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);
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ vs *= pow(10.0, 1.0 * vs_exp);
+ //printf("DEBUG: vs = %lg\n", vs);
+ 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;
+ 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);
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ vs *= pow(10.0, 1.0 * vs_exp);
+ vs_step *= pow(10.0, 1.0 * vs_step_exp);
+ //printf("DEBUG: vs = %lg, vs_step = %lg\n", vs, vs_step);
+ 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;
+ }
+ }
+ }
+ last_read_line++;
+ int *iog_vector = new int[nsph];
+ double *ros_vector = new double[nsph];
+ double **rcf_vector = new double*[nsph];
+ int *nshl_vector = 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_vector + i));
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ //printf("DEBUG: iog [%d] = %d\n", i, iog_vector[i]);
+ }
+ for (int i113 = 1; i113 <= nsph; i113++) {
+ int i_val, nsh;
+ double ros_val;
+ int ros_val_exp;
+ if (iog_vector[i113 - 1] < i113) continue;
+ sscanf(file_lines[++last_read_line].c_str(), " %d %lf D%d", &i_val, &ros_val, &ros_val_exp);
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ 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 * pow(10.0, 1.0 * ros_val_exp);
+ nsh = nshl_vector[i113 - 1];
+ //printf("DEBUG: nshl_vector[%d] = %d, ros_vector[%d] = %lg\n", i113 - 1, i_val, i113 - 1, ros_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;
+ int ns_rcf_exp;
+ sscanf(file_lines[++last_read_line].c_str(), " %lf D%d", &ns_rcf, &ns_rcf_exp);
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ rcf_vector[i113 -1][ns] = ns_rcf * pow(10.0, 1.0 * ns_rcf_exp);
+ //printf("DEBUG: rcf_vector[%d][%d] = %lg\n", i113-1, ns, rcf_vector[i113 -1][ns]);
+ }
+ }
+ //printf("DEBUG: max_ici = %d\n", max_ici);
+ complex<double> ***dc0m = new complex<double>**[max_ici];
+ 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>[nxi];
+ }
+ }
+ 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;
+ int dc0_real_exp, dc0_img_exp;
+ sscanf(file_lines[++last_read_line].c_str(), " (%lf D%d, %lf D%d)", &dc0_real, &dc0_real_exp, &dc0_img, &dc0_img_exp);
+ //printf("DEBUG: %s\n", file_lines[last_read_line].c_str());
+ dc0_real *= pow(10.0, 1.0 * dc0_real_exp);
+ dc0_img *= pow(10.0, 1.0 * dc0_img_exp);
+ //printf("DEBUG: dc0m[%d][%d] = %lg + i(%lg)\n", i157, i162-1, dc0_real, dc0_img);
+ dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * 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
+ );
+ 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) {
+ const double two_pi = acos(0.0) * 4.0;
+ const double evc = 6.5821188e-16;
+ int max_ici = 0;
+ if (mode.compare("LEGACY") == 0) { // Legacy mode was chosen.
+ fstream output;
+ int ies = (use_external_sphere)? 1 : 0;
+ double *xi_vec;
+ if (reference_variable_name.compare("XIV") == 0) xi_vec = scale_vec;
+ else {
+ xi_vec = new double[number_of_scales];
+ if (reference_variable_name.compare("WNS") == 0) {
+ for (int i = 0; i < number_of_scales; i++)
+ xi_vec[i] = 3.0e8 * scale_vec[i] / wp;
+ } else if (reference_variable_name.compare("WLS") == 0) {
+ for (int i = 0; i < number_of_scales; i++) {
+ double wn = two_pi / scale_vec[i];
+ xi_vec[i] = 3.0e8 * wn / wp;
+ }
+ } else if (reference_variable_name.compare("PUS") == 0) {
+ for (int i = 0; i < number_of_scales; i++)
+ xi_vec[i] = scale_vec[i] / wp;
+ } else if (reference_variable_name.compare("EVS") == 0) {
+ for (int i = 0; i < number_of_scales; i++) {
+ double pu = scale_vec[i] / evc;
+ xi_vec[i] = pu / wp;
+ }
+ } else {
+ throw UnrecognizedConfigurationException(
+ "Wrong parameter set: unrecognized scale type "
+ + reference_variable_name
+ );
+ }
+ }
+ 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 *>(&(xi_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;
+ if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
+ 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++) {
+ wn_vec[i] = scale_vec[i];
+ wl_vec[i] = two_pi / wn_vec[i];
+ xi_vec[i] = 3.0e8 * wn_vec[i] / wp;
+ pu_vec[i] = xi_vec[i] * wp;
+ 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++) {
+ wl_vec[i] = scale_vec[i];
+ wn_vec[i] = two_pi / wl_vec[i];
+ xi_vec[i] = 3.0e8 * wn_vec[i] / wp;
+ pu_vec[i] = xi_vec[i] * wp;
+ 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++) {
+ pu_vec[i] = scale_vec[i];
+ xi_vec[i] = pu_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++) {
+ ev_vec[i] = scale_vec[i];
+ pu_vec[i] = ev_vec[i] / evc;
+ xi_vec[i] = pu_vec[i] / wp;
+ 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),
+ 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);
+}
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
new file mode 100644
index 0000000000000000000000000000000000000000..c4605423d251488cd3f9cc12fb2f23cf2008df4e
--- /dev/null
+++ b/src/include/Configuration.h
@@ -0,0 +1,321 @@
+/*! \file Configuration.h
+ */
+
+#ifndef INCLUDE_CONFIGURATION_H_
+#define INCLUDE_CONFIGURATION_H_
+
+#include <complex>
+#include <exception>
+#include <string>
+
+/**
+ * \brief Exception for open file error handlers.
+ */
+class OpenConfigurationFileException: public std::exception {
+protected:
+ //! \brief Name of the file that was accessed.
+ std::string file_name;
+
+public:
+ /**
+ * \brief Exception instance constructor.
+ *
+ * \param name: `string` Name of the file that was accessed.
+ */
+ OpenConfigurationFileException(std::string name) { file_name = name; }
+
+ /**
+ * \brief Exception message.
+ */
+ virtual const char* what() const throw() {
+ std::string message = "Error opening configuration file: " + file_name;
+ return message.c_str();
+ }
+};
+
+/**
+ * \brief Exception for unrecognized configuration data sets.
+ */
+class UnrecognizedConfigurationException: public std::exception {
+protected:
+ //! Description of the problem.
+ std::string message;
+public:
+ /**
+ * \brief Exception instance constructor.
+ *
+ * \param problem: `string` Description of the problem that occurred.
+ */
+ UnrecognizedConfigurationException(std::string problem) { message = problem; }
+ /**
+ * \brief Exception message.
+ */
+ virtual const char* what() const throw() {
+ return message.c_str();
+ }
+};
+
+/**
+ * \brief A class to represent the configuration of the scattering geometry.
+ *
+ * GeometryConfiguration is a class designed to store the necessary configuration
+ * data to describe the scattering geometry, including the distribution of the
+ * particle components, the orientation of the incident and scattered radiation
+ * fields and their polarization properties.
+ */
+class GeometryConfiguration {
+ //! Temporary work-around to allow sphere() peeking in.
+ friend void sphere();
+protected:
+ //! \brief Number of spherical components.
+ int number_of_spheres;
+ //! \brief Maximum expansion order of angular momentum.
+ int l_max;
+ //! \brief Incident field polarization status (0 - linear, 1 - circular).
+ int in_pol;
+ //! \brief Number of transition points. QUESTION: correct?
+ int npnt;
+ //! \brief Transition smoothness. QUESTION: correct?
+ int npntts;
+ //! \brief Type of meridional plane definition.
+ int meridional_type;
+ //! \brief Transition matrix layer ID. QUESTION: correct?
+ int jwtm;
+ //! \brief Incident field initial azimuth.
+ double in_theta_start;
+ //! \brief Incident field azimuth step.
+ double in_theta_step;
+ //! \brief Incident field final azimuth.
+ double in_theta_end;
+ //! \brief Scattered field initial azimuth.
+ double sc_theta_start;
+ //! \brief Scattered field azimuth step.
+ double sc_theta_step;
+ //! \brief Scattered field final azimuth.
+ double sc_theta_end;
+ //! \brief Incident field initial elevation.
+ double in_phi_start;
+ //! \brief Incident field elevation step.
+ double in_phi_step;
+ //! \brief Incident field final elevation.
+ double in_phi_end;
+ //! \brief Scattered field initial elevation.
+ double sc_phi_start;
+ //! \brief Scattered field elevation step.
+ double sc_phi_step;
+ //! \brief Scattered field final elevation.
+ double sc_phi_end;
+ //! \brief Vector of spherical components X coordinates.
+ double *sph_x;
+ //! \brief Vector of spherical components Y coordinates.
+ double *sph_y;
+ //! \brief Vector of spherical components Z coordinates.
+ double *sph_z;
+
+public:
+ /*! \brief Build a scattering geometry configuration structure.
+ *
+ * \param nsph: `int` Number of spheres to be used in calculation.
+ * \param lm: `int` Maximum field angular momentum expansion order.
+ * \param in_pol: `int` Incident field polarization status
+ * \param npnt: `int` Number of transition points. QUESTION: correct?
+ * \param npntts: `int` Transition smoothness. QUESTION: correct?
+ * \param meridional_type: `int` Type of meridional plane definition (<0
+ * for incident angles, 0 if determined by incidence and observation, =1
+ * accross z-axis for incidence and observation, >1 across z-axis as a
+ * function of incidence angles for fixed scattering).
+ * \param x: `double*` Vector of spherical components X coordinates.
+ * \param y: `double*` Vector of spherical components Y coordinates.
+ * \param z: `double*` Vector of spherical components Z coordinates.
+ * \param in_th_start: `double` Incident field starting azimuth angle.
+ * \param in_th_step: `double` Incident field azimuth angle step.
+ * \param in_th_end: `double` Incident field final azimuth angle.
+ * \param sc_th_start: `double` Scattered field starting azimuth angle.
+ * \param sc_th_step: `double` Scattered field azimuth angle step.
+ * \param sc_th_end: `double` Scattered field final azimuth angle.
+ * \param in_ph_start: `double` Incident field starting elevation angle.
+ * \param in_ph_step: `double` Incident field elevation angle step.
+ * \param in_ph_end: `double` Incident field final elevation angle.
+ * \param sc_ph_start: `double` Scattered field starting elevation angle.
+ * \param sc_ph_step: `double` Scattered field elevation angle step.
+ * \param sc_ph_end: `double` Scattered field final elevation angle.
+ * \param jwtm: `int` Transition Matrix layer ID. QUESTION: correct?
+ */
+ GeometryConfiguration(
+ int nsph, int lm, int in_pol, int npnt, int npntts, int meridional_type,
+ 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
+ );
+
+ /*! \brief Destroy a GeometryConfiguration instance.
+ */
+ ~GeometryConfiguration();
+
+ /*! \brief Build geometry configuration from legacy configuration input file.
+ *
+ * To allow for consistency tests and backward compatibility, geometry
+ * configurations can be built from legacy configuration files. This function
+ * replicates the approach implemented by the FORTRAN SPH and CLU codes, but
+ * using a C++ oriented work-flow.
+ *
+ * \param file_name: `string` Name of the legacy configuration data file.
+ * \return config: `GeometryConfiguration*` Pointer to object containing the
+ * configuration data.
+ */
+ static GeometryConfiguration *from_legacy(std::string file_name);
+};
+
+/**
+ * \brief A class to represent scatterer configuration objects.
+ *
+ * ScattererConfiguration is a class designed to store the necessary configuration
+ * data to describe the scatterer properties.
+ */
+class ScattererConfiguration {
+ //! Temporary work-around to allow sphere() peeking in.
+ friend void sphere();
+protected:
+ //! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x LAYERS].
+ std::complex<double> ***dc0_matrix;
+ //! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].
+ double *radii_of_spheres;
+ //! \brief Matrix of fractional transition radii with size [N_SPHERES x LAYERS].
+ double **rcf;
+ //! \brief Vector of sphere ID numbers, with size [N_SPHERES].
+ int *iog_vec;
+ //! \brief Vector of layer numbers for every sphere, with size [N_SPHERES].
+ int *nshl_vec;
+ //! \brief Vector of scale parameters, with size [N_SCALES].
+ double *scale_vec;
+ //! \brief Name of the reference variable type (one of XIV, WNS, WLS, PUS, EVS).
+ std::string reference_variable_name;
+ //! \brief Number of spherical components.
+ int number_of_spheres;
+ //! \brief Number of scales to use in calculation.
+ int number_of_scales;
+ //! \brief Type of dielectric functions (<0 at XIP, =0 as function of XI, >0 contants).
+ int idfc;
+ //! \brief EXDC. QUESTION: better definition?
+ double exdc;
+ //! \brief WP. QUESTION: better definition?
+ double wp;
+ //! \brief Peak XI. QUESTION: correct?
+ double xip;
+ //! \brief Flag to control whether to add an external layer.
+ bool use_external_sphere;
+
+public:
+ /*! \brief Build a scatterer configuration structure.
+ *
+ * Prepare a default configuration structure by allocating the necessary
+ * memory structures.
+ *
+ * \param nsph: `int` The number of spheres in the simulation.
+ * \param scale_vector: `double*` The radiation-particle scale vector.
+ * \param nxi: `int` The number of radiation-particle scalings.
+ * \param variable_name: `string` The name of the radiation-particle scaling type.
+ * \param iog_vector: `int*` Array of sphere identification numbers. QUESTION: correct?
+ * \param ros_vector: `double*` Sphere radius array.
+ * \param nshl_vector: `int*` Array of layer numbers.
+ * \param rcf_vector: `double**` Array of fractional break radii. QUESTION: correct?
+ * \param dielectric_func_type: `int` Type of dielectric function definition (=0 for constant,
+ * \>0 as function of scale parameter, <0 for functions at XIP value and XI is scale factor
+ * for dimensions).
+ * \param dc_matrix: Matrix of reference dielectric constants.
+ * \param has_external: `bool` Flag to set whether to add an external spherical layer.
+ * \param exdc: `double` EXDC
+ * \param wp: `double` wp
+ * \param xip: `double` xip
+ */
+ ScattererConfiguration(
+ int nsph,
+ double *scale_vector,
+ int nxi,
+ std::string variable_name,
+ int *iog_vector,
+ double *ros_vector,
+ int *nshl_vector,
+ double **rcf_vector,
+ int dielectric_func_type,
+ std::complex<double> ***dc_matrix,
+ bool has_external,
+ double exdc,
+ double wp,
+ double xip
+ );
+
+ /*! \brief Destroy a scatterer configuration instance.
+ */
+ ~ScattererConfiguration();
+
+ /*! \brief Build configuration from binary configuration input file.
+ *
+ * The configuration step can save configuration data as a binary file. The original
+ * FORTRAN code used this possibility to manage communication between the configuring
+ * code and the calculation program. This possibility is maintained, in case the
+ * configuration step needs to be separated from the calculation execution. In this
+ * case, `from_binary()` is the class method that restores a ScattererConfiguration
+ * object from a previously saved binary file.
+ *
+ * \param file_name: `string` Name of the binary configuration data file.
+ * \param mode: `string` Binary encoding. Can be one of "LEGACY", ... . Optional
+ * (default is "LEGACY").
+ * \return config: `ScattererConfiguration*` Pointer to object containing the
+ * scatterer configuration data.
+ */
+ static ScattererConfiguration* from_binary(std::string file_name, std::string mode = "LEGACY");
+
+ /*! \brief Build scatterer configuration from legacy configuration input file.
+ *
+ * To allow for consistency tests and backward compatibility, ScattererConfiguration
+ * objects can be built from legacy configuration files. This function replicates
+ * the approach implemented by the FORTRAN EDFB code, but using a C++ oriented
+ * work-flow.
+ *
+ * \param file_name: `string` Name of the legacy configuration data file.
+ * \return config: `ScattererConfiguration*` Pointer to object containing the
+ * scatterer configuration data.
+ */
+ static ScattererConfiguration* from_dedfb(std::string file_name);
+
+ /*! \brief Print the contents of the configuration object to terminal.
+ *
+ * In case of quick debug testing, `ScattererConfiguration.print()` allows printing
+ * a formatted summary of the configuration data to terminal.
+ */
+ void print();
+
+ /*! \brief Write the scatterer configuration data to binary output.
+ *
+ * The execution work-flow may be split in a configuration step and one or more
+ * calculation steps. In case the calculation is not being run all-in-one, it can
+ * be useful to save the configuration data. `ScattererConfiguration.write_binary()`
+ * performs the operation of saving the configuration in binary format. This function
+ * can work in legacy mode, to write backward compatible configuration files, as well
+ * as by wrapping the data into common scientific formats (NB: this last option still
+ * needs to be implemented).
+ *
+ * \param file_name: `string` Name of the file to be written.
+ * \param mode: `string` Binary encoding. Can be one of "LEGACY", ... . Optional
+ * (default is "LEGACY").
+ */
+ void write_binary(std::string file_name, std::string mode = "LEGACY");
+
+ /*! \brief Write the scatterer configuration data to formatted text output.
+ *
+ * Writing configuration to formatted text is an optional operation, which may turn
+ * out to be useful for consistency checks. As a matter of fact, formatted configuration
+ * output is not read back by the FORTRAN code work-flow and it can be safely omitted,
+ * unless there is a specific interest in assessing that the legacy code and the
+ * updated one are doing the same thing.
+ *
+ * \param file_name: `string` Name of the file to be written.
+ */
+ void write_formatted(std::string file_name);
+};
+
+#endif