/* Copyright (C) 2024   INAF - Osservatorio Astronomico di Cagliari

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
   
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
   
   A copy of the GNU General Public License is distributed along with
   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.
 */

/*! \file inclusion.cpp
 *
 * \brief Implementation of the calculation for a sphere with inclusions.
 */
#include <chrono>
#include <cstdio>
#include <exception>
#include <fstream>
#include <hdf5.h>
#include <string>
#ifdef _OPENMP
#include <omp.h>
#endif
#ifdef USE_MPI
#ifndef MPI_VERSION
#include <mpi.h>
#endif
#endif
#ifdef USE_NVTX
#include <nvtx3/nvToolsExt.h>
#endif
#ifdef USE_MAGMA
#include <cuda_runtime.h>
#endif

#ifndef INCLUDE_TYPES_H_
#include "../include/types.h"
#endif

#ifndef INCLUDE_ERRORS_H_
#include "../include/errors.h"
#endif

#ifndef INCLUDE_LOGGING_H_
#include "../include/logging.h"
#endif

#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif

#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
#endif

#ifndef INCLUDE_SPH_SUBS_H_
#include "../include/sph_subs.h"
#endif

#ifndef INCLUDE_CLU_SUBS_H_
#include "../include/clu_subs.h"
#endif

#ifndef INCLUDE_INCLU_SUBS_H_
#include "../include/inclu_subs.h"
#endif

#ifndef INCLUDE_TRANSITIONMATRIX_H_
#include "../include/TransitionMatrix.h"
#endif

#ifndef INCLUDE_ALGEBRAIC_H_
#include "../include/algebraic.h"
#endif

#ifndef INCLUDE_LIST_H_
#include "../include/List.h"
#endif

#ifndef INCLUDE_FILE_IO_H_
#include "../include/file_io.h"
#endif

#ifndef INCLUDE_UTILS_H_
#include "../include/utils.h"
#endif

#ifndef INCLUDE_OUTPUTS_H_
#include "../include/outputs.h"
#endif

using namespace std;

// >>> InclusionIterationData header <<< //
/*! \brief A data structure representing the information used for a single scale
 * of the INCLUSION case.
 */
class InclusionIterationData {
protected:
  //! \brief Vectorized geometric asymmetry parameter components.
  double *vec_zpv;
  
public:
  //! \brief External layer index.
  int nimd;
  //! \brief External layer radius.
  double extr;
  
  //! \brief Pointer to a ParticleDescriptor structure.
  ParticleDescriptor *c1;
  //! \brief Vector of geometric asymmetry factors.
  double *gaps;
    //! \brief Components of extinction contribution to radiation torque on a single sphere along k.
  double **tqse;
  //! \brief Components of polarized extinction contribution to radiation torque on a single sphere along k.
  dcomplex **tqspe;
  //! \brief Components of scattering contribution to radiation torque on a single sphere along k.
  double **tqss;
  //! \brief Components of polarized scattering contribution to radiation torque on a single sphere along k.
  dcomplex **tqsps;
  //! \brief L-dependent coefficients of the geometric asymmetry parameter.
  double ****zpv;
  //! \brief Mean geometric asymmetry parameters.
  double **gapm;
  //! \brief Mean geometric asymmetry parameters referred to polarization plane.
  dcomplex **gappm;
  //! \brief Imaginary part of the harmonic functions argument.
  double *argi;
  //! \brief Argument of the harmonic functions referred to the scattering plane.
  double *args;
  //! \brief Geometric asymmetry parameters.
  double **gap;
  //! \brief Geometric asymmetry parameters referred to polarization plane.
  dcomplex **gapp;
  //! \brief Components of extinction contribution to radiation torque on the cluster along k.
  double **tqce;
  //! \brief Components of extinction contribution to radiation torque on the cluster along k referred to polarization plane.
  dcomplex **tqcpe;
  //! \brief Components of scattering contribution to radiation torque on the cluster along k.
  double **tqcs;
  //! \brief Components of scattering contribution to radiation torque on the cluster along k referred to polarization plane.
  dcomplex **tqcps;
  //! \brief Variation of unitary radiation vector. QUESTION: correct?
  double *duk;
  //! \brief Cluster extinction cross-section components referred to scattering plane.
  double **cextlr;
  //! \brief Cluster extinction cross-section components referred to meridional plane.
  double **cext;
  //! \brief Cluster Mueller Transformation Matrix components referred to scattering plane.
  double **cmullr;
  //! \brief Cluster Mueller Transformation Matrix components referred to meridional plane.
  double **cmul;
  //! \brief Geometric asymmetry parameter components.
  double *gapv;
  //! \brief Radiation extinction torque components.
  double *tqev;
  //! \brief Radiation scattering torque components.
  double *tqsv;
  //! \brief Incident unitary vector components.
  double *u;
  //! \brief Scattered unitary vector components.
  double *us;
  //! \brief Normal unitary vector components.
  double *un;
  //! \brief Normal scattered unitary vector components.
  double *uns;
  //! \brief Incident unitary vector components on polarization plane.
  double *up;
  //! \brief Scattered unitary vector components on polarization plane.
  double *ups;
  //! \brief Mean unitary vector components normal to polarization plane.
  double *unmp;
  //! \brief Mean scattered unitary vector components normal to polarization plane.
  double *unsmp;
  //! \brief Mean incident unitary vector components on polarization plane.
  double *upmp;
  //! \brief Mean scattered unitary vector components on polarization plane.
  double *upsmp;
  //! \brief Scattering angle.
  double scan;
  //! \brief Control parameter on incidence direction referred to meridional plane.
  double cfmp;
  //! \brief Control parameter on scattering direction referred to meridional plane.
  double sfmp;
  //! \brief Control parameter on incidence direction referred to scattering plane.
  double cfsp;
  //! \brief Control parameter on scattering direction referred to scattering plane.
  double sfsp;
  //! \brief SQSFI = XI^-2
  double sqsfi;
  //! \brief Vectorized scattering coefficient matrix.
  dcomplex *am_vector;
  //! \brief Scattering coefficient matrix.
  dcomplex **am;
  //! \brief Argument of harmonic functions. QUESTION: correct?
  dcomplex arg;
  //! \brief Vacuum magnitude of wave vector.
  double vk;
  //! \brief Wave number.
  double wn;
  //! \brief Normalization scale. QUESTION: correct?
  double xip;
  //! \brief Number of scales (wavelengths) to be computed.
  int number_of_scales;
  //! \brief Size of the block of scales handled by the current process.
  int xiblock;
  //! \brief Index of the first scale handled by the current process.
  int firstxi;
  //! \brief Index of the last scale handled by the current process.
  int lastxi;
  //! \brief ID of the GPU used by one MPI process.
  int proc_device;
  //! \brief Refinement mode selction flag.
  int refinemode;
  //! \brief Maximum number of refinement iterations.
  int maxrefiters;
  //! \brief Required accuracy level.
  double accuracygoal;

  /*! \brief `InclusionIterationData` default instance constructor.
   *
   * \param gconf: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` object.
   * \param sconf: `ScattererConfiguration *` Pointer to a `ScattererConfiguration` object.
   * \param mpidata: `mixMPI *` Pointer to a `mixMPI` object.
   * \param device_count: `const int` Number of offload devices available on the system.
   */
  InclusionIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf, const mixMPI *mpidata, const int device_count);
  
  /*! \brief `InclusionIterationData` copy constructor.
   *
   * \param rhs: `const InclusionIterationData &` Reference to the `InclusionIterationData` object to be copied.
   */
  InclusionIterationData(const InclusionIterationData& rhs);

#ifdef MPI_VERSION
  /*! \brief `InclusionIterationData` MPI constructor.
   *
   * \param mpidata: `const mixMPI *` Pointer to a `mixMPI` instance.
   * \param device_count: `const int` Number of offload devices available on the system.
   */
  InclusionIterationData(const mixMPI *mpidata, const int device_count);

  /*! \brief Broadcast over MPI the InclusionIterationData instance from MPI process 0 to all others.
   *
   * When using MPI, the initial InclusionIterationData instance created by MPI process 0
   * needs to be replicated on all other processes. This function sends it using
   * MPI broadcast calls. The MPI broadcast calls in this function must match those
   * in the constructor using the mixMPI pointer.
   *
   * \param mpidata: `mixMPI *` Pointer to the mpi structure used to do the MPI broadcast.
   */
  void mpibcast(const mixMPI *mpidata);
#endif

  /*! \brief `InclusionIterationData` instance destroyer.
   */
  ~InclusionIterationData();
};

// >>> End of InclusionIterationData header <<< //

// >>> InclusionIterationData implementation <<< //
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]();
  tqse = new double*[2];
  tqspe = new dcomplex*[2];
  tqss = new double*[2];
  tqsps = new dcomplex*[2];
  tqce = new double*[2];
  tqcpe = new dcomplex*[2];
  tqcs = new double*[2];
  tqcps = new dcomplex*[2];
  for (int ti = 0; ti < 2; ti++) {
    tqse[ti] = new double[c1->nsph]();
    tqspe[ti] = new dcomplex[c1->nsph]();
    tqss[ti] = new double[c1->nsph]();
    tqsps[ti] = new dcomplex[c1->nsph]();
    tqce[ti] = new double[3]();
    tqcpe[ti] = new dcomplex[3]();
    tqcs[ti] = new double[3]();
    tqcps[ti] = new dcomplex[3]();
  }
  gapv = new double[3]();
  gapp = new dcomplex*[3];
  gappm = new dcomplex*[3];
  gap = new double*[3];
  gapm = new double*[3];
  for (int gi = 0; gi < 3; gi++) {
    gapp[gi] = new dcomplex[2]();
    gappm[gi] = new dcomplex[2]();
    gap[gi] = new double[2]();
    gapm[gi] = new double[2]();
  }
  u = new double[3]();
  us = new double[3]();
  un = new double[3]();
  uns = new double[3]();
  up = new double[3]();
  ups = new double[3]();
  unmp = new double[3]();
  unsmp = new double[3]();
  upmp = new double[3]();
  upsmp = new double[3]();
  argi = new double[1]();
  args = new double[1]();
  duk = new double[3]();
  cextlr = new double*[4];
  cext = new double*[4];
  cmullr = new double*[4];;
  cmul = new double*[4];
  for (int ci = 0; ci < 4; ci++) {
    cextlr[ci] = new double[4]();
    cext[ci] = new double[4]();
    cmullr[ci] = new double[4]();
    cmul[ci] = new double[4]();
  }
  vec_zpv = new double[c1->lm * 12]();
  zpv = new double***[c1->lm];
  for (int zi = 0; zi < c1->lm; zi++) {
    zpv[zi] = new double**[12];
    for (int zj = 0; zj < 3; zj++) {
      zpv[zi][zj] = new double*[4];
      zpv[zi][zj][0] = vec_zpv + (zi * 12) + (zj * 4);
      zpv[zi][zj][1] = vec_zpv + (zi * 12) + (zj * 4) + 2;
    }
  }
  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);
  }
  
  arg = 0.0 + 0.0 * I;
  // These are suspect initializations
  scan = 0.0;
  cfmp = 0.0;
  sfmp = 0.0;
  cfsp = 0.0;
  sfsp = 0.0;
  // End of suspect initializations
  wn = sconf->wp / 3.0e8;
  xip = sconf->xip;
  sqsfi = 1.0;
  vk = 0.0;
  number_of_scales = sconf->number_of_scales;
  xiblock = (int) ceil(((double) (sconf->number_of_scales-1))/((double) mpidata->nprocs));
  lastxi = ((mpidata->rank+1) * xiblock)+1;
  firstxi = lastxi-xiblock+1;
  if (lastxi > sconf->number_of_scales) lastxi = sconf->number_of_scales;

  nimd = c1->nshl[0] + 1;
  c1->rc[0][nimd - 1] = c1->ros[0] * sconf->get_rcf(0, nimd - 1);
  extr = c1->rc[0][nimd - 1];
  const double pig = acos(0.0) * 2.0;
  c1->gcs = pig * extr * extr;
  
#ifdef USE_MAGMA
  proc_device = mpidata->rank % device_count;
#else
  proc_device = 0;
#endif

  // In the first iteration, if refinement is enabled, determine the number of refinement iterations required to arrive at the target accuracy (if achievable in a reasonable number of iterations)
  refinemode = 2;
  // maxrefiters and accuracygoal should be configurable and preferably set somewhere else
  maxrefiters = 20;
  accuracygoal = 1e-6;
}

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]();
  tqsv = new double[3]();
  for (int ti = 0; ti < 3; ti++) {
    tqev[ti] = rhs.tqev[ti];
    tqsv[ti] = rhs.tqsv[ti];
  }
  tqse = new double*[2];
  tqspe = new dcomplex*[2];
  tqss = new double*[2];
  tqsps = new dcomplex*[2];
  tqce = new double*[2];
  tqcpe = new dcomplex*[2];
  tqcs = new double*[2];
  tqcps = new dcomplex*[2];
  for (int ti = 0; ti < 2; ti++) {
    tqse[ti] = new double[c1->nsph]();
    tqspe[ti] = new dcomplex[c1->nsph]();
    tqss[ti] = new double[c1->nsph]();
    tqsps[ti] = new dcomplex[c1->nsph]();
    for (int tj = 0; tj < c1->nsph; tj++) {
      tqse[ti][tj] = rhs.tqse[ti][tj];
      tqspe[ti][tj] = rhs.tqspe[ti][tj];
      tqss[ti][tj] = rhs.tqss[ti][tj];
      tqsps[ti][tj] = rhs.tqsps[ti][tj];
    }
    tqce[ti] = new double[3]();
    tqcpe[ti] = new dcomplex[3]();
    tqcs[ti] = new double[3]();
    tqcps[ti] = new dcomplex[3]();
    for (int tj = 0; tj < 3; tj++) {
      tqce[ti][tj] = rhs.tqce[ti][tj];
      tqcpe[ti][tj] = rhs.tqcpe[ti][tj];
      tqcs[ti][tj] = rhs.tqcs[ti][tj];
      tqcps[ti][tj] = rhs.tqcps[ti][tj];
    }
  }
  gapv = new double[3]();
  gapp = new dcomplex*[3];
  gappm = new dcomplex*[3];
  gap = new double*[3];
  gapm = new double*[3];
  for (int gi = 0; gi < 3; gi++) {
    gapv[gi] = rhs.gapv[gi];
    gapp[gi] = new dcomplex[2]();
    gappm[gi] = new dcomplex[2]();
    gap[gi] = new double[2]();
    gapm[gi] = new double[2]();
    for (int gj = 0; gj < 2; gj++) {
      gapp[gi][gj] = rhs.gapp[gi][gj];
      gappm[gi][gj] = rhs.gappm[gi][gj];
      gap[gi][gj] = rhs.gap[gi][gj];
      gapm[gi][gj] = rhs.gapm[gi][gj];
    }
  }
  u = new double[3]();
  us = new double[3]();
  un = new double[3]();
  uns = new double[3]();
  up = new double[3]();
  ups = new double[3]();
  unmp = new double[3]();
  unsmp = new double[3]();
  upmp = new double[3]();
  upsmp = new double[3]();
  duk = new double[3]();
  for (int ui = 0; ui < 3; ui++) {
    u[ui] = rhs.u[ui];
    us[ui] = rhs.us[ui];
    un[ui] = rhs.un[ui];
    uns[ui] = rhs.uns[ui];
    up[ui] = rhs.up[ui];
    ups[ui] = rhs.ups[ui];
    unmp[ui] = rhs.unmp[ui];
    unsmp[ui] = rhs.unsmp[ui];
    upmp[ui] = rhs.upmp[ui];
    upsmp[ui] = rhs.upsmp[ui];
    duk[ui] = rhs.duk[ui];
  }
  argi = new double[1]();
  args = new double[1]();
  argi[0] = rhs.argi[0];
  args[0] = rhs.args[0];
  cextlr = new double*[4];
  cext = new double*[4];
  cmullr = new double*[4];;
  cmul = new double*[4];
  for (int ci = 0; ci < 4; ci++) {
    cextlr[ci] = new double[4]();
    cext[ci] = new double[4]();
    cmullr[ci] = new double[4]();
    cmul[ci] = new double[4]();
    for (int cj = 0; cj < 4; cj++) {
      cextlr[ci][cj] = rhs.cextlr[ci][cj];
      cext[ci][cj] = rhs.cext[ci][cj];
      cmullr[ci][cj] = rhs.cmullr[ci][cj];
      cmul[ci][cj] = rhs.cmul[ci][cj];
    }
  }
  vec_zpv = new double[c1->lm * 12];
  zpv = new double***[c1->lm];
  for (int zi = 0; zi < c1->lm; zi++) {
    zpv[zi] = new double **[12];
    for (int zj = 0; zj < 3; zj++) {
      zpv[zi][zj] = new double*[4];
      zpv[zi][zj][0] = vec_zpv + (zi * 12) + (zj * 4);
      zpv[zi][zj][1] = vec_zpv + (zi * 12) + (zj * 4) + 2;
      zpv[zi][zj][0][0] = rhs.zpv[zi][zj][0][0];
      zpv[zi][zj][0][1] = rhs.zpv[zi][zj][0][1];
      zpv[zi][zj][1][0] = rhs.zpv[zi][zj][1][0];
      zpv[zi][zj][1][1] = rhs.zpv[zi][zj][1][1];
    }
  }
  am_vector = new dcomplex[c1->ndm * c1->ndm];
  for (int ai = 0; ai < c1->ndm * c1->ndm; ai++) am_vector[ai] = rhs.am_vector[ai];
  am = new dcomplex*[c1->ndm];
  for (int ai = 0; ai < c1->ndm; ai++) {
    am[ai] = (am_vector + ai * c1->ndm);
  }
  
  arg = rhs.arg;
  // These are suspect initializations
  scan = rhs.scan;
  cfmp = rhs.cfmp;
  sfmp = rhs.sfmp;
  cfsp = rhs.cfsp;
  sfsp = rhs.sfsp;
  // End of suspect initializations
  wn = rhs.wn;
  xip = rhs.xip;
  sqsfi = rhs.sqsfi;
  vk = rhs.vk;
  firstxi = rhs.firstxi;
  lastxi = rhs.lastxi;
  xiblock = rhs.xiblock;
  number_of_scales = rhs.number_of_scales;

  nimd = rhs.nimd;
  extr = rhs.extr;
  
  proc_device = rhs.proc_device;
  refinemode = rhs.refinemode;
  maxrefiters = rhs.maxrefiters;
  accuracygoal = rhs.accuracygoal;
}

#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]();
  tqsv = new double[3]();
  MPI_Bcast(tqev, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(tqsv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  tqse = new double*[2];
  tqspe = new dcomplex*[2];
  tqss = new double*[2];
  tqsps = new dcomplex*[2];
  tqce = new double*[2];
  tqcpe = new dcomplex*[2];
  tqcs = new double*[2];
  tqcps = new dcomplex*[2];
  for (int ti = 0; ti < 2; ti++) {
    tqse[ti] = new double[c1->nsph]();
    tqspe[ti] = new dcomplex[c1->nsph]();
    tqss[ti] = new double[c1->nsph]();
    tqsps[ti] = new dcomplex[c1->nsph]();
    MPI_Bcast(tqse[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqspe[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqss[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqsps[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    tqce[ti] = new double[3]();
    tqcpe[ti] = new dcomplex[3]();
    tqcs[ti] = new double[3]();
    tqcps[ti] = new dcomplex[3]();
    MPI_Bcast(tqce[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcpe[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcs[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcps[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  gapv = new double[3]();
  gapp = new dcomplex*[3];
  gappm = new dcomplex*[3];
  gap = new double*[3];
  gapm = new double*[3];
  MPI_Bcast(gapv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  for (int gi = 0; gi < 3; gi++) {
    gapp[gi] = new dcomplex[2]();
    gappm[gi] = new dcomplex[2]();
    gap[gi] = new double[2]();
    gapm[gi] = new double[2]();
    MPI_Bcast(gapp[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(gappm[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(gap[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(gapm[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  }
  u = new double[3]();
  us = new double[3]();
  un = new double[3]();
  uns = new double[3]();
  up = new double[3]();
  ups = new double[3]();
  unmp = new double[3]();
  unsmp = new double[3]();
  upmp = new double[3]();
  upsmp = new double[3]();
  duk = new double[3]();
  MPI_Bcast(u, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(us, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(un, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(uns, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(up, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ups, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(unmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(unsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(upmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(upsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(duk, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  argi = new double[1]();
  args = new double[1]();
  MPI_Bcast(argi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(args, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  cextlr = new double*[4];
  cext = new double*[4];
  cmullr = new double*[4];;
  cmul = new double*[4];
  for (int ci = 0; ci < 4; ci++) {
    cextlr[ci] = new double[4]();
    cext[ci] = new double[4]();
    cmullr[ci] = new double[4]();
    cmul[ci] = new double[4]();
    MPI_Bcast(cextlr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cext[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cmullr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cmul[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  }
  vec_zpv = new double[c1->lm * 12];
  MPI_Bcast(vec_zpv, c1->lm * 12, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  zpv = new double***[c1->lm];
  for (int zi = 0; zi < c1->lm; zi++) {
    zpv[zi] = new double **[12];
    for (int zj = 0; zj < 3; zj++) {
      zpv[zi][zj] = new double*[4];
      zpv[zi][zj][0] = vec_zpv + (zi * 12) + (zj * 4);
      zpv[zi][zj][1] = vec_zpv + (zi * 12) + (zj * 4) + 2;
    }
  }
  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);
    MPI_Bcast(am[ai], c1->ndm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 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(&sfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&wn, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sqsfi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&vk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);
  lastxi = ((mpidata->rank+1) * xiblock)+1;
  firstxi = lastxi-xiblock+1;
  if (lastxi > number_of_scales) lastxi = number_of_scales;

  MPI_Bcast(&nimd, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&extr, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  
#ifdef USE_MAGMA
  proc_device = mpidata->rank % device_count;
#else
  proc_device = 0;
#endif
  MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}

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);
  for (int ti = 0; ti < 2; ti++) {
    MPI_Bcast(tqse[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqspe[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqss[ti], c1->nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqsps[ti], c1->nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqce[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcpe[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcs[ti], 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(tqcps[ti], 3, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(gapv, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  for (int gi = 0; gi < 3; gi++) {
    MPI_Bcast(gapp[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(gappm[gi], 2, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
    MPI_Bcast(gap[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(gapm[gi], 2, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(u, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(us, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(un, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(uns, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(up, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(ups, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(unmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(unsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(upmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(upsmp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(duk, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(argi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(args, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  for (int ci = 0; ci < 4; ci++) {
    MPI_Bcast(cextlr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cext[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cmullr[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(cmul[ci], 4, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(vec_zpv, c1->lm * 12, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  // since MPI expects an int argument for the number of elements to transfer in one go, transfer am one row at a time
  for (int ai = 0; ai < c1->ndm; ai++) {
    MPI_Bcast(am[ai], c1->ndm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
  }
  MPI_Bcast(&arg, 1, MPI_C_DOUBLE_COMPLEX, 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(&sfmp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&cfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sfsp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&wn, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&xip, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sqsfi, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&vk, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&xiblock, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&number_of_scales, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&nimd, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&extr, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&refinemode, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&maxrefiters, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&accuracygoal, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
}
#endif

InclusionIterationData::~InclusionIterationData() {
  const int nsph = c1->nsph;
  delete[] am_vector;
  delete[] am;
  for (int zi = 0; zi < c1->lm; zi++) {
    for (int zj = 0; zj < 3; zj++) {
      delete[] zpv[zi][zj];
    }
    delete[] zpv[zi];
  }
  delete[] zpv;
  delete[] vec_zpv;
  delete c1;
  delete[] gaps;
  for (int ti = 1; ti > -1; ti--) {
    delete[] tqse[ti];
    delete[] tqss[ti];
    delete[] tqspe[ti];
    delete[] tqsps[ti];
    delete[] tqce[ti];
    delete[] tqcpe[ti];
    delete[] tqcs[ti];
    delete[] tqcps[ti];
  }
  delete[] tqse;
  delete[] tqss;
  delete[] tqspe;
  delete[] tqsps;
  delete[] tqce;
  delete[] tqcpe;
  delete[] tqcs;
  delete[] tqcps;
  delete[] tqev;
  delete[] tqsv;
  for (int gi = 2; gi > -1; gi--) {
    delete[] gapp[gi];
    delete[] gappm[gi];
    delete[] gap[gi];
    delete[] gapm[gi];
  }
  delete[] gapp;
  delete[] gappm;
  delete[] gap;
  delete[] gapm;
  delete[] gapv;
  delete[] u;
  delete[] us;
  delete[] un;
  delete[] uns;
  delete[] up;
  delete[] ups;
  delete[] unmp;
  delete[] unsmp;
  delete[] upmp;
  delete[] upsmp;
  delete[] argi;
  delete[] args;
  delete[] duk;
  for (int ci = 3; ci > -1; ci--) {
    delete[] cextlr[ci];
    delete[] cext[ci];
    delete[] cmullr[ci];
    delete[] cmul[ci];
  }
  delete[] cextlr;
  delete[] cext;
  delete[] cmullr;
  delete[] cmul;
}
// >>> End of InclusionIterationData implementation <<< //

/*! \brief Main calculation loop.
 *
 *  The solution of the scattering problem for different wavelengths is an
 *  embarrasingly parallel task. This function, therefore, collects all the
 *  operations that can be independently executed by different processes,
 *  after the configuration stage and the first calculation loop have been
 *  executed.
 *
 *  \param jxi488: `int` Wavelength loop index.
 *  \param sconf: `ScattererConfiguration *` Pointer to a `ScattererConfiguration` object.
 *  \param gconf: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` object.
 *  \param sa: `ScatteringAngles *` Pointer to a `ScatteringAngles` object.
 *  \param cid: `InclusionIterationData *` Pointer to an `InclusionIterationData` object.
 *  \param oi: `InclusionOutputInfo *` Pointer to an `InclusionOutputInfo` object.
 *  \param output_path: `const string &` Path to the output directory.
 *  \param vtppoanp: `VirtualBinaryFile *` Pointer to a `VirtualBinaryFile` object.
 */
int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, InclusionOutputInfo *output, const string& output_path, VirtualBinaryFile *vtppoanp);

/*! \brief C++ implementation of INCLU
 *
 * \param config_file: `string` Name of the configuration file.
 * \param data_file: `string` Name of the input data file.
 * \param output_path: `string` Directory to write the output files in.
 * \param mpidata: `mixMPI *` Pointer to an instance of MPI data settings.
 */
void inclusion(const string& config_file, const string& data_file, const string& output_path, const mixMPI *mpidata) {
  chrono::time_point<chrono::high_resolution_clock> t_start = chrono::high_resolution_clock::now();
  chrono::duration<double> elapsed;
  string message;
  string timing_name;
  FILE *timing_file;
  Logger *time_logger;
  if (mpidata->rank == 0) {
    timing_name = output_path + "/c_timing_mpi"+ to_string(mpidata->rank) +".log";
    timing_file = fopen(timing_name.c_str(), "w");
    time_logger = new Logger(LOG_DEBG, timing_file);
  }
  Logger *logger = new Logger(LOG_DEBG);
  int device_count = 0;
  //===========
  // Initialise MAGMA
  //===========
#ifdef USE_MAGMA
  const magma_int_t d_array_max_size = 32; // TEMPORARY: can become configurable parameter
  magma_device_t *device_array = new magma_device_t[d_array_max_size];
  magma_int_t num_devices;
  magma_getdevices(device_array, d_array_max_size, &num_devices);
  device_count = (int)num_devices;
  delete[] device_array;
  message = "DEBUG: Proc-" + to_string(mpidata->rank) + " found " + to_string(device_count) + " GPU ";
  if (device_count > 1) message += "devices.\n";
  else message += "device.\n";
  logger->log(message, LOG_DEBG);
  logger->log("INFO: Process " + to_string(mpidata->rank) + " initializes MAGMA.\n");
  magma_int_t magma_result = magma_init();
  if (magma_result != MAGMA_SUCCESS) {
    logger->err("ERROR: Process " + to_string(mpidata->rank) + " failed to initilize MAGMA.\n");
    logger->err("PROC-" + to_string(mpidata->rank) + ": MAGMA error code " + to_string(magma_result) + "\n");
    if (mpidata->rank == 0) {
      fclose(timing_file);
      delete time_logger;
    }
    delete logger;
    return;
  }
#endif // end MAGMA initialisation
  
  //===========================
  // the following only happens on MPI process 0
  //===========================
  if (mpidata->rank == 0) {
#ifdef USE_NVTX
    nvtxRangePush("Set up");
#endif
    //=======================
    // Initialise sconf from configuration file
    //=======================
    logger->log("INFO: making legacy configuration...", LOG_INFO);
    ScattererConfiguration *sconf = NULL;
    try {
      sconf = ScattererConfiguration::from_dedfb(config_file);
    } catch(const OpenConfigurationFileException &ex) {
      logger->err("\nERROR: failed to open scatterer configuration file.\n");
      string message = "FILE: " + string(ex.what()) + "\n";
      logger->err(message);
      fclose(timing_file);
      delete time_logger;
      delete logger;
      return;
    }
    sconf->write_formatted(output_path + "/c_OEDFB");
    sconf->write_binary(output_path + "/c_TEDF");
    sconf->write_binary(output_path + "/c_TEDF.hd5", "HDF5");
    // end scatterer initialisation

    //========================
    // Initialise gconf from configuration files
    //========================
    GeometryConfiguration *gconf = NULL;
    try {
      gconf = GeometryConfiguration::from_legacy(data_file);
    } catch (const OpenConfigurationFileException &ex) {
      logger->err("\nERROR: failed to open geometry configuration file.\n");
      string message = "FILE: " + string(ex.what()) + "\n";
      logger->err(message);
      if (sconf) delete sconf;
      fclose(timing_file);
      delete time_logger;
      delete logger;
      return;
    }
    logger->log(" done.\n", LOG_INFO);
    //end gconf initialisation

#ifdef USE_NVTX
    nvtxRangePop();
#endif
    int s_nsph = sconf->number_of_spheres;
    int nsph = gconf->number_of_spheres;
    // Sanity check on number of sphere consistency, should always be verified
    if (s_nsph == nsph) {
      // Shortcuts to variables stored in configuration objects
      ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);
      double wp = sconf->wp;
      // 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);
      thdps(cid->c1->lm, cid->zpv);
      double exdc = sconf->exdc;
      double exri = sqrt(exdc);

      // Create an empty bynary file
      VirtualBinaryFile *vtppoanp = new VirtualBinaryFile();
      string tppoan_name = output_path + "/c_TPPOAN";
#ifdef USE_MAGMA
      logger->log("INFO: using MAGMA calls.\n", LOG_INFO);
#elif defined USE_LAPACK
      logger->log("INFO: using LAPACK calls.\n", LOG_INFO);
#else
      logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);
#endif
      int iavm = gconf->iavm;
      int isam = gconf->isam;
      int inpol = gconf->in_pol;
      int nxi = sconf->number_of_scales;
      int nth = p_scattering_angles->nth;
      int nths = p_scattering_angles->nths;
      int nph = p_scattering_angles->nph;
      int nphs = p_scattering_angles->nphs;
      
      //========================
      // write a block of info to virtual binary file
      //========================
      vtppoanp->append_line(VirtualBinaryLine(iavm));
      vtppoanp->append_line(VirtualBinaryLine(isam));
      vtppoanp->append_line(VirtualBinaryLine(inpol));
      vtppoanp->append_line(VirtualBinaryLine(nxi));
      vtppoanp->append_line(VirtualBinaryLine(nth));
      vtppoanp->append_line(VirtualBinaryLine(nph));
      vtppoanp->append_line(VirtualBinaryLine(nths));
      vtppoanp->append_line(VirtualBinaryLine(nphs));
      if (sconf->idfc < 0) {
	cid->vk = cid->xip * cid->wn;
	p_output->vec_vk[0] = cid->vk;
      }
      
      // do the first iteration on jxi488 separately, since it seems to be different from the others
      int jxi488 = 1;
      int initialmaxrefiters = cid->maxrefiters;
      
      chrono::time_point<chrono::high_resolution_clock> start_iter_1 = chrono::high_resolution_clock::now();
#ifdef USE_NVTX
      nvtxRangePush("First iteration");
#endif
      // use these pragmas, which should have no effect on parallelism, just to push OMP nested levels at the same level also in the first wavelength iteration
      int jer = 0;
#pragma omp parallel
      {
#pragma omp single
	{
	  jer = inclusion_jxi488_cycle(jxi488, sconf, gconf, p_scattering_angles, cid, p_output, output_path, vtppoanp);
	}
      }
#ifdef USE_NVTX
      nvtxRangePop();
#endif
      chrono::time_point<chrono::high_resolution_clock> end_iter_1 = chrono::high_resolution_clock::now();
      elapsed = start_iter_1 - t_start;
      string message = "INFO: Calculation setup took " + to_string(elapsed.count()) + "s.\n";
      logger->log(message);
      time_logger->log(message);
      elapsed = end_iter_1 - start_iter_1;
      message = "INFO: First iteration took " + to_string(elapsed.count()) + "s.\n";
      logger->log(message);
      time_logger->log(message);
      /* for the next iterations, just always do maxiter iterations, assuming the accuracy is good enough */
      cid->refinemode = 0;
      /* add an extra iteration for margin, if this does not exceed initialmaxrefiters */
      // if (cid->maxrefiters < initialmaxrefiters) cid->maxrefiters++;
      if (jer != 0) {
	// First loop failed. Halt the calculation.
	fclose(timing_file);
	delete time_logger;
	delete p_output;
	delete p_scattering_angles;
	delete cid;
	delete logger;
	delete sconf;
	delete gconf;
	return;
      }

      //==================================================
      // do the first outputs here, so that I open here the new files, afterwards I only append
      //==================================================
      vtppoanp->write_to_disk(output_path + "/c_TPPOAN");
      delete vtppoanp;
      
      // here go the calls that send data to be duplicated on other MPI processes from process 0 to others, using MPI broadcasts, but only if MPI is actually used
#ifdef MPI_VERSION
      if (mpidata->mpirunning) {
	gconf->mpibcast(mpidata);
	sconf->mpibcast(mpidata);	    
	cid->mpibcast(mpidata);
	p_scattering_angles->mpibcast(mpidata);
      }	
#endif
      // Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled
      int ompnumthreads = 1;
      // this is for MPI process 0 (or even if we are not using MPI at all)
      int myjxi488startoffset = 0;
      int myMPIstride = ompnumthreads;
      int myMPIblock = ompnumthreads;
      // Define here shared arrays of virtual ascii and binary files, so that thread 0 will be able to access them all later
      InclusionOutputInfo **p_outarray = NULL;
      VirtualBinaryFile **vtppoanarray = NULL;

#ifdef USE_NVTX
      nvtxRangePush("Parallel loop");
#endif

      //===========================================
      // open the OpenMP parallel context, so each thread can initialise its stuff
      //===========================================
#pragma omp parallel
      {
	// Create and initialise this variable here, so that if OpenMP is enabled it is local to the thread, and if OpenMP is not enabled it has a well-defiled value anyway
	int myompthread = 0;

#ifdef _OPENMP
	// If OpenMP is enabled, give actual values to myompthread and ompnumthreads, and open thread-local output files
	myompthread = omp_get_thread_num();
	if (myompthread == 0) ompnumthreads = omp_get_num_threads();
#endif

	if (myompthread == 0) {
	  // Initialise some shared variables only on thread 0
	  p_outarray = new InclusionOutputInfo*[ompnumthreads];
	  vtppoanarray = new VirtualBinaryFile*[ompnumthreads];
	  myMPIblock = ompnumthreads;
	  myMPIstride = myMPIblock;
	}

#ifdef MPI_VERSION
	if (myompthread == 0) {
	  if (mpidata->mpirunning) {
	    // only go through this if MPI has been actually used
	    for (int rr=1; rr<mpidata->nprocs; rr++) {
	      // individually send their respective starting points to other MPI processes: they start immediately after the frequencies computed by previous processes so far
	      int remotejxi488startoffset = myMPIstride;
	      MPI_Send(&remotejxi488startoffset, 1, MPI_INT, rr, 3, MPI_COMM_WORLD);
	      int remoteMPIblock;
	      MPI_Recv(&remoteMPIblock, 1, MPI_INT, rr, 3, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
	      // update myMPIstride to include the ones due to MPI process rr
	      myMPIstride += remoteMPIblock;
	    }
	    // now I know the total myMPIstride, I can send it to all processes
	    MPI_Bcast(&myMPIstride, 1, MPI_INT, 0, MPI_COMM_WORLD);
	  }
	}
#endif
	// add an omp barrier to make sure that the global variables defined by thread 0 are known to all threads below this
#pragma omp barrier

	// To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
	InclusionIterationData *cid_2 = NULL;
	InclusionOutputInfo *p_output_2 = NULL;
	VirtualBinaryFile *vtppoanp_2 = NULL;
	// for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
	if (myompthread == 0) {
	  cid_2 = cid;
	  // OMP thread 0 of MPI process 0 holds the pointer to the full output structure
	  p_output_2 = p_output;
	  p_outarray[0] = p_output_2;
	} else {
	  // this is not thread 0, so do create fresh copies of all local variables
	  cid_2 = new InclusionIterationData(*cid);
	}
	// make sure all threads align here: I don't want the following loop to accidentally start for thread 0, possibly modifying some variables before they are copied by all other threads
	if (myompthread==0) {
	  logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
	}
#pragma omp barrier
	// ok, now I can actually start the parallel calculations
	for (int ixi488=2; ixi488<=cid_2->number_of_scales; ixi488 +=myMPIstride) {
	  // the parallel loop over MPI processes covers a different set of indices for each thread
#pragma omp barrier
	  int myjxi488 = ixi488+myompthread;
	  // each thread opens new virtual files and stores their pointers in the shared array
	  vtppoanp_2 = new VirtualBinaryFile();
	  // each thread puts a copy of the pointers to its virtual files in the shared arrays
	  vtppoanarray[myompthread] = vtppoanp_2;
#pragma omp barrier

	  // each MPI process handles a number of contiguous scales corresponding to its number of OMP threads at this omp level of parallelism
	  if (myjxi488 <= cid_2->number_of_scales) {
	    if (myompthread > 0) {
	      // UPDATE: non-0 threads need to allocate memory for one scale at a time.
	      p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, 1);
	      p_outarray[myompthread] = p_output_2;
	    }
	    int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
	  } else {
	    if (myompthread > 0) {
	      // If there is no input for this thread, set output pointer to NULL.
	      p_outarray[myompthread] = NULL;
	    }
	  }
#pragma omp barrier

#ifdef USE_NVTX
	  nvtxRangePush("Output concatenation");
#endif
#pragma omp barrier
	  // threads different from 0 append their virtual files to the one of thread 0, and delete them
	  if (myompthread == 0) {
	    for (int ti=1; ti<ompnumthreads; ti++) {
	      if (p_outarray[ti] != NULL) {
		p_outarray[0]->insert(*(p_outarray[ti]));
		delete p_outarray[ti];
		p_outarray[ti] = NULL;
	      }
	      vtppoanarray[0]->append(*(vtppoanarray[ti]));
	      delete vtppoanarray[ti];
	    }
	  }
#pragma omp barrier
	  //==============================================
	  // Collect all virtual files on thread 0 of MPI process 0, and append them to disk
	  //==============================================
	  if (myompthread == 0) {
	    // thread 0 writes its virtual files, now including contributions from all threads, to disk, and deletes them
	    // p_outarray[0]->append_to_disk(output_path + "/c_OINCLU");
	    // delete p_outarray[0];
	    vtppoanarray[0]->append_to_disk(output_path + "/c_TPPOAN");
	    delete vtppoanarray[0];

#ifdef MPI_VERSION
	    if (mpidata->mpirunning) {
	      // only go through this if MPI has been actually used
	      for (int rr=1; rr<mpidata->nprocs; rr++) {
		// get the data from process rr by receiving it in total memory structure
		p_outarray[0]->mpireceive(mpidata, rr);
		// get the data from process rr, creating a new virtual ascii file
		// VirtualAsciiFile *p_output = new VirtualAsciiFile(mpidata, rr);
		// append to disk and delete virtual ascii file
		// p_output->append_to_disk(output_path + "/c_OINCLU");
		// delete p_output;
		
		// get the data from process rr, creating a new virtual binary file
		VirtualBinaryFile *vtppoanp = new VirtualBinaryFile(mpidata, rr);
		// append to disk and delete virtual binary file
		vtppoanp->append_to_disk(output_path + "/c_TPPOAN");
		delete vtppoanp;
		int test = MPI_Barrier(MPI_COMM_WORLD);
	      }
	    }
#endif
	  }
	  // end block writing to disk
#ifdef USE_NVTX
	  nvtxRangePop();
#endif
#pragma omp barrier

	} // close strided loop running on MPI processes, ixi488 loop
	// delete the shared arrays I used to make available to thread 0 the virtual files of other threads
#pragma omp barrier
	if (myompthread == 0) {
	  delete[] p_outarray;
	  delete[] vtppoanarray;
	}
	{
	  string message = "INFO: Closing thread-local output files of thread " + to_string(myompthread) + " and syncing threads.\n";
	  logger->log(message);
	}
#ifdef USE_NVTX
	nvtxRangePop();
#endif
	delete cid_2;
      }
      delete p_scattering_angles;
      p_output->write(output_path + "/c_OINCLU.hd5", "HDF5");
      p_output->write(output_path + "/c_OINCLU", "LEGACY");
      delete p_output;
    } // closes s_nsph == nsph check
	  
    else { // Sphere number inconsistency error.
      throw UnrecognizedConfigurationException(
	"Inconsistent geometry and scatterer configurations."
      );
    }

    delete sconf;
    delete gconf;
#ifdef USE_MAGMA
    logger->log("INFO: Process " + to_string(mpidata->rank) + " finalizes MAGMA.\n");
    magma_finalize();
#endif
    chrono::time_point<chrono::high_resolution_clock> t_end = chrono::high_resolution_clock::now();
    elapsed = t_end - t_start;
    string message = "INFO: Calculation lasted " + to_string(elapsed.count()) + "s.\n";
    logger->log(message);
    logger->log("Finished: output written to " + output_path + "/c_OINCLU\n");
    time_logger->log(message);
    fclose(timing_file);
    delete time_logger;
  } // end instructions block of MPI process 0
  
    //===============================
    // instruction block for MPI processes different from 0
    //===============================
#ifdef MPI_VERSION
  else {
    // here go the code for MPI processes other than 0
    // copy gconf, sconf, cid and p_scattering_angles from MPI process 0
    GeometryConfiguration *gconf = new GeometryConfiguration(mpidata);
    ScattererConfiguration *sconf = new ScattererConfiguration(mpidata);
    InclusionIterationData *cid = new InclusionIterationData(mpidata, device_count);
    ScatteringAngles *p_scattering_angles = new ScatteringAngles(mpidata);

    // Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled
    int ompnumthreads = 1;
    InclusionOutputInfo **p_outarray = NULL;
    VirtualBinaryFile **vtppoanarray = NULL;
    int myjxi488startoffset;
    int myMPIstride = ompnumthreads;
    int myMPIblock = ompnumthreads;
      
#pragma omp parallel
    {
      // Create and initialise this variable here, so that if OpenMP is enabled it is local to the thread, and if OpenMP is not enabled it has a well-defiled value anyway
      int myompthread = 0;
#ifdef _OPENMP
      // If OpenMP is enabled, give actual values to myompthread and ompnumthreads, and open thread-local output files
      myompthread = omp_get_thread_num();
      if (myompthread == 0) ompnumthreads = omp_get_num_threads();
#endif
      if (myompthread == 0) {
	// receive the start parameter from MPI process 0
	MPI_Recv(&myjxi488startoffset, 1, MPI_INT, 0, 3, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
	// send my number of omp threads to process 0
	MPI_Send(&ompnumthreads, 1, MPI_INT, 0, 3, MPI_COMM_WORLD);
	// receive myMPIstride sent by MPI process 0 to all processes
	MPI_Bcast(&myMPIstride, 1, MPI_INT, 0, MPI_COMM_WORLD);
	// allocate virtual files for each thread
	p_outarray = new InclusionOutputInfo*[ompnumthreads];
	vtppoanarray = new VirtualBinaryFile*[ompnumthreads];
      }
#pragma omp barrier
      // To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
      InclusionIterationData *cid_2 = NULL;
      InclusionOutputInfo *p_output_2 = NULL;
      VirtualBinaryFile *vtppoanp_2 = NULL;
      // PLACEHOLDER
      // for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
      if (myompthread == 0) {
	cid_2 = cid;
      } else {
	// this is not thread 0, so do create fresh copies of all local variables
	cid_2 = new InclusionIterationData(*cid);
      }
      // make sure all threads align here: I don't want the following loop to accidentally start for thread 0, possibly modifying some variables before they are copied by all other threads
#pragma omp barrier
      // ok, now I can actually start the parallel calculations
      for (int ixi488=2; ixi488<=cid_2->number_of_scales; ixi488 +=myMPIstride) {
	// the parallel loop over MPI processes covers a different set of indices for each thread
#pragma omp barrier
	int myjxi488 = ixi488 + myjxi488startoffset + myompthread;
	// each thread opens new virtual files and stores their pointers in the shared array
	vtppoanp_2 = new VirtualBinaryFile();
	// each thread puts a copy of the pointers to its virtual files in the shared arrays
	vtppoanarray[myompthread] = vtppoanp_2;
#pragma omp barrier
	if (myompthread==0) logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
	// ok, now I can actually start the parallel calculations
	// each MPI process handles a number of contiguous scales corresponding to its number of OMP threads at this omp level of parallelism
	if (myjxi488 <= cid_2->number_of_scales) {
	  if (myompthread > 0) {
	    // UPDATE: non-0 threads need to allocate memory for one scale at a time.
	    p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, 1);
	    p_outarray[myompthread] = p_output_2;
	  } else {
	    // Thread 0 of non-zero MPI processes needs to allocate memory for the
	    // output of all threads.
	    p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, ompnumthreads);
	    p_outarray[0] = p_output_2;
	  }
	  int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
	} else {
	  if (myompthread > 0) {
	    // If there is no input for this thread, set the output pointer to NULL.
	    p_outarray[myompthread] = NULL;
	  }	  
	}

#pragma omp barrier
	// threads different from 0 append their virtual files to the one of thread 0, and delete them
	if (myompthread == 0) {
	  for (int ti=1; ti<ompnumthreads; ti++) {
	    if (p_outarray[ti] != NULL) {
	      p_outarray[0]->insert(*(p_outarray[ti]));
	      delete p_outarray[ti];
	      p_outarray[ti] = NULL;
	    }
	    vtppoanarray[0]->append(*(vtppoanarray[ti]));
	    delete vtppoanarray[ti];
	  }
	  // thread 0 sends the collected virtualfiles to thread 0 of MPI process 0, then deletes them
	  for (int rr=1; rr<mpidata->nprocs; rr++) {
	    if (rr == mpidata->rank) {
	      p_outarray[0]->mpisend(mpidata);
	      delete p_outarray[0];
	      vtppoanarray[0]->mpisend(mpidata);
	      delete vtppoanarray[0];
	    }
	    int test = MPI_Barrier(MPI_COMM_WORLD);
	  }
	}
      } // close strided loop running on MPI processes
      
	// Clean memory
#pragma omp barrier
      if (myompthread == 0) {
	delete[] p_outarray;
	delete[] vtppoanarray;
      }
      delete cid_2;

    } // close pragma omp parallel
    delete p_scattering_angles;
    delete sconf;
    delete gconf;
#endif
#ifdef USE_MAGMA
    logger->log("INFO: Process " + to_string(mpidata->rank) + " finalizes MAGMA.\n");
    magma_finalize();
#endif
    delete logger;
#ifdef MPI_VERSION
  }
#endif
}
	
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;
  char virtual_line[256];
  string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";
  Logger *logger = new Logger(LOG_DEBG);
  logger->log(message);
  chrono::duration<double> elapsed;
  chrono::time_point<chrono::high_resolution_clock> interval_start, interval_end;
  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;
  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;
  int isq, ibf;
  int last_configuration;
  dcomplex ent, entn;
  double enti;
  int num_configs = sconf->configurations;
  int ndirs = sa->nkks;
  int oindex = -1;
  int jindex = jxi488 - output->first_xi + 1;

#ifdef USE_NVTX
  nvtxRangePush("Prepare matrix calculation");
#endif
  double xi = sconf->get_scale(jxi488 - 1);
  double exdc = sconf->exdc;
  double exri = sqrt(exdc);
  int idfc = (int)sconf->idfc;
  double vkarg = 0.0;
  if (idfc >= 0) {
    cid->vk = xi * cid->wn;
    vkarg = cid->vk;
    output->vec_vk[jindex - 1] = cid->vk;
    output->vec_xi[jindex - 1] = xi;
    // goes to 120
  } else { // label 119
    vkarg = xi * cid->vk;
    cid->sqsfi = 1.0 / (xi * xi);
    output->vec_vk[jindex - 1] = cid->vk;
    output->vec_xi[jindex - 1] = xi;
  }
  // label 120
  double sze = vkarg * cid->extr;
  last_configuration = 0;
  for (int i133 = 1; i133 <= cid->c1->nsph; i133++) {
    int iogi = cid->c1->iog[i133 - 1];
    if (iogi != i133) {
      for (int l123 = 1; l123 <= cid->c1->li; l123++) {
	cid->c1->rmi[l123 - 1][i133 - 1] = cid->c1->rmi[l123 - 1][iogi - 1];
	cid->c1->rei[l123 - 1][i133 - 1] = cid->c1->rei[l123 - 1][iogi - 1];
      } // l123 loop
    } else { // label 125
      last_configuration++;
      int nsh = cid->c1->nshl[last_configuration - 1];
      int ici = (nsh + 1) / 2;
      if (i133 == 1) ici++;
      if (idfc == 0) {
	for (int ic = 0; ic < ici; ic++)
	  cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i133 - 1, jxi488 - 1);
	// goes to 129
      } else { // label 127
	if (jxi488 == 1) {
	  for (int ic = 0; ic < ici; ic++) {
	    cid->c1->dc0[ic] = sconf->get_dielectric_constant(ic, i133 - 1, 0);
	  }
	}
      }
      // label 129
      if (i133 == 1) {
	ent = cid->c1->dc0[ici - 1];
	enti = imag(ent);
	entn = csqrt(ent);
	// goes to 131
      } else { // label 130
	if (nsh % 2 == 0) cid->c1->dc0[ici] = ent;
      }
      indme(i133, npnt, npntts, vkarg, ent, enti, entn, jer, lcalc, cid->arg, cid->c1);
      if (jer != 0) {
	output->vec_ier[jindex - 1] = 1;
	output->vec_jxi[jindex - 1] = -jxi488;
	message = "ERROR: indme failed with error code " + to_string(jer) + ".\n";
	logger->log(message, LOG_ERRO);
	delete logger;
	return jer;
	//break;
      }
    }
  } // i133 loop
  ospv(cid->c1, vkarg, sze, exri, entn, enti, jer, lcalc, cid->arg);
  if (jer != 0) {
    output->vec_ier[jindex - 1] = 2;
    output->vec_jxi[jindex - 1] = -jxi488;
    message = "ERROR: ospv failed with error code " + to_string(jer) + ".\n";
    logger->log(message, LOG_ERRO);
    delete logger;
    return jer;
    // break;
  } // i133 loop
#ifdef USE_NVTX
  nvtxRangePop();
#endif
  interval_start = chrono::high_resolution_clock::now();
#ifdef USE_NVTX
  nvtxRangePush("Calculate inverted matrix");
#endif
#ifdef DEBUG_AM
  /* now, before cms, output am to p_outam0 */
  VirtualAsciiFile *outam0 = new VirtualAsciiFile();
  string outam0_name = output_path + "/c_AM0_JXI" + to_string(jxi488) + ".txt";
  sprintf(virtual_line, " AM matrix before CMS\n");
  outam0->append_line(virtual_line);
  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
  outam0->append_line(virtual_line);
  write_dcomplex_matrix(outam0, cid->am, cid->c1->ndm, cid->c1->ndm);
  outam0->write_to_disk(outam0_name);
  delete outam0;
#endif
  incms(cid->am, enti, cid->c1);
#ifdef DEBUG_AM
  VirtualAsciiFile *outam1 = new VirtualAsciiFile();
  string outam1_name = output_path + "/c_AM1_JXI" + to_string(jxi488) + ".txt";
  sprintf(virtual_line, " AM matrix after CMS before LUCIN\n");
  outam1->append_line(virtual_line);
  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
  outam1->append_line(virtual_line);
  write_dcomplex_matrix(outam1, cid->am, cid->c1->ndm, cid->c1->ndm, " %5d %5d (%17.8lE,%17.8lE)\n", 1);
  outam1->write_to_disk(outam1_name);
  delete outam1;
#endif
#ifdef USE_NVTX
  nvtxRangePop();
#endif
  interval_end = chrono::high_resolution_clock::now();
  elapsed = interval_end - interval_start;
  message = "INFO: matrix calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
  logger->log(message);
  interval_start = chrono::high_resolution_clock::now();
#ifdef USE_NVTX
  nvtxRangePush("Invert the matrix");
#endif
  // we the accuracygoal in, get the actual accuracy back out
  double actualaccuracy = cid->accuracygoal;
  invert_matrix(cid->am, cid->c1->ndm, jer, cid->maxrefiters, actualaccuracy, cid->refinemode, output_path, jxi488, mxndm, cid->proc_device);
  // in principle, we should check whether the returned actualaccuracy is indeed lower than the accuracygoal, and do something about it if not
#ifdef USE_REFINEMENT
  if (cid->refinemode==2) {
    message = "INFO: calibration obtained accuracy " + to_string(actualaccuracy) + " (" + to_string(cid->accuracygoal) + " requested) in " + to_string(cid->maxrefiters) + " refinement iterations\n";
    logger->log(message);
    if (actualaccuracy > 1e-2) {
      printf("Accuracy worse than 0.01, stopping");
      exit(1);
    }
  }
#endif // USE_REFINEMENT
#ifdef USE_NVTX
  nvtxRangePop();
#endif
  interval_end = chrono::high_resolution_clock::now();
  elapsed = interval_end - interval_start;
  message = "INFO: matrix inversion for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
  logger->log(message);
  if (jer != 0) {
    message = "ERROR: matrix inversion ended with error code " + to_string(jer) + ".\n";
    logger->err(message);
    delete logger;
    return jer;
    // break; // jxi488 loop: goes to memory clean
  }
  interval_start = chrono::high_resolution_clock::now();
#ifdef USE_NVTX
  nvtxRangePush("Average calculation");
#endif
  exma(cid->am, cid->c1);
#ifdef DEBUG_AM
  VirtualAsciiFile *outam3 = new VirtualAsciiFile();
  string outam3_name = output_path + "/c_AM3_JXI" + to_string(jxi488) + ".txt";
  sprintf(virtual_line, " AM matrix after EXMA\n");
  outam3->append_line(virtual_line);
  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
  outam3->append_line(virtual_line);
  write_dcomplex_matrix(outam3, cid->am, cid->c1->ndm, cid->c1->ndm);
  outam3->write_to_disk(outam3_name);
  delete outam3;
#endif
  if (idfc >= 0) {
    if (jxi488 == jwtm) {
      int nlemt = 2 * cid->c1->nlem;
      string ttms_name = output_path + "/c_TTMS.hd5";
      TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1->am0m, "HDF5");
      ttms_name = output_path + "/c_TTMS";
      TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1->am0m);
    }
  }
  // label 156: continue from here
  last_configuration = 0;
  for (int i168 = 1; i168 <= nsph; i168++) {
    if (cid->c1->iog[i168 - 1] >= i168) {
      int i = i168 - 1;
      last_configuration++;
      oindex = (jindex - 1) * (num_configs + 1) + last_configuration - 1;
      if (cid->c1->nshl[i168 - 1] != 1) {
	output->vec_sphere_ref_indices[oindex] = cc0;
	output->vec_sphere_sizes[oindex] = cid->c1->vsz[i];
      } else {
	output->vec_sphere_ref_indices[oindex] = cid->c1->vkt[i];
	output->vec_sphere_sizes[oindex] = cid->c1->vsz[i];
      }
    } 
  } // i168 loop
  last_configuration++;
  oindex = (jindex - 1) * (num_configs + 1) + last_configuration - 1;
  output->vec_sphere_sizes[oindex] = sze;
  output->vec_sphere_ref_indices[oindex] = entn;
  // label 160
  double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0);
  double csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcs;
  double sqk = cid->vk * cid->vk * exdc;
  vtppoanp->append_line(VirtualBinaryLine(cid->vk));
  pcrsm0(cid->vk, exri, inpol, cid->c1);
  apcra(cid->zpv, cid->c1->le, cid->c1->am0m, inpol, sqk, cid->gapm, cid->gappm);
#ifdef USE_NVTX
  nvtxRangePop();
#endif
  interval_end = chrono::high_resolution_clock::now();
  elapsed = interval_end - interval_start;
  message = "INFO: average calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
  logger->log(message);
  interval_start = chrono::high_resolution_clock::now();
#ifdef USE_NVTX
  nvtxRangePush("Angle loop");
#endif
  double th = sa->th;
  int done_dirs = 0;
  for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable?
    double ph = sa->ph;
    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) {
	upvmp(th, ph, 0, cost, sint, cosp, sinp, cid->u, cid->upmp, cid->unmp);
	if (isam >= 0) {
	  wmamp(
		0, cost, sint, cosp, sinp, inpol, cid->c1->le, 0,
		nsph, cid->argi, cid->u, cid->upmp, cid->unmp, cid->c1
		);
	  // label 182
	  apc(cid->zpv, cid->c1->le, cid->c1->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
	  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
	if (isam >= 0) {
	  // label 182
	  apc(cid->zpv, cid->c1->le, cid->c1->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
	  raba(cid->c1->le, cid->c1->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);
	  jw = 1;
	}
      }
      // label 184
      double thsl = sa->ths;
      double phsph = 0.0;
      for (int jths = 1; jths <= sa->nths; jths++) {
	double ths = thsl;
	int icspnv = 0;
	if (isam > 1) ths += sa->thsca;
	if (isam >= 1) {
	  phsph = 0.0;
	  if (ths < 0.0 || ths > 180.0) phsph = 180.0;
	  if (ths < 0.0) ths *= -1.0;
	  if (ths > 180.0) ths = 360.0 - ths;
	  if (phsph != 0.0) icspnv = 1;
	}
	// label 186
	double phs = sa->phs;
	for (int jphs = 1; jphs <= sa->nphs; jphs++) {
	  double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0;
	  if (isam >= 1) {
	    phs = sa->ph + phsph;
	    if (phs > 360.0) phs -= 360.0;
	  }
	  // label 188
	  bool goto190 = (sa->nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
	  if (!goto190) {
	    upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, cid->us, cid->upsmp, cid->unsmp);
	    if (isam >= 0)
	      wmamp(
		    2, costs, sints, cosps, sinps, inpol, cid->c1->le,
		    0, nsph, cid->args, cid->us, cid->upsmp, cid->unsmp, cid->c1
		    );
	  }
	  // label 190
	  if (sa->nkks != 1 || jxi488 <= 1) {
	    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
		  );
	    if (isam < 0) {
	      wmasp(
		    cost, sint, cosp, sinp, costs, sints, cosps, sinps,
		    cid->u, cid->up, cid->un, cid->us, cid->ups, cid->uns, isq, ibf, inpol, cid->c1->le,
		    0, nsph, cid->argi, cid->args, cid->c1
		    );
	    } else { // label 192
	      for (int i193 = 0; i193 < 3; i193++) {
		cid->up[i193] = cid->upmp[i193];
		cid->un[i193] = cid->unmp[i193];
		cid->ups[i193] = cid->upsmp[i193];
		cid->uns[i193] = cid->unsmp[i193];
	      }
	    }
	  }
	  // label 194
	  if (iavm == 1) crsm1(cid->vk, exri, cid->c1);
	  if (isam < 0) {
	    apc(cid->zpv, cid->c1->le, cid->c1->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);
	    raba(cid->c1->le, cid->c1->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);
	    jw = 1;
	  }
	  // label 196
	  vtppoanp->append_line(VirtualBinaryLine(th));
	  vtppoanp->append_line(VirtualBinaryLine(ph));
	  vtppoanp->append_line(VirtualBinaryLine(ths));
	  vtppoanp->append_line(VirtualBinaryLine(phs));
	  vtppoanp->append_line(VirtualBinaryLine(cid->scan));
	  if (jaw != 0) {
	    jaw = 0;
	    mextc(cid->vk, exri, cid->c1->fsacm, cid->cextlr, cid->cext);
	    // We now have some implicit loops writing to binary
	    for (int i = 0; i < 4; i++) {
	      for (int j = 0; j < 4; j++) {
		double value = cid->cext[i][j];
		vtppoanp->append_line(VirtualBinaryLine(value));
	      }
	    }
	    for (int i = 0; i < 2; i++) {
	      double value = cid->c1->scscm[i];
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = real(cid->c1->scscpm[i]);
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = imag(cid->c1->scscpm[i]);
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = cid->c1->ecscm[i];
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = real(cid->c1->ecscpm[i]);
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = imag(cid->c1->ecscpm[i]);
	      vtppoanp->append_line(VirtualBinaryLine(value));
	    }
	    for (int i = 0; i < 3; i++) {
	      for (int j = 0; j < 2; j++) {
		double value = cid->gapm[i][j];
		vtppoanp->append_line(VirtualBinaryLine(value));
		value = real(cid->gappm[i][j]);
		vtppoanp->append_line(VirtualBinaryLine(value));
		value = imag(cid->gappm[i][j]);
		vtppoanp->append_line(VirtualBinaryLine(value));
	      }
	    }
	    int jlr = 2;
	    for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
	      int ipol = (ilr210 % 2 == 0) ? 1 : -1;
	      if (ilr210 == 2) jlr = 1;
	      double extsm = cid->c1->ecscm[ilr210 - 1];
	      double qextm = extsm * cid->sqsfi / cid->c1->gcs;
	      double scasm = cid->c1->scscm[ilr210 - 1];
	      double albdm = scasm / extsm;
	      double qscam = scasm * cid->sqsfi / cid->c1->gcs;
	      double abssm = extsm - scasm;
	      double qabsm = abssm * cid->sqsfi / cid->c1->gcs;
	      dcomplex s0m = cid->c1->fsacm[ilr210 - 1][ilr210 - 1] * exri;
	      double qschum = imag(s0m) * csch;
	      double pschum = real(s0m) * csch;
	      double s0magm = cabs(s0m) * cs0;
	      // if (inpol == 0) {
	      // sprintf(virtual_line, "   LIN %2d\n", ipol);
	      // output->append_line(virtual_line);
	      // } else { // label 206
	      // sprintf(virtual_line, "  CIRC %2d\n", ipol);
	      // output->append_line(virtual_line);
	      // }
	      // label 208
	      double rapr = cid->c1->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1];
	      double cosav = cid->gapm[2][ilr210 - 1] / cid->c1->scscm[ilr210 - 1];
	      double fz = rapr;
	      if (ipol == -1) {
		output->vec_scs1[jindex - 1] = scasm;
		output->vec_abs1[jindex - 1] = abssm;
		output->vec_exs1[jindex - 1] = extsm;
		output->vec_albeds1[jindex - 1] = albdm;
		output->vec_scsrt1[jindex - 1] = qscam;
		output->vec_absrt1[jindex - 1] = qabsm;
		output->vec_exsrt1[jindex - 1] = qextm;
		output->vec_fsas11[jindex - 1] = cid->c1->fsacm[0][0];
		output->vec_fsas21[jindex - 1] = cid->c1->fsacm[1][0];
		output->vec_qschu1[jindex -1] = qschum;
		output->vec_pschu1[jindex -1] = pschum;
		output->vec_s0mag1[jindex -1] = s0magm;
		output->vec_cosav1[jindex -1] = cosav;
		output->vec_raprs1[jindex -1] = rapr;
		output->vec_fk1[jindex -1] = fz;
	      } else if (ipol == 1) {
		output->vec_scs2[jindex - 1] = scasm;
		output->vec_abs2[jindex - 1] = abssm;
		output->vec_exs2[jindex - 1] = extsm;
		output->vec_albeds2[jindex - 1] = albdm;
		output->vec_scsrt2[jindex - 1] = qscam;
		output->vec_absrt2[jindex - 1] = qabsm;
		output->vec_exsrt2[jindex - 1] = qextm;
		output->vec_fsas22[jindex - 1] = cid->c1->fsacm[1][1];
		output->vec_fsas12[jindex - 1] = cid->c1->fsacm[0][1];
		output->vec_qschu2[jindex -1] = qschum;
		output->vec_pschu2[jindex -1] = pschum;
		output->vec_s0mag2[jindex -1] = s0magm;
		output->vec_cosav2[jindex -1] = cosav;
		output->vec_raprs2[jindex -1] = rapr;
		output->vec_fk2[jindex -1] = fz;
	      }
	    } // ilr210 loop
	  }
	  // label 212
	  output->vec_dir_scand[done_dirs] = cid->scan;
	  output->vec_dir_cfmp[done_dirs] = cid->cfmp;
	  output->vec_dir_sfmp[done_dirs] = cid->sfmp;
	  output->vec_dir_cfsp[done_dirs] = cid->cfsp;
	  output->vec_dir_sfsp[done_dirs] = cid->sfsp;
	  if (isam >= 0) {
	    output->vec_dir_un[3 * done_dirs] = cid->un[0];
	    output->vec_dir_un[3 * done_dirs + 1] = cid->un[1];
	    output->vec_dir_un[3 * done_dirs + 2] = cid->un[2];
	    output->vec_dir_uns[3 * done_dirs] = cid->uns[0];
	    output->vec_dir_uns[3 * done_dirs + 1] = cid->uns[1];
	    output->vec_dir_uns[3 * done_dirs + 2] = cid->uns[2];
	  } else { // label 214
	    output->vec_dir_un[3 * done_dirs] = cid->un[0];
	    output->vec_dir_un[3 * done_dirs + 1] = cid->un[1];
	    output->vec_dir_un[3 * done_dirs + 2] = cid->un[2];
	    output->vec_dir_uns[3 * done_dirs] = 0.0;
	    output->vec_dir_uns[3 * done_dirs + 1] = 0.0;
	    output->vec_dir_uns[3 * done_dirs + 2] = 0.0;
	  }
	  // label 220
	  pcros(cid->vk, exri, cid->c1);
	  mextc(cid->vk, exri, cid->c1->fsac, cid->cextlr, cid->cext);
	  mmulc(cid->c1->vint, cid->cmullr, cid->cmul);
	  if (jw != 0) {
	    jw = 0;
	    // Some implicit loops writing to binary.
	    for (int i = 0; i < 4; i++) {
	      for (int j = 0; j < 4; j++) {
		double value = cid->cext[i][j];
		vtppoanp->append_line(VirtualBinaryLine(value));
	      }
	    }
	    for (int i = 0; i < 2; i++) {
	      double value = cid->c1->scsc[i];
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = real(cid->c1->scscp[i]);
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = imag(cid->c1->scscp[i]);
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = cid->c1->ecsc[i];
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = real(cid->c1->ecscp[i]);
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = imag(cid->c1->ecscp[i]);
	      vtppoanp->append_line(VirtualBinaryLine(value));
	    }
	    for (int i = 0; i < 3; i++) {
	      for (int j = 0; j < 2; j++) {
		double value = cid->gap[i][j];
		vtppoanp->append_line(VirtualBinaryLine(value));
		value = real(cid->gapp[i][j]);
		vtppoanp->append_line(VirtualBinaryLine(value));
		value = imag(cid->gapp[i][j]);
		vtppoanp->append_line(VirtualBinaryLine(value));
	      }
	    }
	    for (int i = 0; i < 2; i++) {
	      for (int j = 0; j < 3; j++) {
		double value = cid->tqce[i][j];
		vtppoanp->append_line(VirtualBinaryLine(value));
		value = real(cid->tqcpe[i][j]);
		vtppoanp->append_line(VirtualBinaryLine(value));
		value = imag(cid->tqcpe[i][j]);
		vtppoanp->append_line(VirtualBinaryLine(value));
	      }
	    }
	    for (int i = 0; i < 2; i++) {
	      for (int j = 0; j < 3; j++) {
		double value = cid->tqcs[i][j];
		vtppoanp->append_line(VirtualBinaryLine(value));
		value = real(cid->tqcps[i][j]);
		vtppoanp->append_line(VirtualBinaryLine(value));
		value = imag(cid->tqcps[i][j]);
		vtppoanp->append_line(VirtualBinaryLine(value));
	      }
	    }
	    for (int i = 0; i < 3; i++) {
	      double value = cid->u[i];
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = cid->up[i];
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = cid->un[i];
	      vtppoanp->append_line(VirtualBinaryLine(value));
	    }
	  }
	  // label 254
	  for (int i = 0; i < 16; i++) {
	    double value = real(cid->c1->vint[i]);
	    vtppoanp->append_line(VirtualBinaryLine(value));
	    value = imag(cid->c1->vint[i]);
	    vtppoanp->append_line(VirtualBinaryLine(value));
	  }
	  for (int i = 0; i < 4; i++) {
	    for (int j = 0; j < 4; j++) {
	      double value = cid->cmul[i][j];
	      vtppoanp->append_line(VirtualBinaryLine(value));
	    }
	  }
	  oindex = (jindex - 1) * ndirs + done_dirs;
	  int jlr = 2;
	  for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
	    int ipol = (ilr290 % 2 == 0) ? 1 : -1;
	    if (ilr290 == 2) jlr = 1;
	    double extsec = cid->c1->ecsc[ilr290 - 1];
	    double qext = extsec * cid->sqsfi / cid->c1->gcs;
	    double scasec = cid->c1->scsc[ilr290 - 1];
	    double albedc = scasec / extsec;
	    double qsca = scasec * cid->sqsfi / cid->c1->gcs;
	    double abssec = extsec - scasec;
	    double qabs = abssec * cid->sqsfi / cid->c1->gcs;
	    dcomplex s0 = cid->c1->fsac[ilr290 - 1][ilr290 - 1] * exri;
	    double qschu = imag(s0) * csch;
	    double pschu = real(s0) * csch;
	    double s0mag = cabs(s0) * cs0;
	    // label 275
	    if (ipol == -1) {
	      output->vec_dir_scs1[oindex] = scasec;
	      output->vec_dir_abs1[oindex] = abssec;
	      output->vec_dir_exs1[oindex] = extsec;
	      output->vec_dir_albeds1[oindex] = albedc;
	      output->vec_dir_scsrt1[oindex] = qsca;
	      output->vec_dir_absrt1[oindex] = qabs;
	      output->vec_dir_exsrt1[oindex] = qext;
	      output->vec_dir_fsas11[oindex] = cid->c1->fsac[0][0];
	      output->vec_dir_fsas21[oindex] = cid->c1->fsac[1][0];
	      output->vec_dir_sas11[oindex] = cid->c1->sac[0][0];
	      output->vec_dir_sas21[oindex] = cid->c1->sac[1][0];
	      output->vec_dir_qschu1[oindex] = qschu;
	      output->vec_dir_pschu1[oindex] = pschu;
	      output->vec_dir_s0mag1[oindex] = s0mag;
	    } else if (ipol == 1) {
	      output->vec_dir_scs2[oindex] = scasec;
	      output->vec_dir_abs2[oindex] = abssec;
	      output->vec_dir_exs2[oindex] = extsec;
	      output->vec_dir_albeds2[oindex] = albedc;
	      output->vec_dir_scsrt2[oindex] = qsca;
	      output->vec_dir_absrt2[oindex] = qabs;
	      output->vec_dir_exsrt2[oindex] = qext;
	      output->vec_dir_fsas22[oindex] = cid->c1->fsac[1][1];
	      output->vec_dir_fsas12[oindex] = cid->c1->fsac[0][1];
	      output->vec_dir_sas22[oindex] = cid->c1->sac[1][1];
	      output->vec_dir_sas12[oindex] = cid->c1->sac[0][1];
	      output->vec_dir_qschu2[oindex] = qschu;
	      output->vec_dir_pschu2[oindex] = pschu;
	      output->vec_dir_s0mag2[oindex] = s0mag;
	    }
	    bool goto290 = isam >= 0 && (jths > 1 || jphs > 1);
	    if (!goto290) {
	      cid->gapv[0] = cid->gap[0][ilr290 - 1];
	      cid->gapv[1] = cid->gap[1][ilr290 - 1];
	      cid->gapv[2] = cid->gap[2][ilr290 - 1];
	      double extins = cid->c1->ecsc[ilr290 - 1];
	      double scatts = cid->c1->scsc[ilr290 - 1];
	      double rapr, cosav, fp, fn, fk, fx, fy, fz;
	      rftr(cid->u, cid->up, cid->un, cid->gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz);
	      if (ipol == -1) {
		output->vec_dir_cosav1[oindex] = cosav;
		output->vec_dir_rapr1[oindex] = rapr;
		output->vec_dir_fl1[oindex] = fp;
		output->vec_dir_fr1[oindex] = fn;
		output->vec_dir_fk1[oindex] = fk;
		output->vec_dir_fx1[oindex] = fx;
		output->vec_dir_fy1[oindex] = fy;
		output->vec_dir_fz1[oindex] = fz;
	      } else if (ipol == 1) {
		output->vec_dir_cosav2[oindex] = cosav;
		output->vec_dir_rapr2[oindex] = rapr;
		output->vec_dir_fl2[oindex] = fp;
		output->vec_dir_fr2[oindex] = fn;
		output->vec_dir_fk2[oindex] = fk;
		output->vec_dir_fx2[oindex] = fx;
		output->vec_dir_fy2[oindex] = fy;
		output->vec_dir_fz2[oindex] = fz;
	      }
	      cid->tqev[0] = cid->tqce[ilr290 - 1][0];
	      cid->tqev[1] = cid->tqce[ilr290 - 1][1];
	      cid->tqev[2] = cid->tqce[ilr290 - 1][2];
	      cid->tqsv[0] = cid->tqcs[ilr290 - 1][0];
	      cid->tqsv[1] = cid->tqcs[ilr290 - 1][1];
	      cid->tqsv[2] = cid->tqcs[ilr290 - 1][2];
	      double tep, ten, tek, tsp, tsn, tsk;
	      tqr(cid->u, cid->up, cid->un, cid->tqev, cid->tqsv, tep, ten, tek, tsp, tsn, tsk);
	      if (ipol == -1) {
		output->vec_dir_tqel1[oindex] = tep;
		output->vec_dir_tqer1[oindex] = ten;
		output->vec_dir_tqek1[oindex] = tek;
		output->vec_dir_tqsl1[oindex] = tsp;
		output->vec_dir_tqsr1[oindex] = tsn;
		output->vec_dir_tqsk1[oindex] = tsk;
		output->vec_dir_tqex1[oindex] = cid->tqce[0][0];
		output->vec_dir_tqey1[oindex] = cid->tqce[0][1];
		output->vec_dir_tqez1[oindex] = cid->tqce[0][2];
		output->vec_dir_tqsx1[oindex] = cid->tqcs[0][0];
		output->vec_dir_tqsy1[oindex] = cid->tqcs[0][1];
		output->vec_dir_tqsz1[oindex] = cid->tqcs[0][2];
	      } else if (ipol == 1) {
		output->vec_dir_tqel2[oindex] = tep;
		output->vec_dir_tqer2[oindex] = ten;
		output->vec_dir_tqek2[oindex] = tek;
		output->vec_dir_tqsl2[oindex] = tsp;
		output->vec_dir_tqsr2[oindex] = tsn;
		output->vec_dir_tqsk2[oindex] = tsk;
		output->vec_dir_tqex2[oindex] = cid->tqce[1][0];
		output->vec_dir_tqey2[oindex] = cid->tqce[1][1];
		output->vec_dir_tqez2[oindex] = cid->tqce[1][2];
		output->vec_dir_tqsx2[oindex] = cid->tqcs[1][0];
		output->vec_dir_tqsy2[oindex] = cid->tqcs[1][1];
		output->vec_dir_tqsz2[oindex] = cid->tqcs[1][2];
	      }
	    }
	  } //ilr290 loop
	  for (int i = 0; i < 4; i++) {
	    oindex = 16 * (jindex - 1) * ndirs + 16 * done_dirs + 4 * i;
	    output->vec_dir_mull[oindex] = cid->cmul[i][0];
	    output->vec_dir_mull[oindex + 1] = cid->cmul[i][1];
	    output->vec_dir_mull[oindex + 2] = cid->cmul[i][2];
	    output->vec_dir_mull[oindex + 3] = cid->cmul[i][3];
	  }
	  for (int i = 0; i < 4; i++) {
	    oindex = 16 * (jindex - 1) * ndirs + 16 * done_dirs + 4 * i;
	    output->vec_dir_mulllr[oindex] = cid->cmullr[i][0];
	    output->vec_dir_mulllr[oindex + 1] = cid->cmullr[i][1];
	    output->vec_dir_mulllr[oindex + 2] = cid->cmullr[i][2];
	    output->vec_dir_mulllr[oindex + 3] = cid->cmullr[i][3];
	  }
	  if (iavm != 0) {
	    mmulc(cid->c1->vintm, cid->cmullr, cid->cmul);
	    // Some implicit loops writing to binary.
	    for (int i = 0; i < 16; i++) {
	      double value;
	      value = real(cid->c1->vintm[i]);
	      vtppoanp->append_line(VirtualBinaryLine(value));
	      value = imag(cid->c1->vintm[i]);
	      vtppoanp->append_line(VirtualBinaryLine(value));
	    }
	    for (int i = 0; i < 4; i++) {
	      for (int j = 0; j < 4; j++) {
		double value = cid->cmul[i][j];
		vtppoanp->append_line(VirtualBinaryLine(value));
	      }
	    }
	    // label 318
	    for (int i = 0; i < 4; i++) {
	      oindex = 16 * (jindex - 1) + 4 * i; // if IAVM fails, try adding directions
	      output->vec_dir_mull[oindex] = cid->cmul[i][0];
	      output->vec_dir_mull[oindex + 1] = cid->cmul[i][1];
	      output->vec_dir_mull[oindex + 2] = cid->cmul[i][2];
	      output->vec_dir_mull[oindex + 3] = cid->cmul[i][3];
	    }
	    for (int i = 0; i < 4; i++) {
	      oindex = 16 * (jindex - 1) + 4 * i; // if IAVM fails, try adding directions
	      output->vec_dir_mulllr[oindex] = cid->cmul[i][0];
	      output->vec_dir_mulllr[oindex + 1] = cid->cmullr[i][1];
	      output->vec_dir_mulllr[oindex + 2] = cid->cmullr[i][2];
	      output->vec_dir_mulllr[oindex + 3] = cid->cmullr[i][3];
	    }
	  }
	  // label 420, continues jphs loop
	  if (isam < 1) phs += sa->phsstp;
	  done_dirs++;
	} // jphs loop, labeled 480
	if (isam <= 1) thsl += sa->thsstp;
      } // jths loop, labeled 482
      ph += sa->phstp;
    } // jph484 loop
    th += sa->thstp;
  } // jth486 loop
#ifdef USE_NVTX
  nvtxRangePop();
#endif
  interval_end = chrono::high_resolution_clock::now();
  elapsed = interval_end - interval_start;
  message = "INFO: angle loop for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
  logger->log(message);
  
  logger->log("INFO: finished scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n");

  delete logger;

  return jer;
}