From 61965407e6d03847c0b016ba687e97a1504ba5fc Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Tue, 20 May 2025 17:40:41 +0200
Subject: [PATCH] Fix dynamic order handling for CLUSTER case
---
src/cluster/cluster.cpp | 26 +++++++++++---------------
src/libnptm/Commons.cpp | 23 +++++++++++------------
2 files changed, 22 insertions(+), 27 deletions(-)
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index ac13261..de8a046 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -695,11 +695,7 @@ 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;
+ bool is_first_scale = (jxi488 == 1);
const int nsph = sconf->number_of_spheres;
const np_int mxndm = gconf->mxndm;
const int iavm = gconf->iavm;
@@ -766,6 +762,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
int new_li = (recommended_li < cid->c1->li) ? recommended_li : cid->c1->li;
int new_le = (recommended_le < cid->c1->le) ? recommended_le : cid->c1->le;
cid->update_orders(sconf->_rcf, new_li, new_le);
+ is_first_scale = true;
+ jaw = 1;
}
cid->refinemode = 2;
}
@@ -799,7 +797,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);
}
@@ -1018,7 +1016,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(
@@ -1030,7 +1028,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);
@@ -1061,7 +1059,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)
@@ -1071,7 +1069,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
@@ -2111,13 +2109,11 @@ ClusterIterationData::~ClusterIterationData() {
int ClusterIterationData::update_orders(double **rcf, int inner_order, int outer_order) {
int result = 0;
- int old_li = c1->li;
- int old_le = c1->le;
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->nsph * c1->nlim;
- const np_int ndit = 2 * ndi;
+ 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++) {
@@ -2145,7 +2141,7 @@ int ClusterIterationData::update_orders(double **rcf, int inner_order, int outer
am_vector = new dcomplex[ndit * ndit]();
am = new dcomplex*[ndit];
for (int ai = 0; ai < ndit; ai++) {
- am[ai] = (am_vector + ai * ndit);
+ am[ai] = am_vector + ai * ndit;
}
return result;
}
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index dddc30f..61e04e4 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -1008,6 +1008,9 @@ int ParticleDescriptorCluster::update_orders(int inner_order, int outer_order) {
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;
@@ -1019,6 +1022,11 @@ int ParticleDescriptorCluster::update_orders(int inner_order, int outer_order) {
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;
@@ -1026,20 +1034,17 @@ int ParticleDescriptorCluster::update_orders(int inner_order, int outer_order) {
_lmtpo = _li + _le + 1;
_lmtpos = _lmtpo * _lmtpo;
_nv3j = (_lm * (_lm + 1) * (2 * _lm + 7)) / 6;
- _nlim = _li * (_le + 2);
- _ndi = _nsph * _nlim;
- _ndit = 2 * _nsph * _nlim;
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[] ind3j;
- ind3j = new int*[_lm + 1];
- for (int ii = 0; ii <= _lm; ii++) ind3j[ii] = vec_ind3j + (_lm * ii);
delete[] rac3j;
rac3j = new double[_lmtpo]();
delete[] vec_gis;
@@ -1059,12 +1064,6 @@ int ParticleDescriptorCluster::update_orders(int inner_order, int outer_order) {
delete[] sam;
sam = new dcomplex*[_ndit];
for (int si = 0; si < _ndit; si++) sam[si] = vec_sam + (si * _nlemt);
- int nllt = (_nlemt == 0) ? 2 * _nsph * _li * (_li + 2) : _nlemt;
- delete[] vec_w;
- delete[] w;
- vec_w = new dcomplex[nllt * 4]();
- w = new dcomplex*[nllt];
- for (int wi = 0; wi < nllt; wi++) w[wi] = vec_w + (4 * wi);
}
return result;
}
--
GitLab