Commit 61965407 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Fix dynamic order handling for CLUSTER case

parent d6e28991

src/cluster/cluster.cpp

+11 −15
Original line number Diff line number Diff line
@@ -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;

}

src/libnptm/Commons.cpp

+11 −12
Original line number Diff line number Diff line
@@ -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;

}