Commit 82df0eee authored by Giacomo Mulas's avatar Giacomo Mulas
Browse files

Merge branch 'layered_cluster' into 'master' 

Layered cluster

See merge request giacomo.mulas/np_tmcode!79
parents 109e56c4 6a927a23

.gitlab-ci.yml

+18 −7
Original line number Diff line number Diff line
@@ -236,19 +236,27 @@ running_stage: 
      - cd build/libnptm

      - export LD_LIBRARY_PATH="$LD_LIBRARY_PATH;$PWD"

      - cd ../sphere

      - echo "Testing configurator for sphere"

      - echo "Testing configurator for SPHERE"

      - ../../src/scripts/model_maker.py ../../test_data/sphere/config_dev.yml

      - echo "Running np_sphere"

      - chmod +x np_sphere

      - OMP_NUM_THREADS=1 ./np_sphere DEDFB DSPH .

      - mv c_OSPH c_OSPH_dev

      - mv c_TEDF c_TEDF_dev

      - mv c_TEDF.hd5 c_TEDF_dev.hd5

      - echo "Testing parallel execution of SPHERE"

      - ../../src/scripts/model_maker.py ../../test_data/sphere/config_181_xi.yml

      - echo "Running parallel np_sphere"

      - OMP_NUM_THREADS=2 ./np_sphere DEDFB_181_xi DSPH_181_xi .

      - mv c_OSPH c_OSPH_181_xi

      - cd ../cluster

      - echo "Testing configurator for cluster"

      - echo "Testing configurator for CLUSTER"

      - ../../src/scripts/model_maker.py ../../test_data/cluster/config_dev.yml

      - echo "Running np_cluster"

      - chmod +x np_cluster

      - OMP_NUM_THREADS=1 ./np_cluster DEDFB DCLU .

      - cd ../inclusion

      - echo "Testing configurator for inclusion"

      - echo "Testing configurator for INCLUSION"

      - ../../src/scripts/model_maker.py ../../test_data/inclusion/config_dev.yml

      - echo "Running np_inclusion"

      - chmod +x np_inclusion
@@ -282,11 +290,14 @@ testing_stage: 
      - cd build/libnptm

      - export LD_LIBRARY_PATH="$LD_LIBRARY_PATH;$PWD"

      - cd ../sphere

      - export FFILE=../../test_data/sphere/OSPH

      - echo "Comparing output of SPHERE"

      - python3 ../../src/scripts/pycompare.py --no-progress --ffile $FFILE --cfile c_OSPH --html

      - export FFILE=../../test_data/sphere/OSPH

      - python3 ../../src/scripts/pycompare.py --no-progress --ffile $FFILE --cfile c_OSPH_dev --html

      - echo "Comparing output of parallel SPHERE"

      - export FFILE=../../test_data/sphere/OSPH_181_xi

      - python3 ../../src/scripts/pycompare.py --no-progress --ffile $FFILE --cfile c_OSPH_181_xi --html

      - echo "Checking consistency among legacy and HDF5 configuration files"

      - ../testing/test_TEDF ../../test_data/sphere/DEDFB c_TEDF c_TEDF.hd5

      - ../testing/test_TEDF ../../test_data/sphere/DEDFB c_TEDF_dev c_TEDF_dev.hd5

      - cd ../inclusion

      - echo "Comparing output of INCLUSION"

      - export FFILE=../../test_data/inclusion/OINCLU
@@ -316,7 +327,7 @@ testing_stage: 
      - export FFILE=../../test_data/cluster/OCLU_24

      - python3 ../../src/scripts/pycompare.py --no-progress --ffile $FFILE --cfile c_OCLU_24

      - rm -rf c_OCLU_24

      - echo "Checking consistency of HDF5 incluson output"

      - echo "Checking consistency of HDF5 inclusion output"

      - chmod u+x test_inclusion_outputs

      - ./test_inclusion_outputs

      - export FFILE=../../test_data/inclusion/OINCLU

src/include/Commons.h

+4 −4
Original line number Diff line number Diff line
@@ -80,12 +80,12 @@ protected: 
  int _nsph;

  //! \brief Maximum internal field expansion order.

  int _li;

  //! \brief Maximum number of layers in known sphere types.

  int _max_layers;

  //! \brief Number of different sphere types.

  int _num_configurations;

  //! \brief Total number of layers from all sphere types.

  int _num_layers;

  //! \brief Space for different sphere types.

  int _nl;

  //! \brief NHSPO = 2 * MAX(NPNT,NPNTTS) - 1

  int _nhspo;

  //! \brief Number of points for numerical integration in layered spheres.
@@ -184,12 +184,12 @@ public: 
  const int &nsph = _nsph;

  //! \brief Read-only view of maximum internal field expansion order.

  const int &li = _li;

  //! \brief Read-only view of the maximum number of layers in types.

  const int &max_layers = _max_layers;

  //! \brief Read-only view of number of different sphere types.

  const int &num_configurations = _num_configurations;

  //! \brief Read-only view of total number of layers from all sphere types.

  const int &num_layers = _num_layers;

  //! \brief Read-only view on the space for different sphere configurations.

  const int &nl = _nl;

  //! \brief Read-only view of NHSPO.

  const int &nhspo = _nhspo;

  //! \brief Read-only view on number of points for numerical integration in layered spheres.

src/include/IterationData.h

+16 −0
Original line number Diff line number Diff line
@@ -374,6 +374,22 @@ protected: 
  double *vec_zpv;

  

public:

  //! \brief Cosine of incident radiation azimuth.

  double cost;

  //! \brief Sine of incident radiation azimuth.

  double sint;

  //! \brief Cosine of incident radiation elevation.

  double cosp;

  //! \brief Sine of incident radiation elevation.

  double sinp;

  //! \brief Cosine of scattered radiation azimuth.

  double costs;

  //! \brief Sine of scattered radiation azimuth.

  double sints;

  //! \brief Cosine of scattered radiation elevation.

  double cosps;

  //! \brief Sine of scattered radiation elevation.

  double sinps;

  //! \brief Vacuum magnitude of wave vector.

  double vk;

  //! \brief Wave number.

src/libnptm/Commons.cpp

+30 −24
Original line number Diff line number Diff line
@@ -84,7 +84,12 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo 
  gcs = 0.0;

  _num_configurations = sconf->configurations;

  _num_layers = (sconf->use_external_sphere) ? 1 : 0;

  for (int nli = 0; nli < num_configurations; nli++) _num_layers += sconf->get_nshl(nli);

  _max_layers = 1;

  for (int nli = 0; nli < num_configurations; nli++) {

    int nl = sconf->get_nshl(nli);

    _num_layers += nl;

    if (nl >  _max_layers) _max_layers = nl;

  }



  vec_rmi = new dcomplex[_li * _nsph]();

  rmi = new dcomplex*[_li];
@@ -131,12 +136,12 @@ 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;

  // _nl = sconf->configurations;

  // if (_nsph == 1 && _nl == 1) _nl = 5;

  ris = new dcomplex[_nhspo]();

  dlri = new dcomplex[_nhspo]();

  vkt = new dcomplex[_nsph]();

  dc0 = new dcomplex[_nl]();

  dc0 = new dcomplex[_max_layers + 1]();

  vsz = new double[_nsph]();

  

  // >>> NEEDED BY SPHERE AND CLUSTER <<<
@@ -219,6 +224,7 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) { 
  _ndit = rhs._ndit;

  _ndm = rhs._ndm;

  gcs = rhs.gcs;

  _max_layers = rhs._max_layers;

  _num_configurations = rhs._num_configurations;

  _num_layers = rhs._num_layers;
@@ -271,7 +277,7 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) { 
  _npnt = rhs._npnt;

  _npntts = rhs._npntts;

  _nhspo = rhs._nhspo;

  _nl = rhs._nl;

  _max_layers = rhs._max_layers;

  ris = new dcomplex[_nhspo]();

  dlri = new dcomplex[_nhspo]();

  for (int ri = 0; ri < _nhspo; ri++) {
@@ -284,8 +290,8 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) { 
    vkt[vi] = rhs.vkt[vi];

    vsz[vi] = rhs.vsz[vi];

  }

  dc0 = new dcomplex[_nl]();

  for (int di = 0; di < _nl; di++) {

  dc0 = new dcomplex[_max_layers]();

  for (int di = 0; di < _max_layers; di++) {

    dc0[di] = rhs.dc0[di];

  }

  // >>> NEEDED BY SPHERE AND CLUSTER <<<
@@ -416,7 +422,7 @@ ParticleDescriptor::ParticleDescriptor(const mixMPI *mpidata) { 
  MPI_Bcast(&_npnt, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_npntts, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_nhspo, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_nl, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_max_layers, 1, MPI_INT, 0, MPI_COMM_WORLD);

  ris = new dcomplex[_nhspo];

  MPI_Bcast(ris, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  dlri = new dcomplex[_nhspo];
@@ -425,8 +431,8 @@ ParticleDescriptor::ParticleDescriptor(const mixMPI *mpidata) { 
  MPI_Bcast(vkt, _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  vsz = new double[_nsph];

  MPI_Bcast(vsz, _nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  dc0 = new dcomplex[_nl];

  MPI_Bcast(dc0, _nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  dc0 = new dcomplex[_max_layers];

  MPI_Bcast(dc0, _max_layers, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

}



void ParticleDescriptor::mpibcast(const mixMPI *mpidata) {
@@ -466,12 +472,12 @@ void ParticleDescriptor::mpibcast(const mixMPI *mpidata) { 
  MPI_Bcast(&_npnt, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_npntts, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_nhspo, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_nl, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&_max_layers, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(ris, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(dlri, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vkt, _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vsz, _nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(dc0, _nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(dc0, _max_layers, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  // >>> NEEDED BY SPHERE AND CLUSTER <<< //

  if (_class_type == SPHERE_TYPE || _class_type == CLUSTER_TYPE) {

    MPI_Bcast(vec_sas, 4 * _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
@@ -1216,8 +1222,8 @@ ParticleDescriptorInclusion::ParticleDescriptorInclusion(const mixMPI *mpidata) 
// >>> ParticleDescriptorSphere class implementation. <<< //

ParticleDescriptorSphere::ParticleDescriptorSphere(GeometryConfiguration *gconf, ScattererConfiguration *sconf) : ParticleDescriptor(gconf, sconf) {

  _class_type = SPHERE_TYPE;

  vec_sas = new dcomplex[4 * (_nsph + 1)]();

  vec_vints = new dcomplex[16 * (_nsph + 1)]();

  vec_sas = new dcomplex[4 * _nsph]();

  vec_vints = new dcomplex[16 * _nsph]();



  fsas = new dcomplex[_nsph];

  sscs = new double[_nsph]();
@@ -1230,7 +1236,7 @@ ParticleDescriptorSphere::ParticleDescriptorSphere(GeometryConfiguration *gconf, 
  sas = new dcomplex**[_nsph];

  vints = new dcomplex*[_nsph];

  for (int vi = 0; vi < _nsph; vi++) {

    vints[vi] = vec_vints + (16 * _nsph);

    vints[vi] = vec_vints + (16 * vi);

    sas[vi] = new dcomplex*[2];

    sas[vi][0] = vec_sas + (4 * vi);

    sas[vi][1] = vec_sas + (4 * vi) + 2;
@@ -1243,6 +1249,14 @@ ParticleDescriptorSphere::ParticleDescriptorSphere(const ParticleDescriptorSpher 
  for (int vsi = 0; vsi < 4 * _nsph; vsi++) vec_sas[vsi] = rhs.vec_sas[vsi];

  vec_vints = new dcomplex[16 * _nsph];

  for (int vvi = 0; vvi < 16 * _nsph; vvi++) vec_vints[vvi] = rhs.vec_vints[vvi];

  sas = new dcomplex**[_nsph];

  vints = new dcomplex*[_nsph];

  for (int vi = 0; vi < _nsph; vi++) {

    vints[vi] = vec_vints + (16 * vi);

    sas[vi] = new dcomplex*[2];

    sas[vi][0] = vec_sas + (4 * vi);

    sas[vi][1] = vec_sas + (4 * vi) + 2;

  }

  

  fsas = new dcomplex[_nsph];

  sscs = new double[_nsph];
@@ -1262,14 +1276,6 @@ ParticleDescriptorSphere::ParticleDescriptorSphere(const ParticleDescriptorSpher 
    sqabs[gi] = rhs.sqabs[gi];

    gcsv[gi] = rhs.gcsv[gi];

  }

  sas = new dcomplex**[_nsph];

  vints = new dcomplex*[_nsph];

  for (int vi = 0; vi < nsph; vi++) {

    vints[vi] = vec_vints + (16 * nsph);

    sas[vi] = new dcomplex*[2];

    sas[vi][0] = vec_sas + (4 * vi);

    sas[vi][1] = vec_sas + (4 * vi) + 2;

  }

}



#ifdef MPI_VERSION
@@ -1299,7 +1305,7 @@ ParticleDescriptorSphere::ParticleDescriptorSphere(const mixMPI *mpidata) : Part 
  sas = new dcomplex**[_nsph];

  vints = new dcomplex*[_nsph];

  for (int vi = 0; vi < nsph; vi++) {

    vints[vi] = vec_vints + (16 * nsph);

    vints[vi] = vec_vints + (16 * vi);

    sas[vi] = new dcomplex*[2];

    sas[vi][0] = vec_sas + (4 * vi);

    sas[vi][1] = vec_sas + (4 * vi) + 2;

src/sphere/sphere.cpp

+60 −27
Original line number Diff line number Diff line
@@ -548,6 +548,7 @@ int sphere_jxi488_cycle( 
  Logger *logger = new Logger(LOG_INFO);

  int oindex = 0;

  int jxi = jxi488 + 1;

  int jxindex = jxi - oi->first_xi;

  int nsph = gconf->number_of_spheres;

  int l_max = gconf->l_max;

  int in_pol = gconf->in_pol;
@@ -570,14 +571,14 @@ int sphere_jxi488_cycle( 
  if (idfc >= 0) {

    vk = xi * wn;

    vkarg = vk;

    oi->vec_vk[jxi488] = vk;

    oi->vec_xi[jxi488] = xi;

    oi->vec_vk[jxindex] = vk;

    oi->vec_xi[jxindex] = xi;

  } else { // IDFC < 0

    vk = sconf->xip * wn;

    vkarg = xi * vk;

    sqsfi = 1.0 / (xi * xi);

    oi->vec_vk[jxi488] = vk;

    oi->vec_xi[jxi488] = xi;

    oi->vec_vk[jxindex] = vk;

    oi->vec_xi[jxindex] = xi;

  }

  vtppoanp->append_line(VirtualBinaryLine(vk));

  double thsca = (gconf->isam > 1) ? sa->ths - sa->th : 0.0;
@@ -607,8 +608,9 @@ int sphere_jxi488_cycle( 
    if (nsh % 2 == 0) sid->c1->dc0[ici] = exdc;

    dme(l_max, i, npnt, npntts, vkarg, exdc, exri, sid->c1, jer, lcalc, sid->arg);

    if (jer != 0) {

      oi->vec_ier[jxi] = 1;

      oi->vec_ier[jxindex] = 1;

      oi->lcalc = lcalc;

      delete logger;

      return jer;

    }

  } // i132
@@ -635,7 +637,7 @@ int sphere_jxi488_cycle( 
    int i = i170 + 1;

    if (sid->c1->iog[i170] >= i) {

      last_configuration++;

      oindex = jxi488 * configurations + last_configuration - 1;

      oindex = jxindex * configurations + last_configuration - 1;

      double albeds = sid->c1->sscs[i170] / sid->c1->sexs[i170];

      sid->c1->sqscs[i170] *= sqsfi;

      sid->c1->sqabs[i170] *= sqsfi;
@@ -698,15 +700,15 @@ int sphere_jxi488_cycle( 
    } // End if iog[i170] >= i

  } // i170 loop

  if (nsph != 1) {

    oi->vec_fsat[jxi488] = sid->tfsas;

    oi->vec_fsat[jxindex] = sid->tfsas;

    double csch = 2.0 * vk * sqsfi / sid->c1->gcs;

    dcomplex s0 = sid->tfsas * exri;

    double qschu = csch * imag(s0);

    double pschu = csch * real(s0);

    double s0mag = cs0 * cabs(s0);

    oi->vec_qschut[jxi488] = qschu;

    oi->vec_pschut[jxi488] = pschu;

    oi->vec_s0magt[jxi488] = s0mag;

    oi->vec_qschut[jxindex] = qschu;

    oi->vec_pschut[jxindex] = pschu;

    oi->vec_s0magt[jxindex] = s0mag;

  }

  double th = sa->th;

  int done_dirs = 0;
@@ -721,12 +723,11 @@ int sphere_jxi488_cycle( 
    for (int jph484 = 0; jph484 < sa->nph; jph484++) {

      int jph = jph484 + 1;

      bool goto182 = (sa->nk == 1) && (jxi > 1);

      double cost, sint, cosp, sinp;

      if (!goto182) {

	upvmp(th, ph, 0, cost, sint, cosp, sinp, sid->u, sid->upmp, sid->unmp);

	upvmp(th, ph, 0, sid->cost, sid->sint, sid->cosp, sid->sinp, sid->u, sid->upmp, sid->unmp);

      }

      if (gconf->isam >= 0) {

	wmamp(0, cost, sint, cosp, sinp, in_pol, l_max, 0, nsph, sid->argi, sid->u, sid->upmp, sid->unmp, sid->c1);

	wmamp(0, sid->cost, sid->sint, sid->cosp, sid->sinp, in_pol, l_max, 0, nsph, sid->argi, sid->u, sid->upmp, sid->unmp, sid->c1);

	for (int i183 = 0; i183 < nsph; i183++) {

	  double rapr = sid->c1->sexs[i183] - sid->gaps[i183];

	  frx = rapr * sid->u[0];
@@ -752,16 +753,15 @@ int sphere_jxi488_cycle( 
	double phs = sa->phs;

	for (int jphs480 = 0; jphs480 < sa->nphs; jphs480++) {

	  int jphs = jphs480 + 1;

	  double costs, sints, cosps, sinps;

	  if (gconf->isam >= 1) {

	    phs = ph + phsph;

	    if (phs >= 360.0) phs -= 360.0;

	  }

	  bool goto190 = (nks == 1) && ((jxi > 1) || (jth > 1) || (jph > 1));

	  if (!goto190) {

	    upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, sid->us, sid->upsmp, sid->unsmp);

	    upvmp(ths, phs, icspnv, sid->costs, sid->sints, sid->cosps, sid->sinps, sid->us, sid->upsmp, sid->unsmp);

	    if (gconf->isam >= 0) {

	      wmamp(2, costs, sints, cosps, sinps, in_pol, l_max, 0, nsph, sid->args, sid->us, sid->upsmp, sid->unsmp, sid->c1);

	      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) {
@@ -772,10 +772,10 @@ int sphere_jxi488_cycle( 
	    );

	    if (gconf->isam < 0) {

	      wmasp(

		cost, sint, cosp, sinp, costs, sints, cosps, sinps,

		sid->u, sid->up, sid->un, sid->us, sid->ups, sid->uns,

		isq, ibf, in_pol, l_max, 0, nsph, sid->argi, sid->args,

		sid->c1

		sid->cost, sid->sint, sid->cosp, sid->sinp, sid->costs,

		sid->sints, sid->cosps, sid->sinps, sid->u, sid->up,

		sid->un, sid->us, sid->ups, sid->uns, isq, ibf, in_pol,

		l_max, 0, nsph, sid->argi, sid->args, sid->c1

	      );

	    }

	    for (int i193 = 0; i193 < 3; i193++) {
@@ -823,24 +823,24 @@ int sphere_jxi488_cycle( 
	  last_configuration = 0;

	  for (int ns226 = 0; ns226 < nsph; ns226++) {

	    int ns = ns226 + 1;

	    oindex = jxi488 * nsph * ndirs + nsph * done_dirs + ns226;

	    oindex = jxindex * nsph * ndirs + nsph * done_dirs + ns226;

	    oi->vec_dir_sas11[oindex] = sid->c1->sas[ns226][0][0];

	    oi->vec_dir_sas21[oindex] = sid->c1->sas[ns226][1][0];

	    oi->vec_dir_sas12[oindex] = sid->c1->sas[ns226][0][1];

	    oi->vec_dir_sas22[oindex] = sid->c1->sas[ns226][1][1];

	    oi->vec_dir_fx[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frx;

	    oi->vec_dir_fy[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = fry;

	    oi->vec_dir_fz[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frz;

	    oi->vec_dir_fx[jxindex * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frx;

	    oi->vec_dir_fy[jxindex * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = fry;

	    oi->vec_dir_fz[jxindex * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frz;

	    for (int i225 = 0; i225 < 16; i225++) sid->c1->vint[i225] = sid->c1->vints[ns226][i225];

	    mmulc(sid->c1->vint, sid->cmullr, sid->cmul);

	    for (int imul = 0; imul < 4; imul++) {

	      int muls_index = 16 * jxi488 * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;

	      int muls_index = 16 * jxindex * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;

	      for (int jmul = 0; jmul < 4; jmul++) {

		oi->vec_dir_muls[muls_index + jmul] = sid->cmul[imul][jmul];

	      }

	    }

	    for (int imul = 0; imul < 4; imul++) {

	      int muls_index = 16 * jxi488 * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;

	      int muls_index = 16 * jxindex * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;

	      for (int jmul = 0; jmul < 4; jmul++) {

		oi->vec_dir_mulslr[muls_index + jmul] = sid->cmullr[imul][jmul];

	      }
@@ -867,8 +867,9 @@ int sphere_jxi488_cycle( 
    } // jph484 loop on elevation

    th += sa->thstp;

  } // jth486 loop on azimuth

  oi->vec_jxi[jxi488] = jxi;

  oi->vec_jxi[jxindex] = jxi;

  logger->log("INFO: finished scale iteration " + to_string(jxi) + " of " + to_string(nxi) + ".\n");

  delete logger;

  return jer;

}
@@ -886,6 +887,14 @@ SphereIterationData::SphereIterationData( 
  c1 = new ParticleDescriptorSphere(gconf, sconf);

  argi = new double[1]();

  args = new double[1]();

  cost = 0.0;

  sint = 0.0;

  cosp = 0.0;

  sinp = 0.0;

  costs = 0.0;

  sints = 0.0;

  cosps = 0.0;

  sinps = 0.0;

  scan = 0.0;

  cfmp = 0.0;

  cfsp = 0.0;
@@ -959,6 +968,14 @@ SphereIterationData::SphereIterationData(const SphereIterationData &rhs) { 
  args = new double[1];

  argi[0] = rhs.argi[0];

  args[0] = rhs.args[0];

  cost = rhs.cost;

  sint = rhs.sint;

  cosp = rhs.cosp;

  sinp = rhs.sinp;

  costs = rhs.costs;

  sints = rhs.sints;

  cosps = rhs.cosps;

  sinps = rhs.sinps;

  scan = rhs.scan;

  cfmp = rhs.cfmp;

  cfsp = rhs.cfsp;
@@ -1095,6 +1112,14 @@ SphereIterationData::SphereIterationData(const mixMPI *mpidata, const int device 
  MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(&s0, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(&tfsas, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(&cost, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&sint, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&cosp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&sinp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&costs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&sints, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&cosps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&sinps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&scan, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&cfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
@@ -1173,6 +1198,14 @@ int SphereIterationData::mpibcast(const mixMPI *mpidata) { 
  MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(&s0, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(&tfsas, 1, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(&cost, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&sint, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&cosp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&sinp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&costs, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&sints, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&cosps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&sinps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&scan, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&cfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);