From 26321e08fffe7b9629311e8878754162d70a3a12 Mon Sep 17 00:00:00 2001 From: Giovanni La Mura <giovanni.lamura@inaf.it> Date: Thu, 19 Oct 2023 12:20:42 +0200 Subject: [PATCH] Create a module for configuration data --- src/Configuration.cpp | 801 ++++++++++++++++++++++++++++++++++++ src/include/Configuration.h | 321 +++++++++++++++ 2 files changed, 1122 insertions(+) create mode 100644 src/Configuration.cpp create mode 100644 src/include/Configuration.h diff --git a/src/Configuration.cpp b/src/Configuration.cpp new file mode 100644 index 00000000..4ece13cc --- /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 00000000..c4605423 --- /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 -- GitLab