diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index f3e9141dd5567e96a0155c03099870a1f6d7fcd8..77a149a4234063dad795003cbf874ccd6f9447ff 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -19,6 +19,7 @@
* \brief Implementation of the calculation for a cluster of spheres.
*/
#include <chrono>
+#include <cmath>
#include <cstdio>
#include <exception>
#include <fstream>
@@ -228,6 +229,10 @@ void cluster(const string& config_file, const string& data_file, const string& o
int nsph = gconf->number_of_spheres;
// Sanity check on number of sphere consistency, should always be verified
if (s_nsph == nsph) {
+ char virtual_line[256];
+ sprintf(virtual_line, "%.5g.\n", sconf->get_particle_radius(gconf));
+ message = "INFO: particle radius is " + (string)virtual_line;
+ logger->log(message);
ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);
double wp = sconf->wp;
// ClusterOutputInfo : Thread 0 of MPI process 0 allocates the memory to
@@ -240,7 +245,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
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");
- str(sconf, cid->c1);
+ str(sconf->_rcf, cid->c1);
thdps(cid->c1->lm, cid->zpv);
double exdc = sconf->exdc;
double exri = sqrt(exdc);
@@ -676,10 +681,14 @@ void cluster(const string& config_file, const string& data_file, const string& o
#endif
}
-int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, ClusterOutputInfo *output, const string& output_path, VirtualBinaryFile *vtppoanp)
-{
+int cluster_jxi488_cycle(
+ int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf,
+ ScatteringAngles *sa, ClusterIterationData *cid, ClusterOutputInfo *output,
+ const string& output_path, VirtualBinaryFile *vtppoanp
+) {
int nxi = sconf->number_of_scales;
const dcomplex cc0 = 0.0 + I * 0.0;
+ const double pi = acos(-1.0);
char virtual_line[256];
string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";
Logger *logger = new Logger(LOG_DEBG);
@@ -689,24 +698,19 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
int jer = 0;
int lcalc = 0;
int jaw = 1;
- int li = gconf->li;
- int le = gconf->le;
- int lm = 0;
- if (le > lm) lm = le;
- if (li > lm) lm = li;
- int nsph = sconf->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->isam;
- int jwtm = gconf->jwtm;
- np_int ndit = 2 * nsph * cid->c1->nlim;
+ bool is_first_scale = (jxi488 == 1);
+ const int nsph = sconf->number_of_spheres;
+ const np_int mxndm = gconf->mxndm;
+ const int iavm = gconf->iavm;
+ const int inpol = gconf->in_pol;
+ const int npnt = gconf->npnt;
+ const int npntts = gconf->npntts;
+ const int isam = gconf->isam;
+ const int jwtm = gconf->jwtm;
int isq, ibf;
int last_configuration;
- int num_configs = sconf->configurations;
- int ndirs = sa->nkks;
+ const int num_configs = sconf->configurations;
+ const int ndirs = sa->nkks;
int oindex = -1;
int jindex = jxi488 - output->first_xi + 1;
@@ -716,7 +720,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
double xi = sconf->get_scale(jxi488 - 1);
double exdc = sconf->exdc;
double exri = sqrt(exdc);
- int idfc = (int)sconf->idfc;
+ const int idfc = (int)sconf->idfc;
double vkarg = 0.0;
if (idfc >= 0) {
cid->vk = xi * cid->wn;
@@ -729,6 +733,49 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
output->vec_vk[jindex - 1] = cid->vk;
output->vec_xi[jindex - 1] = xi;
}
+ // Dynamic order check
+ const int max_li = gconf->li;
+ const int max_le = gconf->le;
+ const double alamb = 2.0 * pi / cid->vk;
+ double size_par_li = 2.0 * pi * sqrt(exdc) * sconf->get_max_radius() / alamb;
+ int recommended_li = 4 + (int)ceil(size_par_li + 4.05 * pow(size_par_li, 1.0 / 3.0));
+ double size_par_le = 2.0 * pi * sqrt(exdc) * sconf->get_particle_radius(gconf) / alamb;
+ int recommended_le = 1 + (int)ceil(size_par_le + 11.0 * pow(size_par_le, 1.0 / 3.0));
+ if (recommended_li != cid->c1->li || recommended_le != cid->c1->le) {
+ if (recommended_li > cid->c1->li) {
+ message = "WARNING: internal order " + to_string(cid->c1->li) + " for scale iteration "
+ + to_string(jxi488) + " too low (recommended order is " + to_string(recommended_li)
+ + ").\n";
+ logger->log(message, LOG_WARN);
+ } else if (recommended_li < cid->c1->li) {
+ message = "INFO: lowering internal order from " + to_string(cid->c1->li) + " to "
+ + to_string(recommended_li) + " for scale iteration " + to_string(jxi488) + ".\n";
+ logger->log(message, LOG_INFO);
+ }
+ if (recommended_le > cid->c1->le) {
+ message = "WARNING: external order " + to_string(cid->c1->le) + " for scale iteration "
+ + to_string(jxi488) + " too low (recommended order is " + to_string(recommended_le)
+ + ").\n";
+ logger->log(message, LOG_WARN);
+ } else if (recommended_le < cid->c1->le) {
+ message = "INFO: lowering external order from " + to_string(cid->c1->le) + " to "
+ + to_string(recommended_le) + " for scale iteration " + to_string(jxi488) + ".\n";
+ logger->log(message, LOG_INFO);
+ }
+ if (recommended_li < max_li || recommended_le < max_le) {
+ int new_li = (recommended_li < max_li) ? recommended_li : max_li;
+ int new_le = (recommended_le < max_le) ? recommended_le : max_le;
+ cid->update_orders(sconf->_rcf, new_li, new_le);
+ is_first_scale = true;
+ jaw = 1;
+ cid->refinemode = 2;
+ }
+ }
+ int li = cid->c1->li;
+ int le = cid->c1->le;
+ int lm = cid->c1->lm;
+ np_int ndit = 2 * nsph * cid->c1->nlim;
+ // End of dynamic order check
hjv(exri, vkarg, jer, lcalc, cid->arg, cid->c1);
if (jer != 0) {
output->vec_ier[jindex - 1] = 1;
@@ -754,7 +801,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
for (int ic = 0; ic < ici; ic++)
cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1);
} else {
- if (jxi488 == 1) {
+ if (is_first_scale) {
for (int ic = 0; ic < ici; ic++)
cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0);
}
@@ -837,6 +884,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
}
}
#endif
+ cid->refinemode = 0;
#ifdef DEBUG_AM
VirtualAsciiFile *outam2 = new VirtualAsciiFile();
string outam2_name = output_path + "/c_AM2_JXI" + to_string(jxi488) + ".txt";
@@ -972,7 +1020,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
for (int jph484 = 1; jph484 <= sa->nph; jph484++) {
int jw = 0;
- if (sa->nk != 1 || jxi488 <= 1) {
+ if (sa->nk != 1 || is_first_scale) {
upvmp(th, ph, 0, cost, sint, cosp, sinp, cid->u, cid->upmp, cid->unmp);
if (isam >= 0) {
wmamp(
@@ -984,7 +1032,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
raba(cid->c1->le, cid->c1->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);
jw = 1;
}
- } else { // label 180, NK == 1 AND JXI488 == 1
+ } else { // label 180, NK == 1 AND JXI488 > 1
if (isam >= 0) {
// label 182
apc(cid->zpv, cid->c1->le, cid->c1->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
@@ -1015,7 +1063,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
if (phs > 360.0) phs -= 360.0;
}
// label 188
- bool goto190 = (sa->nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
+ bool goto190 = (sa->nks == 1 && (!(is_first_scale) || jth486 > 1 || jph484 > 1));
if (!goto190) {
upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, cid->us, cid->upsmp, cid->unsmp);
if (isam >= 0)
@@ -1025,7 +1073,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
);
}
// label 190
- if (sa->nkks != 1 || jxi488 <= 1) {
+ if (sa->nkks != 1 || is_first_scale) {
upvsp(
cid->u, cid->upmp, cid->unmp, cid->us, cid->upsmp, cid->unsmp, cid->up, cid->un, cid->ups, cid->uns,
cid->duk, isq, ibf, cid->scan, cid->cfmp, cid->sfmp, cid->cfsp, cid->sfsp
@@ -2062,4 +2110,43 @@ ClusterIterationData::~ClusterIterationData() {
delete[] cmullr;
delete[] cmul;
}
+
+int ClusterIterationData::update_orders(double **rcf, int inner_order, int outer_order) {
+ int result = 0;
+ int old_lm = c1->lm;
+ np_int old_ndit = 2 * c1->nsph * c1->nlim;
+ ((ParticleDescriptorCluster *)c1)->update_orders(inner_order, outer_order);
+ const int ndi = c1->ndi;
+ const np_int ndit = (np_int)c1->ndit;
+ for (int zi = 0; zi < old_lm; zi++) {
+ for (int zj = 0; zj < 3; zj++) {
+ for (int zk = 0; zk < 2; zk++) {
+ delete[] zpv[zi][zj][zk];
+ }
+ delete[] zpv[zi][zj];
+ }
+ delete[] zpv[zi];
+ }
+ delete[] zpv;
+ zpv = new double***[c1->lm];
+ for (int zi = 0; zi < c1->lm; zi++) {
+ zpv[zi] = new double**[3];
+ for (int zj = 0; zj < 3; zj++) {
+ zpv[zi][zj] = new double*[2];
+ for (int zk = 0; zk < 2; zk++) {
+ zpv[zi][zj][zk] = new double[2]();
+ }
+ }
+ }
+ str(rcf, c1);
+ thdps(c1->lm, zpv);
+ delete[] am;
+ delete[] am_vector;
+ am_vector = new dcomplex[ndit * ndit]();
+ am = new dcomplex*[ndit];
+ for (int ai = 0; ai < ndit; ai++) {
+ am[ai] = am_vector + ai * ndit;
+ }
+ return result;
+}
// >>> END OF ClusterIterationData CLASS IMPLEMENTATION <<<
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 31efcc39feb264b0b62b275076477b9402a5e475..54f062e5d35adc5a808193ed6b1441d23ecac927 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -451,6 +451,14 @@ public:
* \return descriptor_type: `string` The descriptor type name.
*/
std::string get_descriptor_type() override { return "cluster descriptor"; }
+
+ /*! \brief Update the field expansion orders.
+ *
+ * \param inner_order: `int` The new inner expansion order to be set.
+ * \param outer_order: `int` The new outer expansion order to be set.
+ * \return result: `int` An exit code (0 if successful).
+ */
+ int update_orders(int inner_order, int outer_order);
};
/*! \brief The data structure describing a particle model for a sphere with inclusions.
@@ -493,6 +501,14 @@ public:
* \return descriptor_type: `string` The descriptor type name.
*/
std::string get_descriptor_type() override { return "inclusion descriptor"; }
+
+ /*! \brief Update the field expansion orders.
+ *
+ * \param inner_order: `int` The new inner expansion order to be set.
+ * \param outer_order: `int` The new outer expansion order to be set.
+ * \return result: `int` An exit code (0 if successful).
+ */
+ int update_orders(int inner_order, int outer_order);
};
/*! \brief The data structure describing a spherical particle model.
@@ -535,6 +551,13 @@ public:
* \return descriptor_type: `string` The descriptor type name.
*/
std::string get_descriptor_type() override { return "sphere descriptor"; }
+
+ /*! \brief Update the field expansion order.
+ *
+ * \param order: `int` The new field expansion order to be set.
+ * \return result: `int` An exit code (0 if successful).
+ */
+ int update_order(int order);
};
/*! \brief A data structure representing the angles to be evaluated in the problem.
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index d23aa6ff178ee60c8f251ab9b88a86349b6e1762..4235f56d9d42cb4849380282f66b90871b5bf1cb 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -288,11 +288,9 @@ protected:
dcomplex ***_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].
+ //! \brief Vector of layer numbers for every sphere, with size [CONFIGURATIONS].
int *_nshl_vec;
//! \brief Vector of scale parameters, with size [N_SCALES].
double *_scale_vec;
@@ -383,6 +381,8 @@ public:
const int& max_layers = _max_layers;
//! \brief Read only view on flag to control whether to add an external layer.
const bool& use_external_sphere = _use_external_sphere;
+ //! \brief Matrix of fractional transition radii with size [CONFIGURATIONS x LAYERS].
+ double **_rcf;
/*! \brief Build a scatterer configuration structure.
*
@@ -511,6 +511,12 @@ public:
*/
int get_iog(int index) { return _iog_vec[index]; }
+ /*! \brief Get the maximum radius of the sphere components.
+ *
+ * \return radius: `double` The radius of the largest sphere.
+ */
+ double get_max_radius();
+
/*! \brief Get the number of layers for a given configuration.
*
* This is a specialized function to get the number of layers in a specific
@@ -521,6 +527,13 @@ public:
*/
int get_nshl(int index) { return _nshl_vec[index]; }
+ /*! \brief Get the radius of the smallest sphere containing the particle.
+ *
+ * \param gc: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` instance.
+ * \return radius: `double` The radius of the sphere containing the particle.
+ */
+ double get_particle_radius(GeometryConfiguration *gc);
+
/*! \brief Get the radius of a sphere by its index.
*
* This is a specialized function to get the radius of a sphere through its
diff --git a/src/include/IterationData.h b/src/include/IterationData.h
index bf38b033096cfd0329dd53ddf075a2f70f36ec26..4fe1d250ddd847e5a812e0fb1933ec1b367789c2 100644
--- a/src/include/IterationData.h
+++ b/src/include/IterationData.h
@@ -178,6 +178,15 @@ public:
/*! \brief `ClusterIterationData` instance destroyer.
*/
~ClusterIterationData();
+
+ /*! \brief Update field expansion orders.
+ *
+ * \param rcf: `double **` Matrix of sphere fractional radii.
+ * \param inner_order: `int` The new inner expansion order to be set.
+ * \param outer_order: `int` The new outer expansion order to be set.
+ * \return result: `int` An exit code (0 if successful).
+ */
+ int update_orders(double** rcf, int inner_order, int outer_order);
};
// >>> END OF ClusterIterationData CLASS DEFINITION <<<
@@ -345,6 +354,15 @@ public:
/*! \brief `InclusionIterationData` instance destroyer.
*/
~InclusionIterationData();
+
+ /*! \brief Update field expansion orders.
+ *
+ * \param rcf: `double **` Matrix of sphere fractional radii.
+ * \param inner_order: `int` The new inner expansion order to be set.
+ * \param outer_order: `int` The new outer expansion order to be set.
+ * \return result: `int` An exit code (0 if successful).
+ */
+ int update_orders(double** rcf, int inner_order, int outer_order);
};
// >>> END OF InclusionIterationData CLASS DEFINITION <<< //
@@ -503,6 +521,13 @@ public:
/*! \brief `SphereIterationData` instance destroyer.
*/
~SphereIterationData();
+
+ /*! \brief Update field expansion order.
+ *
+ * \param order: `int` The new expansion order to be set.
+ * \return result: `int` An exit code (0 if successful).
+ */
+ int update_order(int order);
};
// >>> END OF SphereIterationData CLASS DEFINITION <<<
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 9c5feffe260086c5b8a42fa7608a5bb82585f35a..40c5f0e322702b1a3e3c7801cdaf400fc792e852 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -369,10 +369,10 @@ void scr2(
* to radial coordinates, then it calls `sphar()` to calculate the vector of spherical
* harmonics of the incident field.
*
- * \param sconf: `ScattererConfiguration *` Pointer to scatterer configuration object.
+ * \param rcf: `double **` Matrix of sphere configuration fractional radii.
* \param c1: `ParticleDescriptor *` Pointer to a ParticleDescriptor instance.
*/
-void str(ScattererConfiguration *sconf, ParticleDescriptor *c1);
+void str(double **rcf, ParticleDescriptor *c1);
/*! \brief Compute radiation torques on particles on a k-vector oriented system.
*
diff --git a/src/include/inclu_subs.h b/src/include/inclu_subs.h
index 64e5a1f3da5d11fc060f709b60ac50d4d2055554..c7af5e26c70594e4a0d9184acf33fffc1de9113e 100644
--- a/src/include/inclu_subs.h
+++ b/src/include/inclu_subs.h
@@ -75,10 +75,10 @@ void indme(
/*! \brief C++ porting of INSTR.
*
- * \param sconf: `ScattererConfiguration *` Pointer to a ScattererConfiguration instance.
+ * \param rcf: `double **` Pointer to matrix of fractional radii.
* \param c1: `ParticleDescriptor *` Pointer to a ParticleDescriptor instance.
*/
-void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1);
+void instr(double **rcf, ParticleDescriptor *c1);
/*! \brief C++ porting of OSPV.
*
diff --git a/src/inclusion/inclusion.cpp b/src/inclusion/inclusion.cpp
index 0b0033c4a20f7609c5fd5e9dc940efdf8f7509df..7d6d4ae38b3403862747b80dc9ac3e310e05d101 100644
--- a/src/inclusion/inclusion.cpp
+++ b/src/inclusion/inclusion.cpp
@@ -19,6 +19,7 @@
* \brief Implementation of the calculation for a sphere with inclusions.
*/
#include <chrono>
+#include <cmath>
#include <cstdio>
#include <exception>
#include <fstream>
@@ -233,12 +234,10 @@ void inclusion(const string& config_file, const string& data_file, const string&
// Open empty virtual ascii file for output
InclusionOutputInfo *p_output = new InclusionOutputInfo(sconf, gconf, mpidata);
InclusionIterationData *cid = new InclusionIterationData(gconf, sconf, mpidata, device_count);
- const np_int ndi = cid->c1->nsph * cid->c1->nlim;
- const np_int ndit = 2 * ndi;
logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)cid->c1->ndm) + " x " + to_string((int64_t)cid->c1->ndm) +".\n");
time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)cid->c1->ndm) + " x " + to_string((int64_t)cid->c1->ndm) +".\n");
- instr(sconf, cid->c1);
+ instr(sconf->_rcf, cid->c1);
thdps(cid->c1->lm, cid->zpv);
double exdc = sconf->exdc;
double exri = sqrt(exdc);
@@ -537,6 +536,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
time_logger->log(message);
fclose(timing_file);
delete time_logger;
+ delete logger;
} // end instructions block of MPI process 0
//===============================
@@ -673,8 +673,9 @@ void inclusion(const string& config_file, const string& data_file, const string&
}
int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, InclusionOutputInfo *output, const string& output_path, VirtualBinaryFile *vtppoanp) {
- const dcomplex cc0 = 0.0 + I * 0.0;
int nxi = sconf->number_of_scales;
+ const dcomplex cc0 = 0.0 + I * 0.0;
+ const double pi = acos(-1.0);
char virtual_line[256];
string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";
Logger *logger = new Logger(LOG_DEBG);
@@ -684,24 +685,24 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
int jer = 0;
int lcalc = 0;
int jaw = 1;
- int li = cid->c1->li;
- int le = cid->c1->le;
- int lm = cid->c1->lm;
- int nsph = cid->c1->nsph;
+ bool is_first_scale = (jxi488 == 1);
+ // int li = cid->c1->li;
+ // int le = cid->c1->le;
+ // int lm = cid->c1->lm;
+ const int nsph = cid->c1->nsph;
np_int mxndm = gconf->mxndm;
- int iavm = gconf->iavm;
- int inpol = gconf->in_pol;
- int npnt = cid->c1->npnt;
- int npntts = cid->c1->npntts;
- int isam = gconf->iavm;
- int jwtm = gconf->jwtm;
- np_int ndit = cid->c1->ndit;
+ const int iavm = gconf->iavm;
+ const int inpol = gconf->in_pol;
+ const int npnt = cid->c1->npnt;
+ const int npntts = cid->c1->npntts;
+ const int isam = gconf->iavm;
+ const int jwtm = gconf->jwtm;
int isq, ibf;
int last_configuration;
dcomplex ent, entn;
double enti;
- int num_configs = sconf->configurations;
- int ndirs = sa->nkks;
+ const int num_configs = sconf->configurations;
+ const int ndirs = sa->nkks;
int oindex = -1;
int jindex = jxi488 - output->first_xi + 1;
@@ -725,6 +726,48 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
output->vec_vk[jindex - 1] = cid->vk;
output->vec_xi[jindex - 1] = xi;
}
+ // Dynamic order check
+ const int max_li = gconf->li;
+ const int max_le = gconf->le;
+ const double alamb = 2.0 * pi / cid->vk;
+ double size_par_li = 2.0 * pi * sqrt(exdc) * sconf->get_max_radius() / alamb;
+ int recommended_li = 4 + (int)ceil(size_par_li + 4.05 * pow(size_par_li, 1.0 / 3.0));
+ double size_par_le = 2.0 * pi * sqrt(exdc) * sconf->get_particle_radius(gconf) / alamb;
+ int recommended_le = 1 + (int)ceil(size_par_le + 11.0 * pow(size_par_le, 1.0 / 3.0));
+ if (recommended_li != cid->c1->li || recommended_le != cid->c1->le) {
+ if (recommended_li > cid->c1->li) {
+ message = "WARNING: internal order " + to_string(cid->c1->li) + " for scale iteration "
+ + to_string(jxi488) + " too low (recommended order is " + to_string(recommended_li)
+ + ").\n";
+ logger->log(message, LOG_WARN);
+ } else if (recommended_li < cid->c1->li) {
+ message = "INFO: lowering internal order from " + to_string(cid->c1->li) + " to "
+ + to_string(recommended_li) + " for scale iteration " + to_string(jxi488) + ".\n";
+ logger->log(message, LOG_INFO);
+ }
+ if (recommended_le > cid->c1->le) {
+ message = "WARNING: external order " + to_string(cid->c1->le) + " for scale iteration "
+ + to_string(jxi488) + " too low (recommended order is " + to_string(recommended_le)
+ + ").\n";
+ logger->log(message, LOG_WARN);
+ } else if (recommended_le < cid->c1->le) {
+ message = "INFO: lowering external order from " + to_string(cid->c1->le) + " to "
+ + to_string(recommended_le) + " for scale iteration " + to_string(jxi488) + ".\n";
+ logger->log(message, LOG_INFO);
+ }
+ if (recommended_li < max_li || recommended_le < max_le) {
+ int new_li = (recommended_li < max_li) ? recommended_li : max_li;
+ int new_le = (recommended_le < max_le) ? recommended_le : max_le;
+ cid->update_orders(sconf->_rcf, new_li, new_le);
+ is_first_scale = true;
+ jaw = 1;
+ cid->refinemode = 2;
+ }
+ }
+ int li = cid->c1->li;
+ int le = cid->c1->le;
+ int lm = cid->c1->lm;
+ // End of dynamic order check
// label 120
double sze = vkarg * cid->extr;
last_configuration = 0;
@@ -745,7 +788,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i133 - 1, jxi488 - 1);
// goes to 129
} else { // label 127
- if (jxi488 == 1) {
+ if (is_first_scale) {
for (int ic = 0; ic < ici; ic++) {
cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i133 - 1, 0);
}
@@ -921,7 +964,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
for (int jph484 = 1; jph484 <= sa->nph; jph484++) {
int jw = 0;
- if (sa->nk != 1 || jxi488 <= 1) {
+ if (sa->nk != 1 || is_first_scale) {
upvmp(th, ph, 0, cost, sint, cosp, sinp, cid->u, cid->upmp, cid->unmp);
if (isam >= 0) {
wmamp(
@@ -964,7 +1007,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
if (phs > 360.0) phs -= 360.0;
}
// label 188
- bool goto190 = (sa->nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
+ bool goto190 = (sa->nks == 1 && (!(is_first_scale) || jth486 > 1 || jph484 > 1));
if (!goto190) {
upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, cid->us, cid->upsmp, cid->unsmp);
if (isam >= 0)
@@ -974,7 +1017,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
);
}
// label 190
- if (sa->nkks != 1 || jxi488 <= 1) {
+ if (sa->nkks != 1 || is_first_scale) {
upvsp(
cid->u, cid->upmp, cid->unmp, cid->us, cid->upsmp, cid->unsmp, cid->up, cid->un, cid->ups, cid->uns,
cid->duk, isq, ibf, cid->scan, cid->cfmp, cid->sfmp, cid->cfsp, cid->sfsp
@@ -1390,8 +1433,6 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
// >>> IMPLEMENTATION OF InclusionIterationData CLASS <<< //
InclusionIterationData::InclusionIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata, const int device_count) {
c1 = new ParticleDescriptorInclusion(gconf, sconf);
- const int ndi = c1->nsph * c1->nlim;
- const np_int ndit = 2 * ndi;
gaps = new double[c1->nsph]();
tqev = new double[3]();
tqsv = new double[3]();
@@ -1502,8 +1543,6 @@ InclusionIterationData::InclusionIterationData(GeometryConfiguration *gconf, Sca
InclusionIterationData::InclusionIterationData(const InclusionIterationData& rhs) {
c1 = new ParticleDescriptorInclusion(reinterpret_cast<ParticleDescriptorInclusion &>(*(rhs.c1)));
- const int ndi = c1->nsph * c1->nlim;
- const np_int ndit = 2 * ndi;
gaps = new double[c1->nsph]();
for (int gi = 0; gi < c1->nsph; gi++) gaps[gi] = rhs.gaps[gi];
tqev = new double[3]();
@@ -1654,8 +1693,6 @@ InclusionIterationData::InclusionIterationData(const InclusionIterationData& rhs
#ifdef MPI_VERSION
InclusionIterationData::InclusionIterationData(const mixMPI *mpidata, const int device_count) {
c1 = new ParticleDescriptorInclusion(mpidata);
- const int ndi = c1->nsph * c1->nlim;
- const np_int ndit = 2 * ndi;
gaps = new double[c1->nsph]();
MPI_Bcast(gaps, c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
tqev = new double[3]();
@@ -1792,8 +1829,6 @@ InclusionIterationData::InclusionIterationData(const mixMPI *mpidata, const int
void InclusionIterationData::mpibcast(const mixMPI *mpidata) {
c1->mpibcast(mpidata);
- const int ndi = c1->nsph * c1->nlim;
- const np_int ndit = 2 * ndi;
MPI_Bcast(gaps, c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(tqev, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Bcast(tqsv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
@@ -1927,4 +1962,40 @@ InclusionIterationData::~InclusionIterationData() {
delete[] cmullr;
delete[] cmul;
}
+
+int InclusionIterationData::update_orders(double **rcf, int inner_order, int outer_order) {
+ int result = 0;
+ int old_lm = c1->lm;
+ ((ParticleDescriptorInclusion *)c1)->update_orders(inner_order, outer_order);
+ for (int zi = 0; zi < old_lm; zi++) {
+ for (int zj = 0; zj < 3; zj++) {
+ for (int zk = 0; zk < 2; zk++) {
+ delete[] zpv[zi][zj][zk];
+ }
+ delete[] zpv[zi][zj];
+ }
+ delete[] zpv[zi];
+ }
+ delete[] zpv;
+ zpv = new double***[c1->lm];
+ for (int zi = 0; zi < c1->lm; zi++) {
+ zpv[zi] = new double**[3];
+ for (int zj = 0; zj < 3; zj++) {
+ zpv[zi][zj] = new double*[2];
+ for (int zk = 0; zk < 2; zk++) {
+ zpv[zi][zj][zk] = new double[2]();
+ }
+ }
+ }
+ instr(rcf, c1);
+ thdps(c1->lm, zpv);
+ delete[] am;
+ delete[] am_vector;
+ am_vector = new dcomplex[c1->ndm * c1->ndm]();
+ am = new dcomplex*[c1->ndm];
+ for (int ai = 0; ai < c1->ndm; ai++) {
+ am[ai] = am_vector + ai * c1->ndm;
+ }
+ return result;
+}
// >>> END OF InclusionIterationData CLASS IMPLEMENTATION <<<
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 29bf18b637b149e7d3d6b9fa71ae8433c7e135cd..b93abab2ce2fa19c341256ed4a9847901eb99f4c 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -83,10 +83,12 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
_ndm = 0;
gcs = 0.0;
_num_configurations = sconf->configurations;
- _num_layers = (sconf->use_external_sphere) ? 1 : 0;
+ // _num_layers = (sconf->use_external_sphere) ? 1 : 0;
+ _num_layers = 0;
_max_layers = 1;
for (int nli = 0; nli < num_configurations; nli++) {
int nl = sconf->get_nshl(nli);
+ if (nli == 0 && sconf->use_external_sphere) nl++;
_num_layers += nl;
if (nl > _max_layers) _max_layers = nl;
}
@@ -136,8 +138,6 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
_npntts = gconf->npntts;
int max_n = (npnt > npntts) ? npnt : npntts;
_nhspo = 2 * max_n - 1;
- // _nl = sconf->configurations;
- // if (_nsph == 1 && _nl == 1) _nl = 5;
ris = new dcomplex[_nhspo]();
dlri = new dcomplex[_nhspo]();
vkt = new dcomplex[_nsph]();
@@ -310,7 +310,6 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) {
tfsas = 0.0 + I * 0.0;
vec_tsas = NULL;
tsas = NULL;
- gcs = 0.0;
scs = 0.0;
ecs = 0.0;
acs = 0.0;
@@ -982,6 +981,91 @@ ParticleDescriptorCluster::ParticleDescriptorCluster(const mixMPI *mpidata) : Pa
MPI_Bcast(rac3j, _lmtpo, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}
#endif // MPI_VERSION
+
+int ParticleDescriptorCluster::update_orders(int inner_order, int outer_order) {
+ int result = 0;
+ bool changed_li = false;
+ bool changed_le = false;
+ if (inner_order != _li) {
+ _li = inner_order;
+ changed_li = true;
+ delete[] vec_rmi;
+ delete[] rmi;
+ vec_rmi = new dcomplex[_li * _nsph]();
+ rmi = new dcomplex*[_li];
+ delete[] vec_rei;
+ delete[] rei;
+ vec_rei = new dcomplex[_li * _nsph]();
+ rei = new dcomplex*[_li];
+ for (int ri = 0; ri < _li; ri++) {
+ rmi[ri] = vec_rmi + (_nsph * ri);
+ rei[ri] = vec_rei + (_nsph * ri);
+ }
+ _litpo = _li + _li + 1;
+ _litpos = _litpo * _litpo;
+ delete[] vh;
+ vh = new dcomplex[_ncou * _litpo]();
+ delete[] vyhj;
+ vyhj = new dcomplex[_ncou * _litpos]();
+ _nlim = _li * (_li + 2);
+ _ndi = _nsph * _nlim;
+ _ndit = 2 * _nsph * _nlim;
+ }
+ if (outer_order != _le) {
+ _le = outer_order;
+ changed_le = true;
+ _nlem = _le * (_le + 2);
+ _nlemt = 2 * _nlem;
+ delete[] vec_am0m;
+ vec_am0m = new dcomplex[_nlemt * _nlemt]();
+ delete[] am0m;
+ am0m = new dcomplex*[_nlemt];
+ for (int ai = 0; ai < _nlemt; ai++) am0m[ai] = vec_am0m + (ai * _nlemt);
+ delete[] vec_w;
+ delete[] w;
+ vec_w = new dcomplex[_nlemt * 4]();
+ w = new dcomplex*[nlemt];
+ for (int wi = 0; wi < nlemt; wi++) w[wi] = vec_w + (4 * wi);
+ }
+ if (changed_li || changed_le) {
+ _lm = (_li > _le) ? _li : _le;
+ _lmpo = _lm + 1;
+ _lmtpo = _li + _le + 1;
+ _lmtpos = _lmtpo * _lmtpo;
+ _nv3j = (_lm * (_lm + 1) * (2 * _lm + 7)) / 6;
+ delete[] vec_ind3j;
+ vec_ind3j = new int[(_lm + 1) * _lm]();
+ delete[] ind3j;
+ ind3j = new int*[_lm + 1];
+ for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
+ delete[] vj0;
+ vj0 = new dcomplex[_nsph * _lmtpo]();
+ delete[] vyj0;
+ vyj0 = new dcomplex[_nsph * _lmtpos]();
+ delete[] v3j0;
+ v3j0 = new double[_nv3j]();
+ delete[] rac3j;
+ rac3j = new double[_lmtpo]();
+ delete[] vec_gis;
+ vec_gis = new dcomplex[_ndi * _nlem]();
+ delete[] gis;
+ gis = new dcomplex*[_ndi];
+ delete[] vec_gls;
+ vec_gls = new dcomplex[_ndi * _nlem]();
+ delete[] gls;
+ gls = new dcomplex*[_ndi];
+ for (int gi = 0; gi < _ndi; gi++) {
+ gis[gi] = vec_gis + (gi * _nlem);
+ gls[gi] = vec_gls + (gi * _nlem);
+ }
+ delete[] vec_sam;
+ vec_sam = new dcomplex[_ndit * _nlemt]();
+ delete[] sam;
+ sam = new dcomplex*[_ndit];
+ for (int si = 0; si < _ndit; si++) sam[si] = vec_sam + (si * _nlemt);
+ }
+ return result;
+}
// >>> End of ParticleDescriptorCluster class implementation. <<< //
// >>> ParticleDescriptorInclusion class implementation. <<< //
@@ -1066,6 +1150,8 @@ ParticleDescriptorInclusion::ParticleDescriptorInclusion(const ParticleDescripto
vec_ind3j = new int[(_lm + 1) * _lm];
np_int vec_ind3j_size = (_lm + 1) * _lm;
for (np_int ii = 0; ii < vec_ind3j_size; ii++) vec_ind3j[ii] = rhs.vec_ind3j[ii];
+ ind3j = new int*[_lm + 1];
+ for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
vh = new dcomplex[_ncou * _litpo];
for (int hi = 0; hi < _ncou * _litpo; hi++) vh[hi] = rhs.vh[hi];
@@ -1106,9 +1192,7 @@ ParticleDescriptorInclusion::ParticleDescriptorInclusion(const ParticleDescripto
ecscm[ci] = rhs.ecscm[ci];
}
v3j0 = new double[_nv3j];
- for (int vj = 0; vj < _nv3j; vj++) v3j0[vj] = rhs.v3j0[_nv3j];
- ind3j = new int*[_lm + 1];
- for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
+ for (int vj = 0; vj < _nv3j; vj++) v3j0[vj] = rhs.v3j0[vj];
rac3j = new double[_lmtpo];
for (int ri = 0; ri < _lmtpo; ri++) rac3j[ri] = rhs.rac3j[ri];
// >>> NEEDED BY INCLU <<< //
@@ -1217,6 +1301,96 @@ ParticleDescriptorInclusion::ParticleDescriptorInclusion(const mixMPI *mpidata)
for (int ai = 0; ai < _nlemt; ai++) at[ai] = vec_at + (ai * _ndm);
}
#endif // MPI_VERSION
+
+int ParticleDescriptorInclusion::update_orders(int inner_order, int outer_order) {
+ int result = 0;
+ bool changed_li = false;
+ bool changed_le = false;
+ if (inner_order != _li) {
+ _li = inner_order;
+ changed_li = true;
+ _nlim = _li * (_li + 2);
+ _litpo = _li + _li + 1;
+ _litpos = _litpo * _litpo;
+ _ndi = _nsph * _nlim;
+ _ndit = 2 * _nsph * _nlim;
+ delete[] rmi;
+ delete[] vec_rmi;
+ delete[] rei;
+ delete[] vec_rei;
+ vec_rmi = new dcomplex[_li * _nsph]();
+ rmi = new dcomplex*[_li];
+ vec_rei = new dcomplex[_li * _nsph]();
+ rei = new dcomplex*[_li];
+ for (int ri = 0; ri < _li; ri++) {
+ rmi[ri] = vec_rmi + (_nsph * ri);
+ rei[ri] = vec_rei + (_nsph * ri);
+ }
+ delete[] vh;
+ vh = new dcomplex[_ncou * _litpo]();
+ delete[] vyhj;
+ vyhj = new dcomplex[_ncou * _litpos]();
+ }
+ if (outer_order != le) {
+ _le = outer_order;
+ changed_le = true;
+ _nlem = _le * (_le + 2);
+ _nlemt = 2 * _nlem;
+ delete[] w;
+ delete[] vec_w;
+ vec_w = new dcomplex[_nlemt * 4]();
+ w = new dcomplex*[_nlemt];
+ for (int wi = 0; wi < _nlemt; wi++) w[wi] = vec_w + (4 * wi);
+ delete[] am0m;
+ delete[] vec_am0m;
+ vec_am0m = new dcomplex[_nlemt * _nlemt]();
+ am0m = new dcomplex*[_nlemt];
+ for (int ai = 0; ai < _nlemt; ai++) am0m[ai] = vec_am0m + (ai * _nlemt);
+ delete[] rm0;
+ rm0 = new dcomplex[_le]();
+ delete[] re0;
+ re0 = new dcomplex[_le]();
+ delete[] rmw;
+ rmw = new dcomplex[_le]();
+ delete[] rew;
+ rew = new dcomplex[_le]();
+ delete[] tm;
+ tm = new dcomplex[_le]();
+ delete[] te;
+ te = new dcomplex[_le]();
+ delete[] tm0;
+ tm0 = new dcomplex[_le]();
+ delete[] te0;
+ te0 = new dcomplex[_le]();
+ }
+ if (changed_li || changed_le) {
+ _lm = (_li > _le) ? _li : _le;
+ _lmpo = _lm + 1;
+ _lmtpo = _li + _le + 1;
+ _lmtpos = _lmtpo * _lmtpo;
+ _nv3j = (_lm * (_lm + 1) * (2 * _lm + 7)) / 6;
+ _ndm = 2 * (_nsph * _nlim + _nlem);
+ delete[] ind3j;
+ delete[] vec_ind3j;
+ vec_ind3j = new int[(_lm + 1) * _lm]();
+ ind3j = new int*[_lm + 1];
+ for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
+ delete[] vj0;
+ vj0 = new dcomplex[_nsph * _lmtpo]();
+ delete[] vyj0;
+ vyj0 = new dcomplex[_nsph * _lmtpos]();
+ delete[] v3j0;
+ v3j0 = new double[_nv3j]();
+ delete[] rac3j;
+ rac3j = new double[_lmtpo]();
+ delete[] at;
+ delete[] vec_at;
+ vec_at = new dcomplex[_nlemt * _ndm]();
+ at = new dcomplex*[_nlemt];
+ for (int ai = 0; ai < _nlemt; ai++) at[ai] = vec_at + (ai * _ndm);
+ }
+ return result;
+}
// >>> End of ParticleDescriptorInclusion class implementation. <<< //
// >>> ParticleDescriptorSphere class implementation. <<< //
@@ -1312,6 +1486,33 @@ ParticleDescriptorSphere::ParticleDescriptorSphere(const mixMPI *mpidata) : Part
}
}
#endif // MPI_VERSION
+
+int ParticleDescriptorSphere::update_order(int order) {
+ int result = 0;
+ if (order != _lm) {
+ _lm = order;
+ _li = order;
+ delete[] rmi;
+ delete[] vec_rmi;
+ delete[] rei;
+ delete[] vec_rei;
+ vec_rmi = new dcomplex[_li * _nsph]();
+ rmi = new dcomplex*[_li];
+ vec_rei = new dcomplex[_li * _nsph]();
+ rei = new dcomplex*[_li];
+ for (int ri = 0; ri < _li; ri++) {
+ rmi[ri] = vec_rmi + (_nsph * ri);
+ rei[ri] = vec_rei + (_nsph * ri);
+ }
+ int nlmmt = 2 * _nsph * _li * (_li + 2);
+ delete[] w;
+ delete[] vec_w;
+ vec_w = new dcomplex[nlmmt * 4]();
+ w = new dcomplex*[nlmmt];
+ for (int wi = 0; wi < nlmmt; wi++) w[wi] = vec_w + (4 * wi);
+ }
+ return result;
+}
// >>> End of ParticleDescriptorSphere class implementation. <<< //
// >>> ScatteringAngles class implementation. <<< //
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 68e076630b26c5eb8bb91b2a1dcd5d8b392392d5..2aa61343aaf91e99294d2d786c72ddae2f290c5b 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -971,6 +971,40 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(const std::string& f
return conf;
}
+double ScattererConfiguration::get_max_radius() {
+ double result = 0.0;
+ for (int ci = 0; ci < _configurations; ci++) {
+ if (_radii_of_spheres[ci] > result)
+ result = _radii_of_spheres[ci];
+ }
+ return result;
+}
+
+double ScattererConfiguration::get_particle_radius(GeometryConfiguration *gc) {
+ double result = 0.0;
+ if (_use_external_sphere) {
+ result = _radii_of_spheres[0] * _rcf[0][_nshl_vec[0] - 1];
+ } else {
+ double x, y, z;
+ int far_index = -1;
+ double dist2, max_dist;
+ double max_dist2 = 0.0;
+ for (int si = 0; si < _number_of_spheres; si++) {
+ x = gc->get_sph_x(si);
+ y = gc->get_sph_y(si);
+ z = gc->get_sph_z(si);
+ dist2 = x * x + y * y + z * z;
+ if (dist2 > max_dist2) {
+ max_dist2 = dist2;
+ far_index = si;
+ }
+ }
+ max_dist = sqrt(max_dist2);
+ result = max_dist + _radii_of_spheres[far_index];
+ }
+ return result;
+}
+
void ScattererConfiguration::print() {
int ies = (_use_external_sphere)? 1 : 0;
printf("### CONFIGURATION DATA ###\n");
diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
index 0199f31a96dd9d759be40f90c07d70cef76c7e0b..3548af3d656cff48a18fddafafb2e61a04542852 100644
--- a/src/libnptm/clu_subs.cpp
+++ b/src/libnptm/clu_subs.cpp
@@ -2442,7 +2442,7 @@ void scr2(
#endif
}
-void str(ScattererConfiguration *sconf, ParticleDescriptor *c1) {
+void str(double **rcf, ParticleDescriptor *c1) {
dcomplex *ylm;
const double pi = acos(-1.0);
int last_configuration;
@@ -2457,7 +2457,7 @@ void str(ScattererConfiguration *sconf, ParticleDescriptor *c1) {
c1->gcsv[i18 - 1] = gcss;
int nsh = c1->nshl[last_configuration - 1];
for (int j16 = 1; j16 <= nsh; j16++) {
- c1->rc[last_configuration - 1][j16 - 1] = sconf->get_rcf(last_configuration - 1, j16 - 1) * c1->ros[last_configuration - 1];
+ c1->rc[last_configuration - 1][j16 - 1] = rcf[last_configuration - 1][j16 - 1] * c1->ros[last_configuration - 1];
} // j16 loop
}
c1->gcs += gcss;
diff --git a/src/libnptm/inclu_subs.cpp b/src/libnptm/inclu_subs.cpp
index e2eede6927f1701ca34656ba8cd1a6051baeecdc..12859ab2a1953eb2387fa5d3c9eb804fbb8c4f0b 100644
--- a/src/libnptm/inclu_subs.cpp
+++ b/src/libnptm/inclu_subs.cpp
@@ -493,7 +493,7 @@ void indme(
delete[] dref;
}
-void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1) {
+void instr(double **rcf, ParticleDescriptor *c1) {
const int ylm_size = (c1->litpos > c1->lmtpos) ? c1->litpos : c1->lmtpos;
dcomplex *ylm = new dcomplex[ylm_size]();
double rx, ry, rz, rr, crth, srth, crph, srph;
@@ -503,7 +503,7 @@ void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1) {
if (c1->iog[i18] >= i) {
int nsh = c1->nshl[i18];
for (int j = 0; j < nsh; j++)
- c1->rc[i18][j] = sconf->get_rcf(i18, j) * c1->ros[i18];
+ c1->rc[i18][j] = rcf[i18][j] * c1->ros[i18];
}
} // i18 loop
int i = 0;
diff --git a/src/libnptm/outputs.cpp b/src/libnptm/outputs.cpp
index 48791dc9a947ef24f9183d974b5855d7ebe24602..8c4ecf7df21e3ed3d04320eeac8ec58666d2c510 100644
--- a/src/libnptm/outputs.cpp
+++ b/src/libnptm/outputs.cpp
@@ -4281,22 +4281,22 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
if (isam >= 0) {
fprintf(
p_outfile, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n",
- vec_dir_un[done_dirs],
- vec_dir_un[done_dirs + 1],
- vec_dir_un[done_dirs + 2]
+ vec_dir_un[3 * done_dirs],
+ vec_dir_un[3 * done_dirs + 1],
+ vec_dir_un[3 * done_dirs + 2]
);
fprintf(
p_outfile, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n",
- vec_dir_uns[done_dirs],
- vec_dir_uns[done_dirs + 1],
- vec_dir_uns[done_dirs + 2]
+ vec_dir_uns[3 * done_dirs],
+ vec_dir_uns[3 * done_dirs + 1],
+ vec_dir_uns[3 * done_dirs + 2]
);
} else { // label 214
fprintf(
p_outfile, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n",
- vec_dir_un[done_dirs],
- vec_dir_un[done_dirs + 1],
- vec_dir_un[done_dirs + 2]
+ vec_dir_un[3 * done_dirs],
+ vec_dir_un[3 * done_dirs + 1],
+ vec_dir_un[3 * done_dirs + 2]
);
}
fprintf(p_outfile, " SINGLE SCATTERER\n");
@@ -5644,7 +5644,7 @@ int SphereOutputInfo::write_legacy(const std::string &file_name) {
int done_dirs = 0;
fprintf(p_outfile, "========== JXI =%3d ====================\n", jxi + 1);
if (idfc >= 0) {
- fprintf(p_outfile, " VK=%15.7lE, XI=%15.7lE\n", vec_xi[jxi], vec_vk[jxi]);
+ fprintf(p_outfile, " VK=%15.7lE, XI=%15.7lE\n", vec_vk[jxi], vec_xi[jxi]);
} else { // IDFC < 0
fprintf(p_outfile, " XI=%15.7lE\n", vec_xi[jxi]);
}
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 288c259f4aee6b48efe3b6980af1e87f01e6dbd8..d123b0ff3927e89da061c1c946190575d19afa9c 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -542,11 +542,13 @@ int sphere_jxi488_cycle(
const double pi = 2.0 * half_pi;
int jer = 0;
Logger *logger = new Logger(LOG_INFO);
+ string message = "INIT";
int oindex = 0;
int jxi = jxi488 + 1;
int jxindex = jxi - oi->first_xi;
+ bool is_first_scale = (jxi == 1);
int nsph = gconf->number_of_spheres;
- int l_max = gconf->l_max;
+ //int l_max = gconf->l_max;
int in_pol = gconf->in_pol;
int npnt = gconf->npnt;
int npntts = gconf->npntts;
@@ -576,6 +578,30 @@ int sphere_jxi488_cycle(
oi->vec_vk[jxindex] = vk;
oi->vec_xi[jxindex] = xi;
}
+ // Dynamic order check
+ const int max_lm = gconf->l_max;
+ int l_max = gconf->l_max;
+ const double alamb = 2.0 * pi / vk;
+ double size_par_lm = 2.0 * pi * sqrt(exdc) * sconf->get_max_radius() / alamb;
+ int recommended_lm = 4 + (int)ceil(size_par_lm + 4.05 * pow(size_par_lm, 1.0 / 3.0));
+ if (recommended_lm != l_max) {
+ if (recommended_lm > max_lm) {
+ message = "WARNING: internal order " + to_string(max_lm) + " for scale iteration "
+ + to_string(jxi488) + " too low (recommended order is " + to_string(recommended_lm)
+ + ").\n";
+ logger->log(message, LOG_WARN);
+ } else if (recommended_lm < max_lm) {
+ int new_lm = recommended_lm;
+ message = "INFO: lowering internal order from " + to_string(max_lm) + " to "
+ + to_string(recommended_lm) + " for scale iteration " + to_string(jxi488) + ".\n";
+ logger->log(message, LOG_INFO);
+ sid->update_order(new_lm);
+ is_first_scale = true;
+ // jw = 1;
+ l_max = new_lm;
+ }
+ }
+ // End of dynamic order check
vtppoanp->append_line(VirtualBinaryLine(vk));
double thsca = (gconf->isam > 1) ? sa->ths - sa->th : 0.0;
for (int i132 = 0; i132 < nsph; i132++) {
@@ -595,7 +621,7 @@ int sphere_jxi488_cycle(
for (int ic = 0; ic < ici; ic++)
sid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132, jxi488); // WARNING: IDFC=0 is not tested!
} else { // IDFC != 0
- if (jxi == 1) {
+ if (is_first_scale) {
for (int ic = 0; ic < ici; ic++) {
sid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i132, jxi488);
}
@@ -718,7 +744,7 @@ int sphere_jxi488_cycle(
double ph = sa->ph;
for (int jph484 = 0; jph484 < sa->nph; jph484++) {
int jph = jph484 + 1;
- bool goto182 = (sa->nk == 1) && (jxi > 1);
+ bool goto182 = (sa->nk == 1) && (!is_first_scale);
if (!goto182) {
upvmp(th, ph, 0, sid->cost, sid->sint, sid->cosp, sid->sinp, sid->u, sid->upmp, sid->unmp);
}
@@ -753,14 +779,14 @@ int sphere_jxi488_cycle(
phs = ph + phsph;
if (phs >= 360.0) phs -= 360.0;
}
- bool goto190 = (nks == 1) && ((jxi > 1) || (jth > 1) || (jph > 1));
+ bool goto190 = (nks == 1) && ((!is_first_scale) || (jth > 1) || (jph > 1));
if (!goto190) {
upvmp(ths, phs, icspnv, sid->costs, sid->sints, sid->cosps, sid->sinps, sid->us, sid->upsmp, sid->unsmp);
if (gconf->isam >= 0) {
wmamp(2, sid->costs, sid->sints, sid->cosps, sid->sinps, in_pol, l_max, 0, nsph, sid->args, sid->us, sid->upsmp, sid->unsmp, sid->c1);
}
}
- if (nkks != 0 || jxi == 1) {
+ if (nkks != 0 || is_first_scale) {
upvsp(
sid->u, sid->upmp, sid->unmp, sid->us, sid->upsmp, sid->unsmp,
sid->up, sid->un, sid->ups, sid->uns, sid->duk, isq, ibf,
@@ -1247,4 +1273,34 @@ int SphereIterationData::mpibcast(const mixMPI *mpidata) {
return 0;
}
#endif // MPI_VERSION
+
+int SphereIterationData::update_order(int order) {
+ int result = 0;
+ int old_lm = _lm;
+ if (order != _lm) {
+ _lm = order;
+ ((ParticleDescriptorSphere *)c1)->update_order(order);
+ for (int zi = 0; zi < old_lm; zi++) {
+ for (int zj = 0; zj < 3; zj++) {
+ delete[] zpv[zi][zj];
+ }
+ delete[] zpv[zi];
+ }
+ delete[] zpv;
+ delete[] vec_zpv;
+ vec_zpv = new double[_lm * 12]();
+ zpv = new double***[_lm];
+ for (int zi = 0; zi < _lm; zi++) {
+ zpv[zi] = new double**[3];
+ for (int zj = 0; zj < 3; zj++) {
+ int vec_index = 12 * zi + 4 * zj;
+ zpv[zi][zj] = new double*[2];
+ zpv[zi][zj][0] = (vec_zpv + vec_index);
+ zpv[zi][zj][1] = (vec_zpv + vec_index + 2);
+ }
+ }
+ thdps(_lm, zpv);
+ }
+ return result;
+}
// >>> END OF SphereIterationData CLASS IMPLEMENTATION <<<