- begin calling the scales iterations with "fresh" variables, to test parallelism

	fprintf(output, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);

	fprintf(output, " \n");

      }

      for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {

	// 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

      // do the first iteration on jxi488 separately, since it seems to be different from the others

      int jxi488 = 1;

      jer = cluster_jxi488_cycle(jxi488, sconf, gconf, c1, c1ao, c2, c3, c4, c6, c9, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, nth, nths, nph, nphs, nk, nks, nkks, argi, args, gap, gapp, tqce, tqcpe, tqcs, tqcps, duk, tppoan, cextlr, cext, cmullr, cmul, gapv, tqev, tqsv, nxi, nsph, mxndm, inpol, iavm, npnt, npntts, isam, lm, th, thstp, thlst, ths, thsstp, thslst, ph, phstp, phlst, phs, phsstp, phslst, th1, ph1, ths1, phs1, thsca, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, lcalc, arg, wn, vk, ndit, am, isq, ibf);

      for (jxi488 = 2; jxi488 <= nxi; jxi488++) {

	// 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

	ScattererConfiguration *sconf_2 = new ScattererConfiguration(*sconf);

	GeometryConfiguration *gconf_2 = new GeometryConfiguration(*gconf);

	C1 *c1_2 = new C1(*c1);

	C1_AddOns *c1ao_2 = new C1_AddOns(*c1ao);

	C2 *c2_2 = new C2(*c2);

	C3 *c3_2 = new C3(*c3);

	C4 *c4_2 = new C4(*c4);

	C6 *c6_2 = new C6(*c6);

	C9 *c9_2 = new C9(*c9);

	jer = cluster_jxi488_cycle(jxi488, sconf_2, gconf_2, c1_2, c1ao_2, c2_2, c3_2, c4_2, c6_2, c9_2, output, output_path, gaps, tqse, tqspe, tqss, tqsps, zpv, gapm, gappm, nth, nths, nph, nphs, nk, nks, nkks, argi, args, gap, gapp, tqce, tqcpe, tqcs, tqcps, duk, tppoan, cextlr, cext, cmullr, cmul, gapv, tqev, tqsv, nxi, nsph, mxndm, inpol, iavm, npnt, npntts, isam, lm, th, thstp, thlst, ths, thsstp, thslst, ph, phstp, phlst, phs, phsstp, phslst, th1, ph1, ths1, phs1, thsca, u, us, un, uns, up, ups, unmp, unsmp, upmp, upsmp, scan, cfmp, sfmp, cfsp, sfsp, sqsfi, exri, lcalc, arg, wn, vk, ndit, am, isq, ibf);

	delete sconf_2;

	delete gconf_2;

	delete c1_2;

	delete c1ao_2;

	delete c2_2;

	delete c3_2;

	delete c4_2;

	delete c6_2;

	delete c9_2;

      } // jxi488 loop

      tppoan.close();

    } else { // In case TPPOAN could not be opened. Should never happen.

			double sc_ph_start, double sc_ph_step, double sc_ph_end,

			int jwtm

			);

  /*! \brief Build a scattering geometry configuration structure copying it from an existing one.

   *

   * \param rhs: `GeometryConfiguration` preexisting object to copy from.

   */

  GeometryConfiguration(const  GeometryConfiguration& rhs);

  /*! \brief Destroy a GeometryConfiguration instance.

   */

  std::string reference_variable_name;

  //! \brief Number of spherical components.

  int number_of_spheres;

  //! \brief Number of configurations.

  int configurations;

  //! \brief Number of scales to use in calculation.

  int number_of_scales;

  //! \brief Type of dielectric functions (<0 at XIP, =0 as function of XI, >0 contants).

  //! \brief Type of dielectric functions (<0 at XIP, =0 as function of XI, >0 constants).

  int idfc;

  //! \brief External medium dielectric constant. QUESTION: correct?

  double exdc;

   */

  ScattererConfiguration(

			 int nsph,

			 int configs,

			 double *scale_vector,

			 int nxi,

			 std::string variable_name,

			 double wp,

			 double xip

  );

  /*! \brief Build a scatterer configuration structure copying its contents from a preexisting one.

   *

   * Prepare a default configuration structure by allocating the necessary

   * memory structures.

   *

   * \param nsph: `int` The number of spheres in the simulation.

   * \param scale_vector: `double*` The radiation-particle scale vector.

   * \param nxi: `int` The number of radiation-particle scalings.

   * \param variable_name: `string` The name of the radiation-particle scaling type.

   * \param iog_vector: `int*` Array of sphere identification numbers. QUESTION: correct?

   * \param ros_vector: `double*` Sphere radius array.

   * \param nshl_vector: `int*` Array of layer numbers.

   * \param rcf_vector: `double**` Array of fractional break radii. QUESTION: correct?

   * \param dielectric_func_type: `int` Type of dielectric function definition (=0 for constant,

   * \>0 as function of scale parameter, <0 for functions at XIP value and XI is scale factor

   * for dimensions).

   * \param dc_matrix: `complex double ***` Matrix of reference dielectric constants.

   * \param has_external: `bool` Flag to set whether to add an external spherical layer.

   * \param exdc: `double` External medium dielectric constant.

   * \param wp: `double` wp

   * \param xip: `double` xip

   */

  ScattererConfiguration(const  ScattererConfiguration& rhs);

  /*! \brief Destroy a scatterer configuration instance.

   */

  }

  w = new dcomplex*[nlmmt];

  for (int wi = 0; wi < nlmmt; wi++) w[wi] = new dcomplex[4]();

  int configurations = 0;

  configurations = 0;

  for (int ci = 1; ci <= nsph; ci++) {

    if (_iog[ci - 1] >= ci) configurations++;

  }

  delete[] rei;

  for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];

  int conf_index = 0;

  for (int ci = 1; ci <= nsph; ci++) {

    if (iog[ci] >= ci) {

      delete[] rc[conf_index];

      conf_index++;

    }

  for (int ci = 0; ci < configurations; ci++) {

    delete[] rc[ci];

  }

  delete[] rc;

  for (int vi = nsph - 1; vi > - 1; vi--) {

  sph_y = y;

  sph_z = z;

}

GeometryConfiguration::GeometryConfiguration(const  GeometryConfiguration& rhs)

{

  number_of_spheres = rhs.number_of_spheres;

  l_max = rhs.l_max;

  in_pol = rhs.in_pol;

  npnt = rhs.npnt;

  npntts = rhs.npntts;

  meridional_type = rhs.meridional_type;

  li = rhs.li;

  le = rhs.le;

  mxndm = rhs.mxndm;

  iavm = rhs.iavm;

  in_theta_start = rhs.in_theta_start;

  in_theta_step = rhs.in_theta_step;

  in_theta_end = rhs.in_theta_end;

  in_phi_start = rhs.in_phi_start;

  in_phi_step = rhs.in_phi_step ;

  in_phi_end = rhs.in_phi_end;

  sc_theta_start = rhs.sc_theta_start;

  sc_theta_step = rhs.sc_theta_step;

  sc_theta_end = rhs.sc_theta_end;

  sc_phi_start = rhs.sc_phi_start;

  sc_phi_step = rhs.sc_phi_step;

  sc_phi_end = rhs.sc_phi_end;

  jwtm = rhs.jwtm;

  sph_x = new double[number_of_spheres]();

  sph_y = new double[number_of_spheres]();

  sph_z = new double[number_of_spheres]();

  for (int ni=0; ni<number_of_spheres; ni++) {

    sph_x[ni] = rhs.sph_x[ni];

    sph_y[ni] = rhs.sph_y[ni];

    sph_z[ni] = rhs.sph_z[ni];

  }

}

GeometryConfiguration::~GeometryConfiguration() {

  delete[] sph_x;

ScattererConfiguration::ScattererConfiguration(

					       int nsph,

					       int configs, 

					       double *scale_vector,

					       int nxi,

					       string variable_name,

					       double x

) {

  number_of_spheres = nsph;

  configurations = configs;

  number_of_scales = nxi;

  reference_variable_name = variable_name;

  iog_vec = iog_vector;

    }

  }

}

ScattererConfiguration::ScattererConfiguration(const  ScattererConfiguration& rhs)

{

  number_of_spheres = rhs.number_of_spheres;

  configurations = rhs.configurations;

  number_of_scales = rhs.number_of_scales;

  reference_variable_name = rhs.reference_variable_name;

  idfc = rhs.idfc; 

  use_external_sphere = rhs.use_external_sphere;

  exdc = rhs.exdc;

  wp = rhs.wp;

  xip = rhs.xip;

  iog_vec = new int[number_of_spheres]();

  radii_of_spheres = new double[configurations]();

  nshl_vec = new int[configurations]();

  rcf = new double*[configurations];

  scale_vec = new double[number_of_scales]();

  dc0_matrix = new dcomplex**[configurations];

  for (int si = 0; si < number_of_scales; si++) scale_vec[si] = rhs.scale_vec[si];

  for (int si=0; si<number_of_spheres; si++) {

    iog_vec[si] = rhs.iog_vec[si];

  }

  int dim3 = (idfc == 0) ? number_of_scales : 1;

  for (int si=0; si<configurations; si++) {

    radii_of_spheres[si] = rhs.radii_of_spheres[si];

    nshl_vec[si] = rhs.nshl_vec[si];

    rcf[si] = new double[nshl_vec[si]]();

    dc0_matrix[si] = new dcomplex*[number_of_spheres];

    for (int sj=0; sj<nshl_vec[si]; sj++) rcf[si][sj] = rhs.rcf[si][sj];

    for (int sj=0; sj<number_of_spheres; sj++) {

      dc0_matrix[si][sj] = new dcomplex[dim3]();

      for (int sk=0; sk<dim3; sk++) dc0_matrix[si][sj][sk] = rhs.dc0_matrix[si][sj][sk];

    }

  }

}

ScattererConfiguration::~ScattererConfiguration() {

  int configurations = 0;

  ScattererConfiguration *config = new ScattererConfiguration(

							      nsph,

							      configurations,

							      variable_vector,

							      nxi,

							      variable_name,

    status = hdf_file->close();

    conf = new ScattererConfiguration(

				      nsph,

				      configuration_count,

				      xi_vec,

				      nxi,

				      "XIV",

  ScattererConfiguration *conf = new ScattererConfiguration(

							    _nsph,

							    configurations, 

							    _xi_vec,

							    _nxi,

							    "XIV",