Something went wrong on our end
Select Git revision
Configuration.cpp
-
Giovanni La Mura authoredGiovanni La Mura authored
Configuration.cpp 24.88 KiB
/*! \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_step, in_th_end,
sc_th_start, sc_th_step, sc_th_end,
in_ph_start, in_ph_step, in_ph_end,
sc_ph_start, sc_ph_step, sc_ph_end,
_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);
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
);
_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);
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);
xi *= pow(10.0, 1.0 * xi_exp);
xi_vector.set(0, xi);
for (int jxi310 = 1; jxi310 < nxi; jxi310++) {
sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
xi *= pow(10.0, 1.0 * xi_exp);
xi_vector.append(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);
xi *= pow(10.0, 1.0 * xi_exp);
xi_step *= pow(10.0, 1.0 * xi_step_exp);
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);
vs *= pow(10.0, 1.0 * vs_exp);
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);
vs *= pow(10.0, 1.0 * vs_exp);
vs_step *= pow(10.0, 1.0 * vs_step_exp);
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));
}
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);
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];
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);
rcf_vector[i113 -1][ns] = ns_rcf * pow(10.0, 1.0 * ns_rcf_exp);
}
}
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);
dc0_real *= pow(10.0, 1.0 * dc0_real_exp);
dc0_img *= pow(10.0, 1.0 * dc0_img_exp);
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);
}