diff --git a/src/include/Commons.h b/src/include/Commons.h
index 260228fcb2469d010e3a22899822e06a2cc0694d..24a345525413c8776004f190329fad7286e195ce 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -452,7 +452,7 @@ public:
*/
std::string get_descriptor_type() override { return "cluster descriptor"; }
- /*! \brief Interface function to update field expansion orders.
+ /*! \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.
@@ -552,12 +552,12 @@ public:
*/
std::string get_descriptor_type() override { return "sphere descriptor"; }
- /*! \brief Interface function to update field expansion orders.
+ /*! \brief Update the field expansion order.
*
- * \param order: `int` The new expansion order to be set.
+ * \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) { return 0; }
+ int update_order(int order);
};
/*! \brief A data structure representing the angles to be evaluated in the problem.
diff --git a/src/include/IterationData.h b/src/include/IterationData.h
index 1a42d0aa610cb73697e088c34674c57ea7af65c1..981ed6518e4d5e31003d2c6c7ebff1014433c68e 100644
--- a/src/include/IterationData.h
+++ b/src/include/IterationData.h
@@ -512,6 +512,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/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 61e04e4358882072ee4e60f2f3f7a474fbd1f139..e66b7fe21e93207b17362a9671894ec280dee2b1 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -1397,6 +1397,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/sphere/sphere.cpp b/src/sphere/sphere.cpp
index d9a1fd3449ae4fe4153619f63301f2ee5f43384c..b38d2c27a847324edce1de0c263eac56ee0a021b 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,12 +578,30 @@ int sphere_jxi488_cycle(
oi->vec_vk[jxindex] = vk;
oi->vec_xi[jxindex] = xi;
}
- // Adaptive definition of L_MAX
- double wavelength = 2.0 * pi / vk;
- double size_param = 2.0 * pi * sconf->get_radius(0) / wavelength;
- int N = int(size_param + 4.05 * pow(size_param, 1.0 / 3.0)) + 2;
- if (N < l_max) l_max = N;
- // End of adaptive definition of L_MAX
+ // 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++) {
@@ -601,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);
}
@@ -766,7 +786,7 @@ int sphere_jxi488_cycle(
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,
@@ -1253,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 <<<