Skip to content
Commits on Source (3)
Expand all Hide whitespace changes
Inline Side-by-side
src/cluster/cluster.cpp
View file @ b07fdabd
......@@ -53,7 +53,7 @@ using namespace std;
// I would like to put it all in a struct, but then I'd have to write a constructor for it, due to members defined as references, creating a worse nightmare than the one I'd like to simplify...
int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, int nxi, int nsph, np_int mxndm, int inpol, int iavm, int npnt, int npntts, int isam, int lm, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double sqsfi, double exri, int lcalc, dcomplex arg, double wn, double vk, np_int ndit, int isq, int ibf, Logger *logger);
int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double sqsfi, dcomplex arg, double wn, double vk, Logger *logger);
/*! \brief C++ implementation of CLU
*
......@@ -93,26 +93,13 @@ void cluster(const string& config_file, const string& data_file, const string& o
exit(1);
}
logger->log(" done.\n", LOG_INFO);
int s_nsph = sconf->get_param("nsph");
int nsph = gconf->get_param("nsph");
int s_nsph = sconf->number_of_spheres;
int nsph = gconf->number_of_spheres;
if (s_nsph == nsph) {
// Shortcuts to variables stored in configuration objects
np_int mxndm = (np_int)gconf->get_param("mxndm");
int inpol = (int)gconf->get_param("in_pol");
int npnt = (int)gconf->get_param("npnt");
int npntts = (int)gconf->get_param("npntts");
int iavm = (int)gconf->get_param("iavm");
int isam = (int)gconf->get_param("meridional_type");
int nxi = (int)sconf->get_param("number_of_scales");
int idfc = (int)sconf->get_param("idfc");
ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);
double wp = sconf->get_param("wp");
double wp = sconf->wp;
// Global variables for CLU
int lm = (int)gconf->get_param("l_max");
int li = (int)gconf->get_param("li");
int le = (int)gconf->get_param("le");
if (li > lm) lm = li;
if (le > lm) lm = le;
C1 *c1 = new C1(gconf, sconf);
C3 *c3 = new C3();
C4 *c4 = new C4(gconf);
......@@ -120,18 +107,15 @@ void cluster(const string& config_file, const string& data_file, const string& o
// End of add-ons initialization
C6 *c6 = new C6(c4->lmtpo);
FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
int jer = 0;
int lcalc = 0;
dcomplex arg = 0.0 + 0.0 * I;
dcomplex ccsam = 0.0 + 0.0 * I;
int configurations = (int)sconf->get_param("configurations");
C2 *c2 = new C2(gconf, sconf);
np_int ndit = 2 * nsph * c4->nlim;
const int ndi = c4->nsph * c4->nlim;
np_int ndit = 2 * ndi;
logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");
time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");
const int ndi = c4->nsph * c4->nlim;
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
double *gaps = new double[nsph]();
double *gaps = new double[gconf->number_of_spheres]();
double *tqev = new double[3]();
double *tqsv = new double[3]();
double **tqse, **tqss, **tqce, **tqcs;
......@@ -145,10 +129,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
tqcs = new double*[2];
tqcps = new dcomplex*[2];
for (int ti = 0; ti < 2; ti++) {
tqse[ti] = new double[nsph]();
tqspe[ti] = new dcomplex[nsph]();
tqss[ti] = new double[nsph]();
tqsps[ti] = new dcomplex[nsph]();
tqse[ti] = new double[gconf->number_of_spheres]();
tqspe[ti] = new dcomplex[gconf->number_of_spheres]();
tqss[ti] = new double[gconf->number_of_spheres]();
tqsps[ti] = new dcomplex[gconf->number_of_spheres]();
tqce[ti] = new double[3]();
tqcpe[ti] = new dcomplex[3]();
tqcs[ti] = new double[3]();
......@@ -192,13 +176,12 @@ void cluster(const string& config_file, const string& data_file, const string& o
cmullr[ci] = new double[4]();
cmul[ci] = new double[4]();
}
int isq, ibf;
int jwtm = (int)gconf->get_param("jwtm");
double scan, cfmp, sfmp, cfsp, sfsp;
// End of global variables for CLU
fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
fprintf(output, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam
nsph, c4->li, c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,
gconf->npntts, gconf->iavm, gconf->iavm
);
fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
for (int ri = 0; ri < nsph; ri++) fprintf(output, "%17.8lE%17.8lE%17.8lE\n",
......@@ -219,7 +202,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
p_scattering_angles->phsstp, p_scattering_angles->phslst
);
fprintf(output, " READ(IR,*)JWTM\n");
fprintf(output, " %5d\n", jwtm);
fprintf(output, " %5d\n", gconf->jwtm);
fprintf(output, " READ(ITIN)NSPHT\n");
fprintf(output, " READ(ITIN)(IOG(I),I=1,NSPH)\n");
fprintf(output, " READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
......@@ -239,9 +222,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
}
}
thdps(c4->lm, zpv);
double exdc = sconf->get_param("exdc");
double exdc = sconf->exdc;
double exri = sqrt(exdc);
double vk = 0.0; // NOTE: Not really sure it should be initialized at 0
double vk = 0.0;
fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
fstream *tppoanp = new fstream;
fstream &tppoan = *tppoanp;
......@@ -253,6 +236,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
#else
logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);
#endif
int iavm = gconf->iavm;
int isam = gconf->isam;
int inpol = gconf->in_pol;
int nxi = sconf->number_of_scales;
int nth = p_scattering_angles->nth;
int nths = p_scattering_angles->nths;
int nph = p_scattering_angles->nph;
......@@ -266,9 +253,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
double wn = wp / 3.0e8;
double xip = sconf->get_param("xip");
double xip = sconf->xip;
double sqsfi = 1.0;
if (idfc < 0) {
if (sconf->idfc < 0) {
vk = xip * wn;
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
fprintf(output, " \n");
......@@ -276,12 +263,11 @@ void cluster(const string& config_file, const string& data_file, const string& o
// do the first iteration on jxi488 separately, since it seems to be different from the others
int jxi488 = 1;
chrono::time_point<chrono::high_resolution_clock> start_iter_1 = chrono::high_resolution_clock::now();
jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, c1, c1ao, c2, c3, c4, c6, c9, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, argi, args, gap, gapp, tqce, tqcpe, tqcs, tqcps, duk, tppoan, cextlr, cext, cmullr, cmul, gapv, tqev, tqsv, nxi, nsph, mxndm, inpol, iavm, npnt, npntts, isam, lm, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, lcalc, arg, wn, vk, ndit, isq, ibf, logger);
int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, c1, c1ao, c2, c3, c4, c6, c9, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, argi, args, gap, gapp, tqce, tqcpe, tqcs, tqcps, duk, tppoan, cextlr, cext, cmullr, cmul, gapv, tqev, tqsv, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, arg, wn, vk, logger);
chrono::time_point<chrono::high_resolution_clock> end_iter_1 = chrono::high_resolution_clock::now();
elapsed = end_iter_1 - start_iter_1;
message = "INFO: First iteration took " + to_string(elapsed.count()) + "s.\n";
logger->log(message);
logger->log("Finished: output written to " + output_path + "/c_OCLU\n");
time_logger->log(message);
// Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled
......@@ -297,8 +283,6 @@ void cluster(const string& config_file, const string& data_file, const string& o
if (myompthread == 0) ompnumthreads = omp_get_num_threads();
#endif
// To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
ScattererConfiguration *sconf_2 = NULL;
GeometryConfiguration *gconf_2 = NULL;
C1 *c1_2 = NULL;
C1_AddOns *c1ao_2 = NULL;
C2 *c2_2 = NULL;
......@@ -332,15 +316,6 @@ void cluster(const string& config_file, const string& data_file, const string& o
double *gapv_2 = NULL;
double *tqev_2 = NULL;
double *tqsv_2 = NULL;
int nxi_2 = nxi;
int nsph_2 = nsph;
np_int mxndm_2 = mxndm;
int inpol_2 = inpol;
int iavm_2 = iavm;
int npnt_2 = npnt;
int npntts_2 = npntts;
int isam_2 = isam;
int lm_2 = lm;
double *u_2 = NULL;
double *us_2 = NULL;
double *un_2 = NULL;
......@@ -357,18 +332,11 @@ void cluster(const string& config_file, const string& data_file, const string& o
double cfsp_2 = cfsp;
double sfsp_2 = sfsp;
double sqsfi_2 = sqsfi;
double exri_2 = exri;
int lcalc_2 = lcalc;
dcomplex arg_2 = arg;
double wn_2 = wn;
double vk_2 = vk;
np_int ndit_2 = ndit;
int isq_2 = isq;
int ibf_2 = ibf;
// for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
if (myompthread == 0) {
sconf_2 = sconf;
gconf_2 = gconf;
c1_2 = c1;
c1ao_2 = c1ao;
c2_2 = c2;
......@@ -415,8 +383,6 @@ void cluster(const string& config_file, const string& data_file, const string& o
}
else {
// this is not thread 0, so do create fresh copies of all local variables
sconf_2 = new ScattererConfiguration(*sconf);
gconf_2 = new GeometryConfiguration(*gconf);
c1_2 = new C1(*c1);
c1ao_2 = new C1_AddOns(*c1ao);
c2_2 = new C2(*c2);
......@@ -545,14 +511,12 @@ void cluster(const string& config_file, const string& data_file, const string& o
// ok, now I can actually start the parallel calculations
#pragma omp for
for (jxi488 = 2; jxi488 <= nxi; jxi488++) {
jer = cluster_jxi488_cycle(jxi488, sconf_2, gconf_2, p_scattering_angles, c1_2, c1ao_2, c2_2, c3_2, c4_2, c6_2, c9_2, output_2, output_path, gaps_2, tqse_2, tqspe_2, tqss_2, tqsps_2, zpv_2, gapm_2, gappm_2, argi_2, args_2, gap_2, gapp_2, tqce_2, tqcpe_2, tqcs_2, tqcps_2, duk_2, tppoan_2, cextlr_2, cext_2, cmullr_2, cmul_2, gapv_2, tqev_2, tqsv_2, nxi_2, nsph_2, mxndm_2, inpol_2, iavm_2, npnt_2, npntts_2, isam_2, lm_2, u_2, us_2, un_2, uns_2, up_2, ups_2, unmp_2, unsmp_2, upmp_2, upsmp_2, scan_2, cfmp_2, sfmp_2, cfsp_2, sfsp_2, sqsfi_2, exri_2, lcalc_2, arg_2, wn_2, vk_2, ndit_2, isq_2, ibf_2, logger);
int jer = cluster_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, c1_2, c1ao_2, c2_2, c3_2, c4_2, c6_2, c9_2, output_2, output_path, gaps_2, tqse_2, tqspe_2, tqss_2, tqsps_2, zpv_2, gapm_2, gappm_2, argi_2, args_2, gap_2, gapp_2, tqce_2, tqcpe_2, tqcs_2, tqcps_2, duk_2, tppoan_2, cextlr_2, cext_2, cmullr_2, cmul_2, gapv_2, tqev_2, tqsv_2, u_2, us_2, un_2, uns_2, up_2, ups_2, unmp_2, unsmp_2, upmp_2, upsmp_2, scan_2, cfmp_2, sfmp_2, cfsp_2, sfsp_2, sqsfi_2, arg_2, wn_2, vk_2, logger);
}
#pragma omp barrier
// only threads different from 0 have to free local copies of variables and close local files
if (myompthread != 0) {
delete sconf_2;
delete gconf_2;
delete c1_2;
delete c1ao_2;
delete c2_2;
......@@ -631,7 +595,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
}
#pragma omp barrier
{
string message = "Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
message = "Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
logger->log(message);
}
} // closes pragma omp parallel
......@@ -642,7 +606,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
for (int ri = 1; ri < ompnumthreads; ri++) {
// Giovanni, please add here in this loop code to reopen the temporary files, reread them and append them respectively to the global output and tppoan, before closing them
string partial_file_name = output_path + "/c_OCLU_" + to_string(ri);
string message = "Copying ASCII output of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
message = "Copying ASCII output of thread " + to_string(ri) + " of " + to_string(ompnumthreads - 1) + "... ";
logger->log(message);
FILE *partial_output = fopen(partial_file_name.c_str(), "r");
char c = fgetc(partial_output);
......@@ -769,18 +733,34 @@ void cluster(const string& config_file, const string& data_file, const string& o
}
int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, int nxi, int nsph, np_int mxndm, int inpol, int iavm, int npnt, int npntts, int isam, int lm, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double sqsfi, double exri, int lcalc, dcomplex arg, double wn, double vk, np_int ndit, int isq, int ibf, Logger *logger)
int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double sqsfi, dcomplex arg, double wn, double vk, Logger *logger)
{
int nxi = sconf->number_of_scales;
logger->log("INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");
int jer = 0;
int lcalc = 0;
int jaw = 1;
int jwtm = (int)gconf->get_param("jwtm");
int li = (int)gconf->get_param("li");
int le = (int)gconf->get_param("le");
int li = gconf->li;
int le = gconf->le;
int lm = 0;
if (le > lm) lm = le;
if (li > lm) lm = li;
int nsph = gconf->number_of_spheres;
np_int mxndm = gconf->mxndm;
int iavm = gconf->iavm;
int inpol = gconf->in_pol;
int npnt = gconf->npnt;
int npntts = gconf->npntts;
int isam = gconf->iavm;
int jwtm = gconf->jwtm;
np_int ndit = 2 * nsph * c4->nlim;
int isq, ibf;
fprintf(output, "========== JXI =%3d ====================\n", jxi488);
double xi = sconf->get_scale(jxi488 - 1);
double exdc = sconf->get_param("exdc");
int idfc = (int)sconf->get_param("idfc");
double exdc = sconf->exdc;
double exri = sqrt(exdc);
int idfc = (int)sconf->idfc;
double vkarg = 0.0;
if (idfc >= 0) {
vk = xi * wn;
......
src/include/Configuration.h
View file @ b07fdabd
......@@ -42,59 +42,114 @@ class GeometryConfiguration {
protected:
//! \brief Number of spherical components.
int number_of_spheres;
int _number_of_spheres;
//! \brief Maximum field expansion order.
int l_max;
int _l_max;
//! \brief Maximum internal field expansion order.
int li;
int _li;
//! \brief Maximum external field expansion order.
int le;
int _le;
//! \brief Maximum dimension of allocated matrix allowance (deprecated).
int mxndm;
np_int _mxndm;
//! \brief QUESTION: definition?
int iavm;
int _iavm;
//! \brief Incident field polarization status (0 - linear, 1 - circular).
int in_pol;
int _in_pol;
//! \brief Number of transition points. QUESTION: correct?
int npnt;
int _npnt;
//! \brief Transition smoothness. QUESTION: correct?
int npntts;
int _npntts;
//! \brief Type of meridional plane definition.
int meridional_type;
int _isam;
//! \brief Transition matrix layer ID. QUESTION: correct?
int jwtm;
int _jwtm;
//! \brief Incident field initial azimuth.
double in_theta_start;
double _in_theta_start;
//! \brief Incident field azimuth step.
double in_theta_step;
double _in_theta_step;
//! \brief Incident field final azimuth.
double in_theta_end;
double _in_theta_end;
//! \brief Scattered field initial azimuth.
double sc_theta_start;
double _sc_theta_start;
//! \brief Scattered field azimuth step.
double sc_theta_step;
double _sc_theta_step;
//! \brief Scattered field final azimuth.
double sc_theta_end;
double _sc_theta_end;
//! \brief Incident field initial elevation.
double in_phi_start;
double _in_phi_start;
//! \brief Incident field elevation step.
double in_phi_step;
double _in_phi_step;
//! \brief Incident field final elevation.
double in_phi_end;
double _in_phi_end;
//! \brief Scattered field initial elevation.
double sc_phi_start;
double _sc_phi_start;
//! \brief Scattered field elevation step.
double sc_phi_step;
double _sc_phi_step;
//! \brief Scattered field final elevation.
double sc_phi_end;
double _sc_phi_end;
//! \brief Vector of spherical components X coordinates.
double *sph_x;
double *_sph_x;
//! \brief Vector of spherical components Y coordinates.
double *sph_y;
double *_sph_y;
//! \brief Vector of spherical components Z coordinates.
double *sph_z;
double *_sph_z;
public:
//! \brief Read only view on number of spherical components.
const int& number_of_spheres = _number_of_spheres;
//! \brief Read only view on maximum field expansion order.
const int& l_max = _l_max;
//! \brief Read only view on maximum internal field expansion order.
const int& li = _li;
//! \brief Read only view on maximum external field expansion order.
const int& le = _le;
//! \brief Read only view on maximum dimension of allocated matrix allowance (deprecated).
const np_int& mxndm = _mxndm;
//! \brief QUESTION: definition?
const int& iavm = _iavm;
//! \brief Read only view on incident field polarization status (0 - linear, 1 - circular).
const int& in_pol = _in_pol;
//! \brief Read only view on number of transition points. QUESTION: correct?
const int& npnt = _npnt;
//! \brief Read only view on transition smoothness. QUESTION: correct?
const int& npntts = _npntts;
//! \brief Read only view on type of meridional plane definition.
const int& isam = _isam;
//! \brief Read only view on transition matrix layer ID. QUESTION: correct?
const int& jwtm = _jwtm;
//! \brief Read only view on incident field initial azimuth.
const double& in_theta_start = _in_theta_start;
//! \brief Read only view on incident field azimuth step.
const double& in_theta_step = _in_theta_step;
//! \brief Read only view on incident field final azimuth.
const double& in_theta_end = _in_theta_end;
//! \brief Read only view on scattered field initial azimuth.
const double& sc_theta_start = _sc_theta_start;
//! \brief Read only view on scattered field azimuth step.
const double& sc_theta_step = _sc_theta_step;
//! \brief Read only view on scattered field final azimuth.
const double& sc_theta_end = _sc_theta_end;
//! \brief Read only view on incident field initial elevation.
const double& in_phi_start = _in_phi_start;
//! \brief Read only view on incident field elevation step.
const double& in_phi_step = _in_phi_step;
//! \brief Read only view on incident field final elevation.
const double& in_phi_end = _in_phi_end;
//! \brief Read only view on scattered field initial elevation.
const double& sc_phi_start = _sc_phi_start;
//! \brief Read only view on scattered field elevation step.
const double& sc_phi_step = _sc_phi_step;
//! \brief Read only view on scattered field final elevation.
const double& sc_phi_end = _sc_phi_end;
/*
//! \brief Read only view on vector of spherical components X coordinates.
const double *sph_x = _sph_x;
//! \brief Read only view on vector of spherical components Y coordinates.
const double *sph_y = _sph_y;
//! \brief Read only view on vector of spherical components Z coordinates.
const double *sph_z = _sph_z;
*/
/*! \brief Build a scattering geometry configuration structure.
*
* \param nsph: `int` Number of spheres to be used in calculation.
......@@ -129,7 +184,7 @@ public:
*/
GeometryConfiguration(
int nsph, int lm, int in_pol, int npnt, int npntts, int meridional_type,
int li, int le, int mxndm, int iavm,
int li, int le, np_int mxndm, int iavm,
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,
......@@ -137,6 +192,7 @@ public:
double sc_ph_start, double sc_ph_step, double sc_ph_end,
int jwtm
);
/*! \brief Build a scattering geometry configuration structure copying it from an existing one.
*
* \param rhs: `GeometryConfiguration` preexisting object to copy from.
......@@ -159,20 +215,6 @@ public:
* configuration data.
*/
static GeometryConfiguration *from_legacy(const std::string& file_name);
/*! \brief Get the value of a parameter by name.
*
* The proper way to access read-only parameters from outside a class is to define
* public methods that return their values. For configuration operations, whose
* optimization is not critical, it is possible to define a single function that
* returns simple scalar values called by name. Access to more complicated data
* structures, on the other hand, require specialized methods which avoid the
* burden of searching the necessary value across the whole array every time.
*
* \param param_name: `string` Name of the parameter to be retrieved.
* \return value: `double` Value of the requested parameter.
*/
double get_param(const std::string& param_name);
/*! \brief Get the X coordinate of a sphere by its index.
*
......@@ -182,7 +224,7 @@ public:
* \param index: `int` Index of the scale to be retrieved.
* \return scale: `double` The X coordinate of the requested sphere.
*/
double get_sph_x(int index) { return sph_x[index]; }
double get_sph_x(int index) { return _sph_x[index]; }
/*! \brief Get the Y coordinate of a sphere by its index.
*
......@@ -192,7 +234,7 @@ public:
* \param index: `int` Index of the scale to be retrieved.
* \return scale: `double` The Y coordinate of the requested sphere.
*/
double get_sph_y(int index) { return sph_y[index]; }
double get_sph_y(int index) { return _sph_y[index]; }
/*! \brief Get the Z coordinate of a sphere by its index.
*
......@@ -202,7 +244,7 @@ public:
* \param index: `int` Index of the scale to be retrieved.
* \return scale: `double` The Z coordinate of the requested sphere.
*/
double get_sph_z(int index) { return sph_z[index]; }
double get_sph_z(int index) { return _sph_z[index]; }
};
/**
......@@ -215,35 +257,35 @@ class ScattererConfiguration {
protected:
//! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].
dcomplex ***dc0_matrix;
dcomplex ***_dc0_matrix;
//! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].
double *radii_of_spheres;
double *_radii_of_spheres;
//! \brief Matrix of fractional transition radii with size [N_SPHERES x LAYERS].
double **rcf;
double **_rcf;
//! \brief Vector of sphere ID numbers, with size [N_SPHERES].
int *iog_vec;
int *_iog_vec;
//! \brief Vector of layer numbers for every sphere, with size [N_SPHERES].
int *nshl_vec;
int *_nshl_vec;
//! \brief Vector of scale parameters, with size [N_SCALES].
double *scale_vec;
double *_scale_vec;
//! \brief Name of the reference variable type (one of XIV, WNS, WLS, PUS, EVS).
std::string reference_variable_name;
std::string _reference_variable_name;
//! \brief Number of spherical components.
int number_of_spheres;
int _number_of_spheres;
//! \brief Number of configurations.
int configurations;
int _configurations;
//! \brief Number of scales to use in calculation.
int number_of_scales;
int _number_of_scales;
//! \brief Type of dielectric functions (<0 at XIP, =0 as function of XI, >0 constants).
int idfc;
int _idfc;
//! \brief External medium dielectric constant. QUESTION: correct?
double exdc;
double _exdc;
//! \brief WP. QUESTION: better definition?
double wp;
double _wp;
//! \brief Peak XI. QUESTION: correct?
double xip;
double _xip;
//! \brief Flag to control whether to add an external layer.
bool use_external_sphere;
bool _use_external_sphere;
/*! \brief Build configuration from a HDF5 binary input file.
*
......@@ -291,6 +333,25 @@ protected:
*/
void write_legacy(const std::string& file_name);
public:
//! \brief Read only view on name of the reference variable type.
const std::string& reference_variable_name = _reference_variable_name;
//! \brief Read only view on number of spherical components.
const int& number_of_spheres = _number_of_spheres;
//! \brief Read only view on number of configurations.
const int& configurations = _configurations;
//! \brief Read only view on number of scales to use in calculation.
const int& number_of_scales = _number_of_scales;
//! \brief Read only view on type of dielectric functions.
const int& idfc = _idfc;
//! \brief Read only view on external medium dielectric constant.
const double& exdc = _exdc;
//! \brief Read only view on WP.
const double& wp = _wp;
//! \brief Read only view on peak XI.
const double& xip = _xip;
//! \brief Read only view on flag to control whether to add an external layer.
const bool& use_external_sphere = _use_external_sphere;
/*! \brief Build a scatterer configuration structure.
*
* Prepare a default configuration structure by allocating the necessary
......@@ -330,29 +391,16 @@ public:
double wp,
double xip
);
/*! \brief Build a scatterer configuration structure copying its contents from a preexisting one.
*
* 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: `complex double ***` Matrix of reference dielectric constants.
* \param has_external: `bool` Flag to set whether to add an external spherical layer.
* \param exdc: `double` External medium dielectric constant.
* \param wp: `double` wp
* \param xip: `double` xip
* \param rhs: `ScattererConfiguration&` Reference to the ScattereConfiguration
* object to be copied.
*/
ScattererConfiguration(const ScattererConfiguration& rhs);
ScattererConfiguration(const ScattererConfiguration& rhs);
/*! \brief Destroy a scatterer configuration instance.
*/
......@@ -399,7 +447,7 @@ public:
* \param k: `int` Index of the current scale.
* \return radius: `dcomplex` The requested dielectric constant.
*/
dcomplex get_dielectric_constant(int i, int j, int k) { return dc0_matrix[i][j][k]; }
dcomplex get_dielectric_constant(int i, int j, int k) { return _dc0_matrix[i][j][k]; }
/*! \brief Get the ID of a configuration from the index of the sphere.
*
......@@ -409,7 +457,7 @@ public:
* \param index: `int` Index of the sphere.
* \return ID: `int` ID of the configuration to be applied.
*/
int get_iog(int index) { return iog_vec[index]; }
int get_iog(int index) { return _iog_vec[index]; }
/*! \brief Get the number of layers for a given configuration.
*
......@@ -419,8 +467,9 @@ public:
* \param index: `int` Index of the configuration.
* \return nl: `int` The number of layers for the given configuration.
*/
int get_nshl(int index) { return nshl_vec[index]; }
int get_nshl(int index) { return _nshl_vec[index]; }
/*
/*! \brief Get the value of a parameter by name.
*
* The proper way to access read-only parameters from outside a class is to define
......@@ -432,8 +481,9 @@ public:
*
* \param param_name: `string` Name of the parameter to be retrieved.
* \return value: `double` Value of the requested parameter.
*/
double get_param(const std::string& param_name);
*/
/*! \brief Get the radius of a sphere by its index.
*
......@@ -443,7 +493,7 @@ public:
* \param index: `int` Index of the ID to be retrieved.
* \return radius: `double` The requested sphere radius.
*/
double get_radius(int index) { return radii_of_spheres[index]; }
double get_radius(int index) { return _radii_of_spheres[index]; }
/*! \brief Get the value of a scale by its index.
*
......@@ -453,18 +503,8 @@ public:
* \param index: `int` Index of the scale to be retrieved.
* \return scale: `double` The desired scale.
*/
double get_rcf(int row, int column) { return rcf[row][column]; }
/*! \brief Get the reference variable name.
*
* This is a specialized function to get the name of the reference variable as a
* string.
*
* \return name: `string` The name of the variable used to calculate wavelength /
* size scaling.
*/
std::string get_reference() { return reference_variable_name; }
double get_rcf(int row, int column) { return _rcf[row][column]; }
/*! \brief Get the value of a scale by its index.
*
* This is a specialized function to access a scale (generally a wavelength),
......@@ -473,7 +513,7 @@ public:
* \param index: `int` Index of the scale to be retrieved.
* \return scale: `double` The desired scale.
*/
double get_scale(int index) { return scale_vec[index]; }
double get_scale(int index) { return _scale_vec[index]; }
/*! \brief Print the contents of the configuration object to terminal.
*
......
src/libnptm/Commons.cpp
View file @ b07fdabd
......@@ -17,13 +17,13 @@
#endif
C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
lm = (int)gconf->get_param("l_max");
int li = (int)gconf->get_param("li");
int le = (int)gconf->get_param("le");
lm = gconf->l_max;
int li = gconf->li;
int le = gconf->le;
if (lm == 0) {
lm = (li > le) ? li : le;
}
nsph = (int)gconf->get_param("nsph");
nsph = gconf->number_of_spheres;
nlmmt = 2 * (lm * (lm + 2));
vec_rmi = new dcomplex[nsph * lm]();
......@@ -37,7 +37,7 @@ C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
vec_w = new dcomplex[4 * nlmmt]();
w = new dcomplex*[nlmmt];
for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);
configurations = (int)sconf->get_param("configurations");
configurations = sconf->configurations;
vint = new dcomplex[16]();
vec_vints = new dcomplex[nsph * 16]();
vints = new dcomplex*[nsph];
......@@ -359,12 +359,12 @@ C1_AddOns::~C1_AddOns() {
}
C2::C2(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
nsph = (int)gconf->get_param("nsph");
int npnt = (int)gconf->get_param("npnt");
int npntts = (int)gconf->get_param("npntts");
nsph = gconf->number_of_spheres;
int npnt = gconf->npnt;
int npntts = gconf->npntts;
int max_n = (npnt > npntts) ? npnt : npntts;
nhspo = 2 * max_n - 1;
nl = (int)sconf->get_param("configurations");
nl = sconf->configurations;
if (nsph == 1 && nl == 1) nl = 5;
ris = new dcomplex[nhspo]();
dlri = new dcomplex[nhspo]();
......@@ -432,11 +432,11 @@ C3::~C3() {
}
C4::C4(GeometryConfiguration *gconf) {
li = (int)gconf->get_param("li");
le = (int)gconf->get_param("le");
li = gconf->li;
le = gconf->le;
lm = (li > le) ? li : le;
nv3j = (lm * (lm + 1) * (2 * lm + 7)) / 6;
nsph = (int)gconf->get_param("nsph");
nsph = gconf->number_of_spheres;
// The following is needed to initialize C1_AddOns
litpo = li + li + 1;
litpos = litpo * litpo;
......@@ -530,19 +530,19 @@ C9::~C9() {
}
ScatteringAngles::ScatteringAngles(GeometryConfiguration *gconf) {
int isam = (int)gconf->get_param("meridional_type");
_th = gconf->get_param("in_theta_start");
_thstp = gconf->get_param("in_theta_step");
_thlst = gconf->get_param("in_theta_end");
_ths = gconf->get_param("sc_theta_start");
_thsstp = gconf->get_param("sc_theta_step");
_thslst = gconf->get_param("sc_theta_end");
_ph = gconf->get_param("in_phi_start");
_phstp = gconf->get_param("in_phi_step");
_phlst = gconf->get_param("in_phi_end");
_phs = gconf->get_param("sc_phi_start");
_phsstp = gconf->get_param("sc_phi_step");
_phslst = gconf->get_param("sc_phi_end");
int isam = gconf->isam;
_th = gconf->in_theta_start;
_thstp = gconf->in_theta_step;
_thlst = gconf->in_theta_end;
_ths = gconf->sc_theta_start;
_thsstp = gconf->sc_theta_step;
_thslst = gconf->sc_theta_end;
_ph = gconf->in_phi_start;
_phstp = gconf->in_phi_step;
_phlst = gconf->in_phi_end;
_phs = gconf->sc_phi_start;
_phsstp = gconf->sc_phi_step;
_phslst = gconf->sc_phi_end;
double small = 1.0e-3;
_nth = 0;
_nph = 0;
......
src/libnptm/Configuration.cpp
View file @ b07fdabd
This diff is collapsed.
src/sphere/sphere.cpp
View file @ b07fdabd
......@@ -74,8 +74,8 @@ void sphere(const string& config_file, const string& data_file, const string& ou
delete logger;
exit(1);
}
int s_nsph = (int)sconf->get_param("nsph");
int nsph = (int)gconf->get_param("nsph");
int s_nsph = sconf->number_of_spheres;
int nsph = gconf->number_of_spheres;
if (s_nsph == nsph) {
int isq, ibf;
double *argi, *args, *gaps;
......@@ -114,25 +114,25 @@ void sphere(const string& config_file, const string& data_file, const string& ou
double frx = 0.0, fry = 0.0, frz = 0.0;
double cfmp, cfsp, sfmp, sfsp;
int jw;
int l_max = gconf->get_param("l_max");
int l_max = gconf->l_max;
C1 *c1 = new C1(gconf, sconf);
int npnt = (int)gconf->get_param("npnt");
int npntts = (int)gconf->get_param("npntts");
int in_pol = (int)gconf->get_param("in_pol");
int meridional_type = (int)gconf->get_param("meridional_type");
int jwtm = (int)gconf->get_param("jwtm");
double in_theta_start = gconf->get_param("in_theta_start");
double in_theta_step = gconf->get_param("in_theta_step");
double in_theta_end = gconf->get_param("in_theta_end");
double sc_theta_start = gconf->get_param("sc_theta_start");
double sc_theta_step = gconf->get_param("sc_theta_step");
double sc_theta_end = gconf->get_param("sc_theta_end");
double in_phi_start = gconf->get_param("in_phi_start");
double in_phi_step = gconf->get_param("in_phi_step");
double in_phi_end = gconf->get_param("in_phi_end");
double sc_phi_start = gconf->get_param("sc_phi_start");
double sc_phi_step = gconf->get_param("sc_phi_step");
double sc_phi_end = gconf->get_param("sc_phi_end");
int npnt = gconf->npnt;
int npntts = gconf->npntts;
int in_pol = gconf->in_pol;
int meridional_type = gconf->iavm;
int jwtm = gconf->jwtm;
double in_theta_start = gconf->in_theta_start;
double in_theta_step = gconf->in_theta_step;
double in_theta_end = gconf->in_theta_end;
double sc_theta_start = gconf->sc_theta_start;
double sc_theta_step = gconf->sc_theta_step;
double sc_theta_end = gconf->sc_theta_end;
double in_phi_start = gconf->in_phi_start;
double in_phi_step = gconf->in_phi_step;
double in_phi_end = gconf->in_phi_end;
double sc_phi_start = gconf->sc_phi_start;
double sc_phi_step = gconf->sc_phi_step;
double sc_phi_end = gconf->sc_phi_end;
C2 *c2 = new C2(gconf, sconf);
argi = new double[1];
args = new double[1];
......@@ -238,7 +238,7 @@ void sphere(const string& config_file, const string& data_file, const string& ou
}
}
thdps(l_max, zpv);
double exdc = sconf->get_param("exdc");
double exdc = sconf->exdc;
double exri = sqrt(exdc);
fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
fstream tppoan;
......@@ -261,13 +261,13 @@ void sphere(const string& config_file, const string& data_file, const string& ou
if (in_pol == 0) fprintf(output, " LIN\n");
else fprintf(output, " CIRC\n");
fprintf(output, " \n");
double wp = sconf->get_param("wp");
double xip = sconf->get_param("xip");
double wp = sconf->wp;
double xip = sconf->xip;
double wn = wp / 3.0e8;
double sqsfi = 1.0;
double vk, vkarg;
int idfc = (int)sconf->get_param("idfc");
int nxi = (int)sconf->get_param("number_of_scales");
int idfc = sconf->idfc;
int nxi = sconf->number_of_scales;
if (idfc < 0) {
vk = xip * wn;
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
......
src/trapping/cfrfme.cpp
View file @ b07fdabd
......@@ -176,13 +176,13 @@ void frfme(string data_file, string output_path) {
ScattererConfiguration *tedf = ScattererConfiguration::from_binary(tedf_name, "HDF5");
if (tedf != NULL) {
int iduml, idum;
iduml = (int)tedf->get_param("nsph");
iduml = tedf->number_of_spheres;
idum = tedf->get_iog(iduml - 1);
exdc = tedf->get_param("exdc");
wp = tedf->get_param("wp");
xip = tedf->get_param("xip");
idfc = (int)tedf->get_param("idfc");
nxi = (int)tedf->get_param("nxi");
exdc = tedf->exdc;
wp = tedf->wp;
xip = tedf->xip;
idfc = tedf->idfc;
nxi = tedf->number_of_scales;
if (idfc >= 0) {
if (ixi <= nxi) {
xi = tedf->get_scale(ixi - 1);
......