Skip to content
Snippets Groups Projects
Commit 26321e08 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Create a module for configuration data

parent 2e23f636
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/Configuration.cpp 0 → 100644
+ 801
− 0
View file @ 26321e08
/*! \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);
}
src/include/Configuration.h 0 → 100644
+ 321
− 0
View file @ 26321e08
/*! \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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment