diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index c2fc61cd85b5acb7c5f79c93d094187341d74a99..77a149a4234063dad795003cbf874ccd6f9447ff 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -707,10 +707,9 @@ int cluster_jxi488_cycle(
const int npntts = gconf->npntts;
const int isam = gconf->isam;
const int jwtm = gconf->jwtm;
- // np_int ndit = 2 * nsph * cid->c1->nlim;
int isq, ibf;
int last_configuration;
- int num_configs = sconf->configurations;
+ const int num_configs = sconf->configurations;
const int ndirs = sa->nkks;
int oindex = -1;
int jindex = jxi488 - output->first_xi + 1;
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 24a345525413c8776004f190329fad7286e195ce..54f062e5d35adc5a808193ed6b1441d23ecac927 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -502,13 +502,13 @@ public:
*/
std::string get_descriptor_type() override { return "inclusion 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.
* \return result: `int` An exit code (0 if successful).
*/
- int update_orders(int inner_order, int outer_order) { return 0; }
+ int update_orders(int inner_order, int outer_order);
};
/*! \brief The data structure describing a spherical particle model.
diff --git a/src/include/IterationData.h b/src/include/IterationData.h
index 981ed6518e4d5e31003d2c6c7ebff1014433c68e..4fe1d250ddd847e5a812e0fb1933ec1b367789c2 100644
--- a/src/include/IterationData.h
+++ b/src/include/IterationData.h
@@ -354,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 <<< //
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 e66b7fe21e93207b17362a9671894ec280dee2b1..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;
@@ -1151,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];
@@ -1191,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 <<< //
@@ -1302,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. <<< //
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 216d431d1264ed7add96226c795b6a32f651f38c..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");