Replace get_param configuration functions with read-only parameter views

    exit(1);

  }

  logger->log(" done.\n", LOG_INFO);

  int s_nsph = sconf->get_param("nsph");

  int nsph = gconf->get_param("nsph");

  int s_nsph = sconf->number_of_spheres;

  int nsph = gconf->number_of_spheres;

  if (s_nsph == nsph) {

    // Shortcuts to variables stored in configuration objects

    ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);

    double wp = sconf->get_param("wp");

    double wp = sconf->wp;

    // Global variables for CLU

    C1 *c1 = new C1(gconf, sconf);

    C3 *c3 = new C3();

    logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");

    time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");

    C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);

    double *gaps = new double[(int)gconf->get_param("nsph")]();

    double *gaps = new double[gconf->number_of_spheres]();

    double *tqev = new double[3]();

    double *tqsv = new double[3]();

    double **tqse, **tqss, **tqce, **tqcs;

    tqcs = new double*[2];

    tqcps = new dcomplex*[2];

    for (int ti = 0; ti < 2; ti++) {

      tqse[ti] = new double[(int)gconf->get_param("nsph")]();

      tqspe[ti] = new dcomplex[(int)gconf->get_param("nsph")]();

      tqss[ti] = new double[(int)gconf->get_param("nsph")]();

      tqsps[ti] = new dcomplex[(int)gconf->get_param("nsph")]();

      tqse[ti] = new double[gconf->number_of_spheres]();

      tqspe[ti] = new dcomplex[gconf->number_of_spheres]();

      tqss[ti] = new double[gconf->number_of_spheres]();

      tqsps[ti] = new dcomplex[gconf->number_of_spheres]();

      tqce[ti] = new double[3]();

      tqcpe[ti] = new dcomplex[3]();

      tqcs[ti] = new double[3]();

    // End of global variables for CLU

    fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");

    fprintf(output, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",

	    nsph, c4->li, c4->le, (np_int)gconf->get_param("mxndm"),

	    (int)gconf->get_param("in_pol"), (int)gconf->get_param("npnt"),

	    (int)gconf->get_param("npntts"), (int)gconf->get_param("iavm"),

	    (int)gconf->get_param("meridional_type")

	    nsph, c4->li, c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,

	    gconf->npntts, gconf->iavm, gconf->iavm

	    );

    fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");

    for (int ri = 0; ri < nsph; ri++) fprintf(output, "%17.8lE%17.8lE%17.8lE\n",

	    p_scattering_angles->phsstp, p_scattering_angles->phslst

	    );

    fprintf(output, " READ(IR,*)JWTM\n");

    fprintf(output, " %5d\n", (int)gconf->get_param("jwtm"));

    fprintf(output, " %5d\n", gconf->jwtm);

    fprintf(output, "  READ(ITIN)NSPHT\n");

    fprintf(output, "  READ(ITIN)(IOG(I),I=1,NSPH)\n");

    fprintf(output, "  READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");

      }

    }

    thdps(c4->lm, zpv);

    double exdc = sconf->get_param("exdc");

    double exdc = sconf->exdc;

    double exri = sqrt(exdc);

    double vk = 0.0;

    fprintf(output, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);

#else

      logger->log("INFO: using fall-back lucin() calls.\n", LOG_INFO);

#endif

      int iavm = (int)gconf->get_param("iavm");

      int isam = (int)gconf->get_param("meridional_type");

      int inpol = (int)gconf->get_param("in_pol");

      int nxi = (int)sconf->get_param("nxi");

      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;

      tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));

      tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));

      double wn = wp / 3.0e8;

      double xip = sconf->get_param("xip");

      double xip = sconf->xip;

      double sqsfi = 1.0;

      if (sconf->get_param("idfc") < 0.0) {

      if (sconf->idfc < 0) {

	vk = xip * wn;

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

	fprintf(output, " \n");

int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, C1 *c1, C1_AddOns *c1ao, C2 *c2, C3 *c3, C4 *c4, C6 *c6, C9 *c9, FILE *output, string output_path, double *gaps, double **tqse, dcomplex **tqspe, double **tqss, dcomplex **tqsps, double ****zpv, double **gapm, dcomplex **gappm, double *argi, double *args, double **gap, dcomplex **gapp, double **tqce, dcomplex **tqcpe, double **tqcs, dcomplex **tqcps, double *duk, fstream &tppoan, double **cextlr, double **cext, double **cmullr, double **cmul, double *gapv, double *tqev, double *tqsv, double *u, double *us, double *un, double *uns, double *up, double *ups, double *unmp, double *unsmp, double *upmp, double *upsmp, double &scan, double &cfmp, double &sfmp, double &cfsp, double &sfsp, double sqsfi, dcomplex arg, double wn, double vk, Logger *logger)

{

  int nxi = (int)sconf->get_param("nxi");

  int nxi = sconf->number_of_scales;

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

  int jer = 0;

  int lcalc = 0;

  int jaw = 1;

  int li = (int)gconf->get_param("li");

  int le = (int)gconf->get_param("le");

  int li = gconf->li;

  int le = gconf->le;

  int lm = 0;

  if (le > lm) lm = le;

  if (li > lm) lm = li;

  int nsph = gconf->get_param("nsph");

  int mxndm = gconf->get_param("mxndm");

  int iavm = gconf->get_param("iavm");

  int inpol = gconf->get_param("in_pol");

  int npnt = gconf->get_param("npnt");

  int npntts = gconf->get_param("npntts");

  int isam = gconf->get_param("meridional_type");

  int jwtm = (int)gconf->get_param("jwtm");

  int nsph = gconf->number_of_spheres;

  np_int mxndm = gconf->mxndm;

  int iavm = gconf->iavm;

  int inpol = gconf->in_pol;

  int npnt = gconf->npnt;

  int npntts = gconf->npntts;

  int isam = gconf->iavm;

  int jwtm = gconf->jwtm;

  np_int ndit = 2 * nsph * c4->nlim;

  int isq, ibf;

  fprintf(output, "========== JXI =%3d ====================\n", jxi488);

  double xi = sconf->get_scale(jxi488 - 1);

  double exdc = sconf->get_param("exdc");

  double exdc = sconf->exdc;

  double exri = sqrt(exdc);

  int idfc = (int)sconf->get_param("idfc");

  int idfc = (int)sconf->idfc;

  double vkarg = 0.0;

  if (idfc >= 0) {

    vk = xi * wn;

protected:

  //! \brief Number of spherical components.

  int number_of_spheres;

  int _number_of_spheres;

  //! \brief Maximum field expansion order.

  int l_max;

  int _l_max;

  //! \brief Maximum internal field expansion order.

  int li;

  int _li;

  //! \brief Maximum external field expansion order.

  int le;

  int _le;

  //! \brief Maximum dimension of allocated matrix allowance (deprecated).

  int mxndm;

  np_int _mxndm;

  //! \brief QUESTION: definition?

  int iavm;

  int _iavm;

  //! \brief Incident field polarization status (0 - linear, 1 - circular).

  int in_pol;

  int _in_pol;

  //! \brief Number of transition points. QUESTION: correct?

  int npnt;

  int _npnt;

  //! \brief Transition smoothness. QUESTION: correct?

  int npntts;

  int _npntts;

  //! \brief Type of meridional plane definition.

  int meridional_type;

  int _isam;

  //! \brief Transition matrix layer ID. QUESTION: correct?

  int jwtm;

  int _jwtm;

  //! \brief Incident field initial azimuth.

  double in_theta_start;

  double _in_theta_start;

  //! \brief Incident field azimuth step.

  double in_theta_step;

  double _in_theta_step;

  //! \brief Incident field final azimuth.

  double in_theta_end;

  double _in_theta_end;

  //! \brief Scattered field initial azimuth.

  double sc_theta_start;

  double _sc_theta_start;

  //! \brief Scattered field azimuth step.

  double sc_theta_step;

  double _sc_theta_step;

  //! \brief Scattered field final azimuth.

  double sc_theta_end;

  double _sc_theta_end;

  //! \brief Incident field initial elevation.

  double in_phi_start;

  double _in_phi_start;

  //! \brief Incident field elevation step.

  double in_phi_step;

  double _in_phi_step;

  //! \brief Incident field final elevation.

  double in_phi_end;

  double _in_phi_end;

  //! \brief Scattered field initial elevation.

  double sc_phi_start;

  double _sc_phi_start;

  //! \brief Scattered field elevation step.

  double sc_phi_step;

  double _sc_phi_step;

  //! \brief Scattered field final elevation.

  double sc_phi_end;

  double _sc_phi_end;

  //! \brief Vector of spherical components X coordinates.

  double *sph_x;

  double *_sph_x;

  //! \brief Vector of spherical components Y coordinates.

  double *sph_y;

  double *_sph_y;

  //! \brief Vector of spherical components Z coordinates.

  double *sph_z;

  double *_sph_z;

public:

  //! \brief Read only view on number of spherical components.

  const int& number_of_spheres = _number_of_spheres;

  //! \brief Read only view on maximum field expansion order.

  const int& l_max = _l_max;

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

  const int& li = _li;

  //! \brief Read only view on maximum external field expansion order.

  const int& le = _le;

  //! \brief Read only view on maximum dimension of allocated matrix allowance (deprecated).

  const np_int& mxndm = _mxndm;

  //! \brief QUESTION: definition?

  const int& iavm = _iavm;

  //! \brief Read only view on incident field polarization status (0 - linear, 1 - circular).

  const int& in_pol = _in_pol;

  //! \brief Read only view on number of transition points. QUESTION: correct?

  const int& npnt = _npnt;

  //! \brief Read only view on transition smoothness. QUESTION: correct?

  const int& npntts = _npntts;

  //! \brief Read only view on type of meridional plane definition.

  const int& isam = _isam;

  //! \brief Read only view on transition matrix layer ID. QUESTION: correct?

  const int& jwtm = _jwtm;

  //! \brief Read only view on incident field initial azimuth.

  const double& in_theta_start = _in_theta_start;

  //! \brief Read only view on incident field azimuth step.

  const double& in_theta_step = _in_theta_step;

  //! \brief Read only view on incident field final azimuth.

  const double& in_theta_end = _in_theta_end;

  //! \brief Read only view on scattered field initial azimuth.

  const double& sc_theta_start = _sc_theta_start;

  //! \brief Read only view on scattered field azimuth step.

  const double& sc_theta_step = _sc_theta_step;

  //! \brief Read only view on scattered field final azimuth.

  const double& sc_theta_end = _sc_theta_end;

  //! \brief Read only view on incident field initial elevation.

  const double& in_phi_start = _in_phi_start;

  //! \brief Read only view on incident field elevation step.

  const double& in_phi_step = _in_phi_step;

  //! \brief Read only view on incident field final elevation.

  const double& in_phi_end = _in_phi_end;

  //! \brief Read only view on scattered field initial elevation.

  const double& sc_phi_start = _sc_phi_start;

  //! \brief Read only view on scattered field elevation step.

  const double& sc_phi_step = _sc_phi_step;

  //! \brief Read only view on scattered field final elevation.

  const double& sc_phi_end = _sc_phi_end;

  /*

  //! \brief Read only view on vector of spherical components X coordinates.

  const double *sph_x = _sph_x;

  //! \brief Read only view on vector of spherical components Y coordinates.

  const double *sph_y = _sph_y;

  //! \brief Read only view on vector of spherical components Z coordinates.

  const double *sph_z = _sph_z;

  */

  /*! \brief Build a scattering geometry configuration structure.

   *

   * \param nsph: `int` Number of spheres to be used in calculation.

   */

  GeometryConfiguration(

			int nsph, int lm, int in_pol, int npnt, int npntts, int meridional_type,

			int li, int le, int mxndm, int iavm,

			int li, int le, np_int mxndm, int iavm,

			double *x, double *y, double *z,

			double in_th_start, double in_th_step, double in_th_end,

			double sc_th_start, double sc_th_step, double sc_th_end,

			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.

   */

  static GeometryConfiguration *from_legacy(const std::string& file_name);

  /*! \brief Get the value of a parameter by name.

   *

   * The proper way to access read-only parameters from outside a class is to define

   * public methods that return their values. For configuration operations, whose

   * optimization is not critical, it is possible to define a single function that

   * returns simple scalar values called by name. Access to more complicated data

   * structures, on the other hand, require specialized methods which avoid the

   * burden of searching the necessary value across the whole array every time.

   *

   * \param param_name: `string` Name of the parameter to be retrieved.

   * \return value: `double` Value of the requested parameter.

   */

  double get_param(const std::string& param_name);

  /*! \brief Get the X coordinate of a sphere by its index.

   *

   * This is a specialized function to access the X coordinate of a sphere through

   * \param index: `int` Index of the scale to be retrieved.

   * \return scale: `double` The X coordinate of the requested sphere.

   */

  double get_sph_x(int index) { return sph_x[index]; }

  double get_sph_x(int index) { return _sph_x[index]; }

  /*! \brief Get the Y coordinate of a sphere by its index.

   *

   * \param index: `int` Index of the scale to be retrieved.

   * \return scale: `double` The Y coordinate of the requested sphere.

   */

  double get_sph_y(int index) { return sph_y[index]; }

  double get_sph_y(int index) { return _sph_y[index]; }

  /*! \brief Get the Z coordinate of a sphere by its index.

   *

   * \param index: `int` Index of the scale to be retrieved.

   * \return scale: `double` The Z coordinate of the requested sphere.

   */

  double get_sph_z(int index) { return sph_z[index]; }

  double get_sph_z(int index) { return _sph_z[index]; }

};

/**

protected:

  //! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].

  dcomplex ***dc0_matrix;

  dcomplex ***_dc0_matrix;

  //! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].

  double *radii_of_spheres;

  double *_radii_of_spheres;

  //! \brief Matrix of fractional transition radii with size [N_SPHERES x LAYERS].

  double **rcf;

  double **_rcf;

  //! \brief Vector of sphere ID numbers, with size [N_SPHERES].

  int *iog_vec;

  int *_iog_vec;

  //! \brief Vector of layer numbers for every sphere, with size [N_SPHERES].

  int *nshl_vec;

  int *_nshl_vec;

  //! \brief Vector of scale parameters, with size [N_SCALES].

  double *scale_vec;

  double *_scale_vec;

  //! \brief Name of the reference variable type (one of XIV, WNS, WLS, PUS, EVS).

  std::string reference_variable_name;

  std::string _reference_variable_name;

  //! \brief Number of spherical components.

  int number_of_spheres;

  int _number_of_spheres;

  //! \brief Number of configurations.

  int configurations;

  int _configurations;

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

  int number_of_scales;

  int _number_of_scales;

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

  int idfc;

  int _idfc;

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

  double exdc;

  double _exdc;

  //! \brief WP. QUESTION: better definition?

  double wp;

  double _wp;

  //! \brief Peak XI. QUESTION: correct?

  double xip;

  double _xip;

  //! \brief Flag to control whether to add an external layer.

  bool use_external_sphere;

  bool _use_external_sphere;

  /*! \brief Build configuration from a HDF5 binary input file.

   *

   */

  void write_legacy(const std::string& file_name);

public:

  //! \brief Read only view on name of the reference variable type.

  const std::string& reference_variable_name = _reference_variable_name;

  //! \brief Read only view on number of spherical components.

  const int& number_of_spheres = _number_of_spheres;

  //! \brief Read only view on number of configurations.

  const int& configurations = _configurations;

  //! \brief Read only view on number of scales to use in calculation.

  const int& number_of_scales = _number_of_scales;

  //! \brief Read only view on type of dielectric functions.

  const int& idfc = _idfc;

  //! \brief Read only view on external medium dielectric constant.

  const double& exdc = _exdc;

  //! \brief Read only view on WP.

  const double& wp = _wp;

  //! \brief Read only view on peak XI.

  const double& xip = _xip;

  //! \brief Read only view on flag to control whether to add an external layer.

  const bool& use_external_sphere = _use_external_sphere;

  /*! \brief Build a scatterer configuration structure.

   *

   * Prepare a default configuration structure by allocating the necessary

			 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

   * \param rhs: `ScattererConfiguration&` Reference to the ScattereConfiguration

   * object to be copied.

   */

  ScattererConfiguration(const ScattererConfiguration& rhs);

   * \param k: `int` Index of the current scale.

   * \return radius: `dcomplex` The requested dielectric constant.

   */

  dcomplex get_dielectric_constant(int i, int j, int k) { return dc0_matrix[i][j][k]; }

  dcomplex get_dielectric_constant(int i, int j, int k) { return _dc0_matrix[i][j][k]; }

  /*! \brief Get the ID of a configuration from the index of the sphere.

   *

   * \param index: `int` Index of the sphere.

   * \return ID: `int` ID of the configuration to be applied.

   */

  int get_iog(int index) { return iog_vec[index]; }

  int get_iog(int index) { return _iog_vec[index]; }

  /*! \brief Get the number of layers for a given configuration.

   *

   * \param index: `int` Index of the configuration.

   * \return nl: `int` The number of layers for the given configuration.

   */

  int get_nshl(int index) { return nshl_vec[index]; }

  int get_nshl(int index) { return _nshl_vec[index]; }

  /*

  /*! \brief Get the value of a parameter by name.

   *

   * The proper way to access read-only parameters from outside a class is to define

   *

   * \param param_name: `string` Name of the parameter to be retrieved.

   * \return value: `double` Value of the requested parameter.

   */

  double get_param(const std::string& param_name);

  */

  /*! \brief Get the radius of a sphere by its index.

   *

   * \param index: `int` Index of the ID to be retrieved.

   * \return radius: `double` The requested sphere radius.

   */

  double get_radius(int index) { return radii_of_spheres[index]; }

  double get_radius(int index) { return _radii_of_spheres[index]; }

  /*! \brief Get the value of a scale by its index.

   *

   * \param index: `int` Index of the scale to be retrieved.

   * \return scale: `double` The desired scale.

   */

  double get_rcf(int row, int column) { return rcf[row][column]; }

  /*! \brief Get the reference variable name.

   *

   * This is a specialized function to get the name of the reference variable as a

   * string.

   *

   * \return name: `string` The name of the variable used to calculate wavelength /

   * size scaling.

   */

  std::string get_reference() { return reference_variable_name; }

  double get_rcf(int row, int column) { return _rcf[row][column]; }

  /*! \brief Get the value of a scale by its index.

   *

   * \param index: `int` Index of the scale to be retrieved.

   * \return scale: `double` The desired scale.

   */

  double get_scale(int index) { return scale_vec[index]; }

  double get_scale(int index) { return _scale_vec[index]; }

  /*! \brief Print the contents of the configuration object to terminal.

   *

#endif

C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {

  lm = (int)gconf->get_param("l_max");

  int li = (int)gconf->get_param("li");

  int le = (int)gconf->get_param("le");

  lm = gconf->l_max;

  int li = gconf->li;

  int le = gconf->le;

  if (lm == 0) {

    lm = (li > le) ? li : le;

  }

  nsph = (int)gconf->get_param("nsph");

  nsph = gconf->number_of_spheres;

  nlmmt = 2 * (lm * (lm + 2));

  vec_rmi = new dcomplex[nsph * lm]();

  vec_w = new dcomplex[4 * nlmmt]();

  w = new dcomplex*[nlmmt];

  for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);

  configurations = (int)sconf->get_param("configurations");

  configurations = sconf->configurations;

  vint = new dcomplex[16]();

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

  vints = new dcomplex*[nsph];

}

C2::C2(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {

  nsph = (int)gconf->get_param("nsph");

  int npnt = (int)gconf->get_param("npnt");

  int npntts = (int)gconf->get_param("npntts");

  nsph = gconf->number_of_spheres;

  int npnt = gconf->npnt;

  int npntts = gconf->npntts;

  int max_n = (npnt > npntts) ? npnt : npntts;

  nhspo = 2 * max_n - 1;

  nl = (int)sconf->get_param("configurations");

  nl = sconf->configurations;

  if (nsph == 1 && nl == 1) nl = 5;

  ris = new dcomplex[nhspo]();

  dlri = new dcomplex[nhspo]();

}

C4::C4(GeometryConfiguration *gconf) {

  li = (int)gconf->get_param("li");

  le = (int)gconf->get_param("le");