diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index eae84671732204af8319b4d06b362643121458b1..db5836a479602b39c686fe28f14f82d25e0b294d 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -337,9 +337,32 @@ void cluster(string config_file, string data_file, string output_path) {
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
fprintf(output, " \n");
}
- for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {
+ // do the first iteration on jxi488 separately, since it seems to be different from the others
+ int jxi488 = 1;
+ jer = cluster_jxi488_cycle(jxi488, sconf, gconf, c1, c1ao, c2, c3, c4, c6, c9, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, nth, nths, nph, nphs, nk, nks, nkks, 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, th, thstp, thlst, ths, thsstp, thslst, ph, phstp, phlst, phs, phsstp, phslst, th1, ph1, ths1, phs1, thsca, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, lcalc, arg, wn, vk, ndit, am, isq, ibf);
+ for (jxi488 = 2; jxi488 <= nxi; jxi488++) {
// 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
- jer = cluster_jxi488_cycle(jxi488, sconf, gconf, c1, c1ao, c2, c3, c4, c6, c9, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, nth, nths, nph, nphs, nk, nks, nkks, 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, th, thstp, thlst, ths, thsstp, thslst, ph, phstp, phlst, phs, phsstp, phslst, th1, ph1, ths1, phs1, thsca, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, lcalc, arg, wn, vk, ndit, am, isq, ibf);
+ ScattererConfiguration *sconf_2 = new ScattererConfiguration(*sconf);
+ GeometryConfiguration *gconf_2 = new GeometryConfiguration(*gconf);
+ C1 *c1_2 = new C1(*c1);
+ C1_AddOns *c1ao_2 = new C1_AddOns(*c1ao);
+ C2 *c2_2 = new C2(*c2);
+ C3 *c3_2 = new C3(*c3);
+ C4 *c4_2 = new C4(*c4);
+ C6 *c6_2 = new C6(*c6);
+ C9 *c9_2 = new C9(*c9);
+
+ jer = cluster_jxi488_cycle(jxi488, sconf_2, gconf_2, c1_2, c1ao_2, c2_2, c3_2, c4_2, c6_2, c9_2, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, nth, nths, nph, nphs, nk, nks, nkks, 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, th, thstp, thlst, ths, thsstp, thslst, ph, phstp, phlst, phs, phsstp, phslst, th1, ph1, ths1, phs1, thsca, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, lcalc, arg, wn, vk, ndit, am, isq, ibf);
+
+ delete sconf_2;
+ delete gconf_2;
+ delete c1_2;
+ delete c1ao_2;
+ delete c2_2;
+ delete c3_2;
+ delete c4_2;
+ delete c6_2;
+ delete c9_2;
} // jxi488 loop
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index 7df4b61534af9a87a1da9cace318160d385b88e0..59bd0a1c60eb5b9e41eff8a667e76c62d7300618 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -152,6 +152,11 @@ 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.
+ */
+ GeometryConfiguration(const GeometryConfiguration& rhs);
/*! \brief Destroy a GeometryConfiguration instance.
*/
@@ -199,9 +204,11 @@ protected:
std::string reference_variable_name;
//! \brief Number of spherical components.
int number_of_spheres;
+ //! \brief Number of configurations.
+ int configurations;
//! \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).
+ //! \brief Type of dielectric functions (<0 at XIP, =0 as function of XI, >0 constants).
int idfc;
//! \brief External medium dielectric constant. QUESTION: correct?
double exdc;
@@ -282,6 +289,7 @@ public:
*/
ScattererConfiguration(
int nsph,
+ int configs,
double *scale_vector,
int nxi,
std::string variable_name,
@@ -296,6 +304,29 @@ 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
+ */
+ ScattererConfiguration(const ScattererConfiguration& rhs);
/*! \brief Destroy a scatterer configuration instance.
*/
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 24648680f66e17f88ae5297782e34907128bcc35..fbfb01f70ad4542d1e58c6d5c94d8a4dd7d3ba53 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -33,7 +33,7 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
}
w = new dcomplex*[nlmmt];
for (int wi = 0; wi < nlmmt; wi++) w[wi] = new dcomplex[4]();
- int configurations = 0;
+ configurations = 0;
for (int ci = 1; ci <= nsph; ci++) {
if (_iog[ci - 1] >= ci) configurations++;
}
@@ -149,11 +149,8 @@ C1::~C1() {
delete[] rei;
for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
int conf_index = 0;
- for (int ci = 1; ci <= nsph; ci++) {
- if (iog[ci] >= ci) {
- delete[] rc[conf_index];
- conf_index++;
- }
+ for (int ci = 0; ci < configurations; ci++) {
+ delete[] rc[ci];
}
delete[] rc;
for (int vi = nsph - 1; vi > - 1; vi--) {
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 1bf78b7db3e571858bee6ad893f0c575dbb33644..9036ec3ded6a1beff018f99ef3c03c6a4df7f4de 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -76,6 +76,40 @@ GeometryConfiguration::GeometryConfiguration(
sph_y = y;
sph_z = z;
}
+GeometryConfiguration::GeometryConfiguration(const GeometryConfiguration& rhs)
+{
+ number_of_spheres = rhs.number_of_spheres;
+ l_max = rhs.l_max;
+ in_pol = rhs.in_pol;
+ npnt = rhs.npnt;
+ npntts = rhs.npntts;
+ meridional_type = rhs.meridional_type;
+ li = rhs.li;
+ le = rhs.le;
+ mxndm = rhs.mxndm;
+ iavm = rhs.iavm;
+ in_theta_start = rhs.in_theta_start;
+ in_theta_step = rhs.in_theta_step;
+ in_theta_end = rhs.in_theta_end;
+ in_phi_start = rhs.in_phi_start;
+ in_phi_step = rhs.in_phi_step ;
+ in_phi_end = rhs.in_phi_end;
+ sc_theta_start = rhs.sc_theta_start;
+ sc_theta_step = rhs.sc_theta_step;
+ sc_theta_end = rhs.sc_theta_end;
+ sc_phi_start = rhs.sc_phi_start;
+ sc_phi_step = rhs.sc_phi_step;
+ sc_phi_end = rhs.sc_phi_end;
+ jwtm = rhs.jwtm;
+ sph_x = new double[number_of_spheres]();
+ sph_y = new double[number_of_spheres]();
+ sph_z = new double[number_of_spheres]();
+ for (int ni=0; ni<number_of_spheres; ni++) {
+ sph_x[ni] = rhs.sph_x[ni];
+ sph_y[ni] = rhs.sph_y[ni];
+ sph_z[ni] = rhs.sph_z[ni];
+ }
+}
GeometryConfiguration::~GeometryConfiguration() {
delete[] sph_x;
@@ -199,6 +233,7 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(string file_name) {
ScattererConfiguration::ScattererConfiguration(
int nsph,
+ int configs,
double *scale_vector,
int nxi,
string variable_name,
@@ -214,6 +249,7 @@ ScattererConfiguration::ScattererConfiguration(
double x
) {
number_of_spheres = nsph;
+ configurations = configs;
number_of_scales = nxi;
reference_variable_name = variable_name;
iog_vec = iog_vector;
@@ -241,6 +277,40 @@ ScattererConfiguration::ScattererConfiguration(
}
}
}
+ScattererConfiguration::ScattererConfiguration(const ScattererConfiguration& rhs)
+{
+ number_of_spheres = rhs.number_of_spheres;
+ configurations = rhs.configurations;
+ number_of_scales = rhs.number_of_scales;
+ reference_variable_name = rhs.reference_variable_name;
+ idfc = rhs.idfc;
+ use_external_sphere = rhs.use_external_sphere;
+ exdc = rhs.exdc;
+ wp = rhs.wp;
+ xip = rhs.xip;
+ iog_vec = new int[number_of_spheres]();
+ radii_of_spheres = new double[configurations]();
+ nshl_vec = new int[configurations]();
+ rcf = new double*[configurations];
+ scale_vec = new double[number_of_scales]();
+ dc0_matrix = new dcomplex**[configurations];
+ for (int si = 0; si < number_of_scales; si++) scale_vec[si] = rhs.scale_vec[si];
+ for (int si=0; si<number_of_spheres; si++) {
+ iog_vec[si] = rhs.iog_vec[si];
+ }
+ int dim3 = (idfc == 0) ? number_of_scales : 1;
+ for (int si=0; si<configurations; si++) {
+ radii_of_spheres[si] = rhs.radii_of_spheres[si];
+ nshl_vec[si] = rhs.nshl_vec[si];
+ rcf[si] = new double[nshl_vec[si]]();
+ dc0_matrix[si] = new dcomplex*[number_of_spheres];
+ for (int sj=0; sj<nshl_vec[si]; sj++) rcf[si][sj] = rhs.rcf[si][sj];
+ for (int sj=0; sj<number_of_spheres; sj++) {
+ dc0_matrix[si][sj] = new dcomplex[dim3]();
+ for (int sk=0; sk<dim3; sk++) dc0_matrix[si][sj][sk] = rhs.dc0_matrix[si][sj][sk];
+ }
+ }
+}
ScattererConfiguration::~ScattererConfiguration() {
int configurations = 0;
@@ -519,6 +589,7 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
ScattererConfiguration *config = new ScattererConfiguration(
nsph,
+ configurations,
variable_vector,
nxi,
variable_name,
@@ -614,6 +685,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
status = hdf_file->close();
conf = new ScattererConfiguration(
nsph,
+ configuration_count,
xi_vec,
nxi,
"XIV",
@@ -702,6 +774,7 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
ScattererConfiguration *conf = new ScattererConfiguration(
_nsph,
+ configurations,
_xi_vec,
_nxi,
"XIV",