src/cluster/np_cluster.cpp

@@ -31,7 +31,7 @@
 
 using namespace std;
 
-extern void cluster(string config_file, string data_file, string output_path);
+extern void cluster(const string& config_file, const string& data_file, const string& output_path);
 
 /*! \brief Main program entry point.
  *

src/include/Commons.h

@@ -99,12 +99,10 @@ public:
 
   /*! \brief C1 instance constructor.
    *
-   * \param ns: `int` Number of spheres.
-   * \param l_max: `int` Maximum order of field expansion.
-   * \param nshl: `int *` Array of number of layers in spheres.
-   * \param iog: `int *` Vector of spherical units ID numbers.
+   * \param gconf: `GeometryConfiguration *` Pointer to a geometry configuration object.
+   * \param sconf: `ScattererConfiguration *` Pointer to a scatterer configuration object.
    */
-  C1(int ns, int l_max, int *nshl, int *iog);
+  C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf);
 
   /*! \brief C1 instance constructor copying all contents from a preexisting template
    *
@@ -142,12 +140,10 @@ public:
 
   /*! \brief C2 instance constructor.
    *
-   * \param ns: `int` Number of spheres.
-   * \param nl: `int`
-   * \param npnt: `int`
-   * \param npntts: `int`
+   * \param gconf: `GeometryConfiguration*` Pointer to a GeometryConfiguration instance.
+   * \param sconf: `ScattererConfiguration*` Pointer to a ScattererConfiguration instance.
    */
-  C2(int ns, int nl, int npnt, int npntts);
+  C2(GeometryConfiguration *gconf, ScattererConfiguration *sconf);
 
   /*! \brief C2 instance constructor copying its contents from preexisting instance.
    *
@@ -193,39 +189,44 @@ public:
  */
 class C4 {
 public:
-  //! \brief QUESTION: definition?
+  //! \brief LITPO = 2 * LI + 1.
   int litpo;
-  //! \brief QUESTION: definition?
+  //! \brief LITPOS = LITPO * LITPO
   int litpos;
-  //! \brief Maximum field expansion order plus one. QUESTION: correct?
+  //! \brief LMPO = LM + 1.
   int lmpo;
-  //! \brief Twice maximum field expansion order plus one. QUESTION: correct?
+  //! \brief LMTPO = 2 * LM + 1.
   int lmtpo;
-  //! \brief Square of `lmtpo`.
+  //! \brief LMTPOS = LMTPO * LMTPO.
   int lmtpos;
-  //! \brief QUESTION: definition?
+  //! \brief Internal field expansion order.
   int li;
   //! \brief QUESTION: definition?
   int nlim;
-  //! \brief QUESTION: definition?
+  //! \brief External field expansion order.
   int le;
   //! \brief QUESTION: definition?
   int nlem;
-  //! \brief Maximum field expansion order. QUESTION: correct?
+  //! \brief Maximum field expansion order.
   int lm;
   //! \brief Number of spheres.
   int nsph;
   //! \brief QUESTION: definition?
   int nv3j;
 
-  /*! \brief C3 instance constructor.
+  /*! \brief C4 instance constructor.
+   *
+   * \param gconf: `GeometryConfiguration*` Pointer to a GeometryConfiguration instance.
    */
-  C4(int li, int le, int nsph);
-  /*! \brief C3 instance constructor copying its contents from a preexisting object.
+  C4(GeometryConfiguration *gconf);
+  
+  /*! \brief C4 instance constructor copying its contents from a preexisting object.
+   *
+   * \param rhs: `C4&` Reference of the object to be copied.
    */
   C4(const C4& rhs);
 
-  /*! \brief C3 instance destroyer.
+  /*! \brief C4 instance destroyer.
    */
   ~C4();
 };
@@ -268,6 +269,8 @@ public:
   //! \brief QUESTION: definition?
   dcomplex **am0m;
   //! \brief QUESTION: definition?
+  dcomplex *am0v;
+  //! \brief QUESTION: definition?
   dcomplex *vintm;
   //! \brief QUESTION: definition?
   dcomplex *vintt;
@@ -382,4 +385,183 @@ public:
    */
   ~C9();
 };
+
+/*! \brief A data structure representing the information used for a single scale
+ * of the CLUSTER case.
+ */
+class ClusterIterationData {
+public:
+  //! \brief Pointer to a C1 structure.
+  C1 *c1;
+  //! \brief Pointer to a C1_AddOns structure.
+  C1_AddOns *c1ao;
+  //! \brief Pointer to a C2 structure.
+  C2 *c2;
+  //! \brief Pointer to a C3 structure.
+  C3 *c3;
+  //! brief Pointer to a C4 structure.
+  C4 *c4;
+  //! \brief Pointer to a C6 structure.
+  C6 *c6;
+  //! \brief Pointer to a C9 structure.
+  C9 *c9;
+  //! \brief Pointer to a formatted output file.
+  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;
+  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 qsfi;
+  double sqsfi;
+  dcomplex *am_vector;
+  dcomplex **am;
+  dcomplex arg;
+  //! \brief Wave vector.
+  double vk;
+  //! \brief Wave number.
+  double wn;
+  double xip;
+
+  ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf);
+
+  ClusterIterationData(const ClusterIterationData& rhs);
+
+  ~ClusterIterationData();
+};
+
+/*! \brief A data structure representing the angles to be evaluated in the problem.
+ *
+ */
+class ScatteringAngles {
+protected:
+  //! \brief Number of incident field azimuth angles.
+  int _nth;
+  //! \brief Number of scattered field azimuth angles.
+  int _nths;
+  //! \brief Number of incident field elevation angles.
+  int _nph;
+  //! \brief Number of scattered field elevation angles.
+  int _nphs;
+  //! \brief Number of incident field propagation angles.
+  int _nk;
+  //! \brief Number of scattered field propagation angles.
+  int _nks;
+  //! \brief Total number of field propagation angles.
+  int _nkks;
+  //! \brief First incident field azimuth angle.
+  double _th;
+  //! \brief Incident field azimuth angle increment.
+  double _thstp;
+  //! \brief Last incident field azimuth angle.
+  double _thlst;
+  //! \brief First scattered field azimuth angle.
+  double _ths;
+  //! \brief Scattered field azimuth angle increment.
+  double _thsstp;
+  //! \brief Last scattered field azimuth angle.
+  double _thslst;
+  //! \brief First incident field elevation angle.
+  double _ph;
+  //! \brief Incident field elevation angle increment.
+  double _phstp;
+  //! \brief Last incident field elevation angle.
+  double _phlst;
+  //! \brief First scattered field elevation angle.
+  double _phs;
+  //! \brief Scattered field elevation angle increment.
+  double _phsstp;
+  //! \brief Last scattered field elevation angle.
+  double _phslst;
+  //! \brief Azimuth scattering deflection.
+  double _thsca;
+
+public:
+  //! \brief Read only view of `_nth`.
+  const int& nth = _nth;
+  //! \brief Read only view of `_nths`.
+  const int& nths = _nths;
+  //! \brief Read only view of `_nph`.
+  const int& nph = _nph;
+  //! \brief Read only view of `_nphs`.
+  const int& nphs = _nphs;
+  //! \brief Read only view of `_nk`.
+  const int& nk = _nk;
+  //! \brief Read only view of `_nks`.
+  const int& nks = _nks;
+  //! \brief Read only view of `_nkks`.
+  const int& nkks = _nkks;
+  //! \brief Read only view of `_th`.
+  const double& th = _th;
+  //! \brief Read only view of `_thstp`.
+  const double& thstp = _thstp;
+  //! \brief Read only view of `_thlst`.
+  const double& thlst = _thlst;
+  //! \brief Read only view of `_ths`.
+  const double& ths = _ths;
+  //! \brief Read only view of `_thsstp`.
+  const double& thsstp = _thsstp;
+  //! \brief Read only view of `_thslst`.
+  const double& thslst = _thslst;
+  //! \brief Read only view of `_ph`.
+  const double& ph = _ph;
+  //! \brief Read only view of `_phstp`.
+  const double& phstp = _phstp;
+  //! \brief Read only view of `_phlst`.
+  const double& phlst = _phlst;
+  //! \brief Read only view of `_phs`.
+  const double& phs = _phs;
+  //! \brief Read only view of `_phsstp`.
+  const double& phsstp = _phsstp;
+  //! \brief Read only view of `_phslst`.
+  const double& phslst = _phslst;
+  //! \brief Read only view of `_thsca`.
+  const double& thsca = _thsca;
+
+  /*! \brief ScatteringAngles instance constructor.
+   *
+   * \param gconf: `GeometryConfiguration*` Pointer to a GeometryConfiguration object.
+   */
+  ScatteringAngles(GeometryConfiguration *gconf);
+  
+  /*! \brief ScatteringAngles copy constructor.
+   *
+   * \param rhs: `ScatteringAngles&` Reference to the ScatteringAngles object to be copied.
+   */
+  ScatteringAngles(const ScatteringAngles &rhs);
+};
+
 #endif

src/include/Parsers.h

@@ -32,6 +32,6 @@
  * \return array_lines `string*` An array of strings, one for each input
  * file line.
  */
-std::string *load_file(std::string file_name, int *count);
+std::string *load_file(const std::string& file_name, int *count);
 
 #endif /* INCLUDE_PARSERS_H_ */

src/include/TransitionMatrix.h

@@ -36,7 +36,7 @@ class TransitionMatrix {
    * \return config: `TransitionMatrix *` Pointer to object containing the
    * transition matrix data.
    */
-  static TransitionMatrix *from_hdf5(std::string file_name);
+  static TransitionMatrix *from_hdf5(const std::string& file_name);
   
   /*! \brief Build transition matrix from a legacy binary input file.
    *
@@ -44,23 +44,102 @@ class TransitionMatrix {
    * \return config: `TransitionMatrix *` Pointer to object containing the
    * transition matrix data.
    */
-  static TransitionMatrix *from_legacy(std::string file_name);
+  static TransitionMatrix *from_legacy(const std::string& file_name);
 
   /*! \brief Write the Transition Matrix to HDF5 binary output.
    *
-   * This function takes care of the specific task of building a transition
-   * matrix memory data structure from a binary input file formatted according
-   * to the structure used by the original FORTRAN code.
+   * This function takes care of the specific task of writing the transition
+   * matrix memory data structure to a binary output file formatted according
+   * to the HDF5 standard.
    *
    * \param file_name: `string` Name of the binary configuration data file.
    */
-  void write_hdf5(std::string file_name);
+  void write_hdf5(const std::string& file_name);
+  
+  /*! \brief Write transition matrix data to HDF5 binary output.
+   *
+   * This function takes care of the specific task of writing the transition
+   * matrix memory data structure to a binary output file formatted according
+   * to the HDF5 standard without a pre-existing instance. It is designed to
+   * work for the case of a cluster of spheres.
+   *
+   * \param file_name: `string` Name of the binary configuration data file.
+   * \param _nlemt: `np_int` Size of the matrix (2 * LE * (LE + 2)).
+   * \param _lm: `int` Maximum field expansion order.
+   * \param _vk: `double` Wave number in scale units.
+   * \param _exri: `double` External medium refractive index.
+   * \param _am0m: `complex double **`
+   */
+  static void write_hdf5(
+			 const std::string& file_name, np_int _nlemt, int _lm, double _vk,
+			 double _exri, dcomplex **_am0m
+  );
   
+  /*! \brief Write transition matrix data to HDF5 binary output.
+   *
+   * This function takes care of the specific task of writing the transition
+   * matrix memory data structure to a binary output file formatted according
+   * to the HDF5 standard without a pre-existing instance. It is designed to
+   * work for the case of a single sphere.
+   *
+   * \param file_name: `string` Name of the binary configuration data file.
+   * \param _lm: `int` Maximum field expansion order.
+   * \param _vk: `double` Wave number in scale units.
+   * \param _exri: `double` External medium refractive index.
+   * \param _rmi: `complex double **`
+   * \param _rei: `complex double **`
+   * \param _sphere_radius: `double` Radius of the sphere.
+   */
+  static void write_hdf5(
+			 const std::string& file_name, int _lm, double _vk, double _exri,
+			 dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
+  );
+
   /*! \brief Write the Transition Matrix to legacy binary output.
    *
    * \param file_name: `string` Name of the binary configuration data file.
    */
-  void write_legacy(std::string file_name);
+  void write_legacy(const std::string& file_name);
+
+  /*! \brief Write transition matrix data to binary output using legacy format.
+   *
+   * This function takes care of the specific task of writing the transition
+   * matrix memory data structure to a binary output file formatted according
+   * to the format used by the legacy FORTRAN code. It is designed to work for
+   * the case of clusters of spheres.
+   *
+   * \param file_name: `string` Name of the binary configuration data file.
+   * \param _nlemt: `np_int` Size of the matrix (2 * LE * (LE + 2)).
+   * \param _lm: `int` Maximum field expansion order.
+   * \param _vk: `double` Wave number in scale units.
+   * \param _exri: `double` External medium refractive index.
+   * \param _am0m: `complex double **`
+   */
+  static void write_legacy(
+			   const std::string& file_name, np_int _nlemt, int _lm, double _vk,
+			   double _exri, dcomplex **_am0m
+  );
+  
+  /*! \brief Write transition matrix data to binary output using legacy format.
+   *
+   * This function takes care of the specific task of writing the transition
+   * matrix memory data structure to a binary output file formatted according
+   * to the original code structure without a pre-existing instance. It is designed
+   * to work for the case of a single sphere.
+   *
+   * \param file_name: `string` Name of the binary configuration data file.
+   * \param _lm: `int` Maximum field expansion order.
+   * \param _vk: `double` Wave number in scale units.
+   * \param _exri: `double` External medium refractive index.
+   * \param _rmi: `complex double **`
+   * \param _rei: `complex double **`
+   * \param _sphere_radius: `double` Radius of the sphere.
+   */
+  static void write_legacy(
+			   const std::string& file_name, int _lm, double _vk, double _exri,
+			   dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
+  );
+
  public:
   /*! \brief Default Transition Matrix instance constructor.
    *
@@ -126,7 +205,7 @@ class TransitionMatrix {
    * \return config: `TransitionMatrix *` Pointer to object containing the transition
    * matrix data.
    */
-  static TransitionMatrix* from_binary(std::string file_name, std::string mode="LEGACY");
+  static TransitionMatrix* from_binary(const std::string& file_name, const std::string& mode="LEGACY");
 
   /*! \brief Write the Transition Matrix to a binary file.
    *
@@ -140,8 +219,64 @@ class TransitionMatrix {
    * \param mode: `string` Binary encoding. Can be one of ["LEGACY", "HDF5"] . Optional
    * (default is "LEGACY").
    */
-  void write_binary(std::string file_name, std::string mode="LEGACY");
+  void write_binary(const std::string& file_name, const std::string& mode="LEGACY");
   
+  /*! \brief Write a cluster Transition Matrix to a binary file without instanciating it.
+   *
+   * Transition Matrix data can take a large amount of memory. For such reason, attempts
+   * to create TransitionMatrix instances only for writing purposes can create
+   * unnecessary resource consumption and computing time to duplicate the data into
+   * the output buffer. This function offers output to file as a static method. It
+   * takes the arguments of a constructor together with the usual arguments to specify
+   * the output file name and format, to write the required data directly to a file,
+   * without creating a new TransitionMatrix instance. The implementation works for
+   * TransitionMatrix objects built for the CLUSTER case. It belongs to the public class
+   * interface and it calls the proper versions of `write_legacy()` and `write_hdf5()`,
+   * depending on the requested output format.
+   * 
+   * \param file_name: `string` Name of the file to be written.
+   * \param _nlemt: `np_int` Size of the matrix (2 * LE * (LE + 2)).
+   * \param _lm: `int` Maximum field expansion order.
+   * \param _vk: `double` Wave number in scale units.
+   * \param _exri: `double` External medium refractive index.
+   * \param _am0m: `complex double **`
+   * \param mode: `string` Binary encoding. Can be one of ["LEGACY", "HDF5"] . Optional
+   * (default is "LEGACY").
+   */
+  static void write_binary(
+			   const std::string& file_name, np_int _nlemt, int _lm, double _vk,
+			   double _exri, dcomplex **_am0m, const std::string& mode="LEGACY"
+  );
+  
+  /*! \brief Write a single sphere Transition Matrix to a binary file without instanciating it.
+   *
+   * Transition Matrix data can take a large amount of memory. For such reason, attempts
+   * to create TransitionMatrix instances only for writing purposes can create
+   * unnecessary resource consumption and computing time to duplicate the data into
+   * the output buffer. This function offers output to file as a static method. It
+   * takes the arguments of a constructor together with the usual arguments to specify
+   * the output file name and format, to write the required data directly to a file,
+   * without creating a new TransitionMatrix instance. The implementation works for
+   * TransitionMatrix objects built for the single sphere case. It belongs to the public
+   * class interface and it calls the proper versions of `write_legacy()` and `write_hdf5()`,
+   * depending on the requested output format.
+   *
+   * \param file_name: `string` Name of the file to be written.
+   * \param _lm: `int` Maximum field expansion order.
+   * \param _vk: `double` Wave number in scale units.
+   * \param _exri: `double` External medium refractive index.
+   * \param _rmi: `complex double **`
+   * \param _rei: `complex double **`
+   * \param _sphere_radius: `double` Radius of the sphere.
+   * \param mode: `string` Binary encoding. Can be one of ["LEGACY", "HDF5"] . Optional
+   * (default is "LEGACY").
+   */
+  static void write_binary(
+			   const std::string& file_name, int _lm, double _vk, double _exri,
+			   dcomplex **_rmi, dcomplex **_rei, double _sphere_radius,
+			   const std::string& mode="LEGACY"
+  );
+
   /*! \brief Test whether two instances of TransitionMatrix are equal.
    *
    * Transition matrices can be the result of a calculation or of a file input operation,

src/include/clu_subs.h

@@ -344,14 +344,14 @@ void scr2(
  * to radial coordinates, then it calls `sphar()` to calculate the vector of spherical
  * harmonics of the incident field.
  *
- * \param rcf: `double **` Matrix of double.
+ * \param sconf: `ScattererConfiguration *` Pointer to scatterer configuration object.
  * \param c1: `C1 *` Pointer to a C1 instance.
  * \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.
  * \param c3: `C3 *` Pointer to a C3 instance.
  * \param c4: `C4 *` Pointer to a C4 structure.
  * \param c6: `C6 *` Pointer to a C6 instance.
  */
-void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6);
+void str(ScattererConfiguration *sconf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6);
 
 /*! \brief Compute radiation torques on particles on a k-vector oriented system.
  *

src/include/errors.h

@@ -64,7 +64,7 @@ public:
    *
    * \param name: `string` Name of the file that was accessed.
    */
-  OpenConfigurationFileException(std::string name) { file_name = name; }
+  OpenConfigurationFileException(const std::string& name) { file_name = name; }
 
   /**
    * \brief Exception message.
@@ -86,7 +86,7 @@ public:
    *
    * \param problem: `string` Description of the problem that occurred.
    */
-  MatrixOutOfBoundsException(std::string problem) { message = problem; }
+  MatrixOutOfBoundsException(const std::string& problem) { message = problem; }
   /**
    * \brief Exception message.
    */
@@ -107,7 +107,7 @@ public:
    *
    * \param problem: `string` Description of the problem that occurred.
    */
-  UnrecognizedConfigurationException(std::string problem) { message = problem; }
+  UnrecognizedConfigurationException(const std::string& problem) { message = problem; }
   /**
    * \brief Exception message.
    */
@@ -128,7 +128,7 @@ public:
    *
    * \param problem: `string` Description of the problem that occurred.
    */
-  UnrecognizedFormatException(std::string problem) { message = problem; }
+  UnrecognizedFormatException(const std::string& problem) { message = problem; }
   /**
    * \brief Exception message.
    */
@@ -149,7 +149,7 @@ public:
    *
    * \param problem: `string` Description of the problem that occurred.
    */
-  UnrecognizedParameterException(std::string problem) { message = problem; }
+  UnrecognizedParameterException(const std::string& problem) { message = problem; }
   /**
    * \brief Exception message.
    */

src/include/file_io.h

@@ -31,7 +31,7 @@ class FileSchema {
    * \param rec_types: `string *` Description of the records in the file.
    * \param rec_names: `string *` Names of the records in the file.
    */
-  FileSchema(int num_rec, std::string *rec_types, std::string *rec_names=NULL);
+  FileSchema(int num_rec, const std::string *rec_types, const std::string *rec_names=NULL);
 
   /*! \brief FileSchema instance destroyer.
    */
@@ -43,12 +43,13 @@ class FileSchema {
    */
   int get_record_number() { return num_records; }
 
-  /*! \brief Get a copy of the record types.
+  /*! \brief Get a copy of the record names.
    *
    * \return rec_types: `string *` A new vector of strings with description of records.
    */
   std::string *get_record_names();
-  /*! \brief Get a copy of the record names.
+  
+  /*! \brief Get a copy of the record types.
    *
    * \return rec_names: `string *` A new vector of strings with record names.
    */
@@ -83,7 +84,7 @@ class HDFFile {
    * \param fapl_id: `hid_t` File access property list identifier (default is `H5P_DEFAULT`).
    */
   HDFFile(
-	  std::string name, unsigned int flags=H5F_ACC_EXCL,
+	  const std::string& name, unsigned int flags=H5F_ACC_EXCL,
 	  hid_t fcpl_id=H5P_DEFAULT, hid_t fapl_id=H5P_DEFAULT
   );
 
@@ -105,7 +106,7 @@ class HDFFile {
    * \return hdf_file: `HDFFile *` Pointer to a new, open HDF5 file.
    */
   static HDFFile* from_schema(
-			      FileSchema &schema, std::string name, unsigned int flags=H5F_ACC_EXCL,
+			      FileSchema &schema, const std::string& name, unsigned int flags=H5F_ACC_EXCL,
 			      hid_t fcpl_id=H5P_DEFAULT, hid_t fapl_id=H5P_DEFAULT
   );
 
@@ -133,7 +134,7 @@ class HDFFile {
    * \return status: `herr_t` Exit status of the operation.
    */
   herr_t read(
-	       std::string dataset_name, std::string data_type, void *buffer,
+	       const std::string& dataset_name, const std::string& data_type, void *buffer,
 	       hid_t mem_space_id=H5S_ALL, hid_t file_space_id=H5S_ALL,
 	       hid_t dapl_id=H5P_DEFAULT, hid_t dxpl_id=H5P_DEFAULT
   );
@@ -150,7 +151,7 @@ class HDFFile {
    * \return status: `herr_t` Exit status of the operation.
    */
   herr_t write(
-	       std::string dataset_name, std::string data_type, const void *buffer,
+	       const std::string& dataset_name, const std::string& data_type, const void *buffer,
 	       hid_t mem_space_id=H5S_ALL, hid_t file_space_id=H5S_ALL,
 	       hid_t dapl_id=H5P_DEFAULT, hid_t dxpl_id=H5P_DEFAULT
   );

src/include/logging.h

@@ -67,7 +67,7 @@ class Logger {
    *
    * \param message: `string` The message to be printed.
    */
-  void err(std::string message);
+  void err(const std::string& message);
 
   /*! \brief Print a summary of recurrent messages and clear the stack.
    *
@@ -80,7 +80,7 @@ class Logger {
    * \param message: `string` The message to be printed.
    * \param level: `int` The priority level (default is `LOG_INFO = 1`).
    */
-  void log(std::string message, int level=LOG_INFO);
+  void log(const std::string& message, int level=LOG_INFO);
   
   /*! \brief Push a recurrent message to the message stack.
    *
@@ -92,7 +92,7 @@ class Logger {
    *
    * \param message: `string` The message to be stacked.
    */
-  void push(std::string message);
+  void push(const std::string& message);
 };
 
 #endif

src/include/tfrfme.h

@@ -27,26 +27,26 @@ protected:
    * \param file_name: `string` Name of the file to be loaded.
    * \return instance: `Swap1 *` Pointer to a new Swap1 instance.
    */
-  static Swap1 *from_hdf5(std::string file_name);
+  static Swap1 *from_hdf5(const std::string& file_name);
 
   /*! \brief Load a Swap1 instance from a legacy binary file.
    *
    * \param file_name: `string` Name of the file to be loaded.
    * \return instance: `Swap1 *` Pointer to a new Swap1 instance.
    */
-  static Swap1 *from_legacy(std::string file_name);
+  static Swap1 *from_legacy(const std::string& file_name);
 
   /*! \brief Save a Swap1 instance to a HDF5 binary file.
    *
    * \param file_name: `string` Name of the file to be written.
    */
-  void write_hdf5(std::string file_name);
+  void write_hdf5(const std::string& file_name);
 
   /*! \brief Save a Swap1 instance to a legacy binary file.
    *
    * \param file_name: `string` Name of the file to be written.
    */
-  void write_legacy(std::string file_name);
+  void write_legacy(const std::string& file_name);
 
 public:
   /*! \brief Swap1 instance constructor.
@@ -73,7 +73,7 @@ public:
    * or "LGEACY". Default is "LEGACY").
    * \return instance: `Swap1 *` Pointer to a newly created Swap1 instance.
    */
-  static Swap1* from_binary(std::string file_name, std::string mode="LEGACY");
+  static Swap1* from_binary(const std::string& file_name, const std::string& mode="LEGACY");
 
   /*! \brief Calculate the necessary amount of memory to create a new instance.
    *
@@ -99,7 +99,7 @@ public:
    * \param mode: `string` Format of the file (can be either "HDF5"
    * or "LGEACY". Default is "LEGACY").
    */
-  void write_binary(std::string file_name, std::string mode="LEGACY");
+  void write_binary(const std::string& file_name, const std::string& mode="LEGACY");
 
   /*! \brief Test whether two instances of Swap1 are equal.
    *
@@ -157,26 +157,26 @@ protected:
    * \param file_name: `string` Name of the file to be loaded.
    * \return instance: `Swap2 *` Pointer to a new Swap2 instance.
    */
-  static Swap2 *from_hdf5(std::string file_name);
+  static Swap2 *from_hdf5(const std::string& file_name);
 
   /*! \brief Load a Swap2 instance from a legacy binary file.
    *
    * \param file_name: `string` Name of the file to be loaded.
    * \return instance: `Swap2 *` Pointer to a new Swap2 instance.
    */
-  static Swap2 *from_legacy(std::string file_name);
+  static Swap2 *from_legacy(const std::string& file_name);
 
   /*! \brief Save a Swap2 instance to a HDF5 binary file.
    *
    * \param file_name: `string` Name of the file to be written.
    */
-  void write_hdf5(std::string file_name);
+  void write_hdf5(const std::string& file_name);
 
   /*! \brief Save a Swap2 instance to a legacy binary file.
    *
    * \param file_name: `string` Name of the file to be written.
    */
-  void write_legacy(std::string file_name);
+  void write_legacy(const std::string& file_name);
 
 public:
   /*! \brief Swap2 instance constructor.
@@ -196,7 +196,7 @@ public:
    * or "LGEACY". Default is "LEGACY").
    * \return instance: `Swap2 *` Pointer to a newly created Swap2 instance.
    */
-  static Swap2* from_binary(std::string file_name, std::string mode="LEGACY");
+  static Swap2* from_binary(const std::string& file_name, const std::string& mode="LEGACY");
 
   /*! \brief Get the pointer to the VKZM matrix.
    *
@@ -216,7 +216,7 @@ public:
    * \param param_name: `string` Name of the parameter.
    * \return value: `double` The value of the requested parameter.
    */
-  double get_param(std::string param_name);
+  double get_param(const std::string& param_name);
 
   /*! \brief Get the pointer to the VKV vector.
    *
@@ -249,7 +249,7 @@ public:
    * \param param_name: `string` Name of the parameter.
    * \param value: `double` The value of the parameter.
    */
-  void set_param(std::string param_name, double value);
+  void set_param(const std::string& param_name, double value);
 
   /*! \brief Write a Swap2 instance to binary file.
    *
@@ -257,7 +257,7 @@ public:
    * \param mode: `string` Format of the file (can be either "HDF5"
    * or "LGEACY". Default is "LEGACY").
    */
-  void write_binary(std::string file_name, std::string mode="LEGACY");
+  void write_binary(const std::string& file_name, const std::string& mode="LEGACY");
 
   /*! \brief Test whether two instances of Swap2 are equal.
    *
@@ -321,7 +321,7 @@ protected:
    * \return instance: `TFRFME *` Pointer to a new trapping configuration
    * instance.
    */
-  static TFRFME *from_hdf5(std::string file_name);
+  static TFRFME *from_hdf5(const std::string& file_name);
 
   /*! \brief Load a configuration instance from a legacy binary file.
    *
@@ -329,19 +329,19 @@ protected:
    * \return instance: `TFRFME *` Pointer to a new trapping configuration
    * instance.
    */
-  static TFRFME *from_legacy(std::string file_name);
+  static TFRFME *from_legacy(const std::string& file_name);
 
   /*! \brief Save a configuration instance to a HDF5 binary file.
    *
    * \param file_name: `string` Name of the file to be written.
    */
-  void write_hdf5(std::string file_name);
+  void write_hdf5(const std::string& file_name);
 
   /*! \brief Save a configuration instance to a legacy binary file.
    *
    * \param file_name: `string` Name of the file to be written.
    */
-  void write_legacy(std::string file_name);
+  void write_legacy(const std::string& file_name);
 
 public:
   /*! \brief Trapping configuration instance constructor.
@@ -369,7 +369,7 @@ public:
    * \return instance: `TFRFME *` Pointer to a newly created configuration
    * instance.
    */
-  static TFRFME* from_binary(std::string file_name, std::string mode="LEGACY");
+  static TFRFME* from_binary(const std::string& file_name, const std::string& mode="LEGACY");
 
   /*! \brief Get the pointer to the WSUM matrix.
    *
@@ -397,7 +397,7 @@ public:
    * \param param_name: `string` Name of the parameter.
    * \return value: `double` The value of the requested parameter.
    */
-  double get_param(std::string param_name);
+  double get_param(const std::string& param_name);
 
   /*! \brief Get the pointer to the X coordinates vector.
    *
@@ -422,7 +422,7 @@ public:
    * \param param_name: `string` Name of the parameter.
    * \param value: `double` Value to be stored as parameter.
    */
-  void set_param(std::string param_name, double value);
+  void set_param(const std::string& param_name, double value);
 
   /*! \brief Write a trapping configuration instance to binary file.
    *
@@ -430,7 +430,7 @@ public:
    * \param mode: `string` Format of the file (can be either "HDF5"
    * or "LGEACY". Default is "LEGACY").
    */
-  void write_binary(std::string file_name, std::string mode="LEGACY");
+  void write_binary(const std::string& file_name, const std::string& mode="LEGACY");
   
   /*! \brief Test whether two instances of configuration are equal.
    *
@@ -439,6 +439,6 @@ public:
    * \return result: `bool` True, if the two instances are equal,
    * false otherwise.
    */
-  bool operator ==(TFRFME &other);
+  bool operator ==(const TFRFME& other);
 };
 #endif

src/libnptm/Commons.cpp

@@ -2,28 +2,29 @@
 
 /*! \file Commons.cpp
  *
- * DEVELOPMENT NOTE:
- * The construction of common blocks requires some information
- * that is stored in configuration objects and is needed to compute
- * the allocation size of vectors and matrices. Currently, this
- * information is passed as arguments to the constructors of the
- * common blocks. A simpler and more logical way to operate would
- * be to design the constructors to take as arguments only pointers
- * to the configuration objects. These, on their turn, need to
- * expose methods to access the relevant data in read-only mode.
+ * \brief Implementation of the common data structures.
  */
 #ifndef INCLUDE_TYPES_H_
 #include "../include/types.h"
 #endif
 
+#ifndef INCLUDE_CONFIGURATION_H_
+#include "../include/Configuration.h"
+#endif
+
 #ifndef INCLUDE_COMMONS_H
 #include "../include/Commons.h"
 #endif
 
-C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
-  nlmmt = 2 * (l_max * (l_max + 2));
-  nsph = ns;
-  lm = l_max;
+C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
+  lm = gconf->l_max;
+  int li = gconf->li;
+  int le = gconf->le;
+  if (lm == 0) {
+    lm = (li > le) ? li : le;
+  }
+  nsph = gconf->number_of_spheres;
+  nlmmt = 2 * (lm * (lm + 2));
 
   vec_rmi = new dcomplex[nsph * lm]();
   vec_rei = new dcomplex[nsph * lm]();
@@ -36,23 +37,20 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
   vec_w = new dcomplex[4 * nlmmt]();
   w = new dcomplex*[nlmmt];
   for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);
-  configurations = 0;
-  for (int ci = 1; ci <= nsph; ci++) {
-    if (_iog[ci - 1] >= ci) configurations++;
-  }
+  configurations = sconf->configurations;
   vint = new dcomplex[16]();
   vec_vints = new dcomplex[nsph * 16]();
   vints = new dcomplex*[nsph];
-  rc = new double*[nsph];
-  nshl = new int[nsph]();
+  rc = new double*[configurations];
+  nshl = new int[configurations]();
   iog = new int[nsph]();
   int conf_index = 0;
   for (int vi = 0; vi < nsph; vi++) {
     vints[vi] = &(vec_vints[vi * 16]);
-    iog[vi] = _iog[vi];
+    iog[vi] = sconf->get_iog(vi);
     if (iog[vi] >= vi + 1) {
-      nshl[conf_index] = _nshl[conf_index];
-      rc[conf_index] = new double[_nshl[conf_index]]();
+      nshl[conf_index] = sconf->get_nshl(conf_index);
+      rc[conf_index] = new double[nshl[conf_index]]();
       conf_index++;
     }
   }
@@ -67,7 +65,17 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
   rxx = new double[nsph]();
   ryy = new double[nsph]();
   rzz = new double[nsph]();
-  ros = new double[nsph]();
+  if (nsph > 1) {
+    for (int c1i = 0; c1i < nsph; c1i++) {
+      rxx[c1i] = gconf->get_sph_x(c1i);
+      ryy[c1i] = gconf->get_sph_y(c1i);
+      rzz[c1i] = gconf->get_sph_z(c1i);
+    }
+  }
+  ros = new double[configurations]();
+  for (int ri = 0; ri < configurations; ri++) {
+    ros[ri] = sconf->get_radius(ri);
+  }
 
   sas = new dcomplex**[nsph];
   for (int si = 0; si < nsph; si++) {
@@ -209,9 +217,10 @@ C1_AddOns::C1_AddOns(C4 *c4) {
   vj = new dcomplex[1]();
   vyhj = new dcomplex[(nsph * nsph - 1) * litpos]();
   vyj0 = new dcomplex[nsph * lmtpos]();
+  am0v = new dcomplex[nlemt * nlemt]();
   am0m = new dcomplex*[nlemt];
   for (int ai = 0; ai < nlemt; ai++) {
-    am0m[ai] = new dcomplex[nlemt]();
+    am0m[ai] = (am0v + nlemt * ai);
   }
   vintm = new dcomplex[16]();
   vintt = new dcomplex[16]();
@@ -263,10 +272,11 @@ C1_AddOns::C1_AddOns(const C1_AddOns& rhs) {
   int vyj0size = nsph * lmtpos;
   vyj0 = new dcomplex[vyj0size]();
   for (int hi=0; hi<vyj0size; hi++) vyj0[hi] = rhs.vyj0[hi];
+  am0v = new dcomplex[nlemt * nlemt]();
   am0m = new dcomplex*[nlemt];
   for (int ai = 0; ai < nlemt; ai++) {
-    am0m[ai] = new dcomplex[nlemt]();
-    for (int aj = 0; aj < nlemt; aj++) am0m[ai][aj] = rhs.am0m[ai][aj];
+    for (int aj = 0; aj < nlemt; aj++) am0v[nlemt * ai + aj] = rhs.am0v[nlemt * ai + aj];
+    am0m[ai] = (am0v + nlemt * ai);
   }
   vintm = new dcomplex[16]();
   vintt = new dcomplex[16]();
@@ -326,10 +336,8 @@ C1_AddOns::~C1_AddOns() {
   delete[] vh;
   delete[] vj0;
   delete[] vj;
-  for (int ai = nlemt - 1; ai > -1; ai--) {
-    delete[] am0m[ai];
-  }
   delete[] am0m;
+  delete[] am0v;
   delete[] vintm;
   delete[] vintt;
   for (int fi = 1; fi > -1; fi--) {
@@ -350,11 +358,14 @@ C1_AddOns::~C1_AddOns() {
   delete[] ecsc;
 }
 
-C2::C2(int ns, int _nl, int npnt, int npntts) {
-  nsph = ns;
+C2::C2(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
+  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 = _nl;
+  nl = sconf->configurations;
+  if (nsph == 1 && nl == 1) nl = 5;
   ris = new dcomplex[nhspo]();
   dlri = new dcomplex[nhspo]();
   vkt = new dcomplex[nsph]();
@@ -420,20 +431,20 @@ C3::~C3() {
   delete[] tsas;
 }
 
-C4::C4(int _li, int _le, int _nsph) {
-    li = _li;
-    le = _le;
-    lm = (li > le) ? li : le;
-    nv3j = (lm * (lm  + 1) * (2 * lm + 7)) / 6;
-    nsph = _nsph;
-    // The following is needed to initialize C1_AddOns
-    litpo = li + li + 1;
-    litpos = litpo * litpo;
-    lmtpo = li + le + 1;
-    lmtpos = lmtpo * lmtpo;
-    nlim = li * (li + 2);
-    nlem = le * (le + 2);
-    lmpo = lm + 1;
+C4::C4(GeometryConfiguration *gconf) {
+  li = gconf->li;
+  le = gconf->le;
+  lm = (li > le) ? li : le;
+  nv3j = (lm * (lm  + 1) * (2 * lm + 7)) / 6;
+  nsph = gconf->number_of_spheres;
+  // The following is needed to initialize C1_AddOns
+  litpo = li + li + 1;
+  litpos = litpo * litpo;
+  lmtpo = li + le + 1;
+  lmtpos = lmtpo * lmtpo;
+  nlim = li * (li + 2);
+  nlem = le * (le + 2);
+  lmpo = lm + 1;
 }
 
 C4::C4(const C4& rhs) {
@@ -517,3 +528,387 @@ C9::~C9() {
   for (int si = sam_size_0 - 1; si > -1; si--) delete[] sam[si];
   delete[] sam;
 }
+
+ClusterIterationData::ClusterIterationData(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
+  c1 = new C1(gconf, sconf);
+  c2 = new C2(gconf, sconf);
+  c3 = new C3();
+  c4 = new C4(gconf);
+  c1ao = new C1_AddOns(c4);
+  c6 = new C6(c4->lmtpo);
+  const int ndi = c4->nsph * c4->nlim;
+  const np_int ndit = 2 * ndi;
+  c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
+  gaps = new double[c4->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[c4->nsph]();
+    tqspe[ti] = new dcomplex[c4->nsph]();
+    tqss[ti] = new double[c4->nsph]();
+    tqsps[ti] = new dcomplex[c4->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]();
+  }
+  zpv = new double***[c4->lm];
+  for (int zi = 0; zi < c4->lm; zi++) {
+    zpv[zi] = new double**[3];
+    for (int zj = 0; zj < 3; zj++) {
+      zpv[zi][zj] = new double*[2];
+      for (int zk = 0; zk < 2; zk++) {
+	zpv[zi][zj][zk] = new double[2]();
+      }
+    }
+  }
+  am_vector = new dcomplex[ndit * ndit]();
+  am = new dcomplex*[ndit];
+  for (int ai = 0; ai < ndit; ai++) {
+    am[ai] = (am_vector + ai * ndit);
+  }
+  
+  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;
+  qsfi = 0.0;
+  // End of suspect initializations
+  wn = sconf->wp / 3.0e8;
+  xip = sconf->xip;
+  sqsfi = 1.0;
+  vk = 0.0;
+}
+
+ClusterIterationData::ClusterIterationData(const ClusterIterationData& rhs) {
+  c1 = new C1(*(rhs.c1));
+  c2 = new C2(*(rhs.c2));
+  c3 = new C3(*(rhs.c3));
+  c4 = new C4(*(rhs.c4));
+  c1ao = new C1_AddOns(*(rhs.c1ao));
+  c6 = new C6(*(rhs.c6));
+  const int ndi = c4->nsph * c4->nlim;
+  const np_int ndit = 2 * ndi;
+  c9 = new C9(*(rhs.c9));
+  gaps = new double[c4->nsph]();
+  for (int gi = 0; gi < c4->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[c4->nsph]();
+    tqspe[ti] = new dcomplex[c4->nsph]();
+    tqss[ti] = new double[c4->nsph]();
+    tqsps[ti] = new dcomplex[c4->nsph]();
+    for (int tj = 0; tj < c4->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];
+    }
+  }
+  zpv = new double***[c4->lm];
+  for (int zi = 0; zi < c4->lm; zi++) {
+    zpv[zi] = new double**[3];
+    for (int zj = 0; zj < 3; zj++) {
+      zpv[zi][zj] = new double*[2];
+      for (int zk = 0; zk < 2; zk++) {
+	zpv[zi][zj][zk] = new double[2]();
+	zpv[zi][zj][zk][0] = rhs.zpv[zi][zj][zk][0];
+	zpv[zi][zj][zk][1] = rhs.zpv[zi][zj][zk][1];
+      }
+    }
+  }
+  am_vector = new dcomplex[ndit * ndit]();
+  for (np_int ai = 0; ai < ndit * ndit; ai++) am_vector[ai] = rhs.am_vector[ai];
+  am = new dcomplex*[ndit];
+  for (np_int ai = 0; ai < ndit; ai++) {
+    am[ai] = (am_vector + ai * ndit);
+  }
+  
+  arg = rhs.arg;
+  // These are suspect initializations
+  scan = rhs.scan;
+  cfmp = rhs.cfmp;
+  sfmp = rhs.sfmp;
+  cfsp = rhs.cfsp;
+  sfsp = rhs.sfsp;
+  qsfi = rhs.qsfi;
+  // End of suspect initializations
+  wn = rhs.wn;
+  xip = rhs.xip;
+  sqsfi = rhs.sqsfi;
+  vk = rhs.vk;
+}
+
+ClusterIterationData::~ClusterIterationData() {
+  const int nsph = c4->nsph;
+  delete[] am_vector;
+  delete[] am;
+  for (int zi = c4->lm - 1; zi > -1; zi--) {
+    for (int zj = 2; zj > -1; zj--) {
+      delete[] zpv[zi][zj][1];
+      delete[] zpv[zi][zj][0];
+      delete[] zpv[zi][zj];
+    }
+    delete[] zpv[zi];
+  }
+  delete[] zpv;
+  delete c1;
+  delete c1ao;
+  delete c2;
+  delete c3;
+  delete c4;
+  delete c6;
+  delete c9;
+  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;
+}
+
+ScatteringAngles::ScatteringAngles(GeometryConfiguration *gconf) {
+  int isam = gconf->isam;
+  _th = gconf->in_theta_start;
+  _thstp = gconf->in_theta_step;
+  _thlst = gconf->in_theta_end;
+  _ths = gconf->sc_theta_start;
+  _thsstp = gconf->sc_theta_step;
+  _thslst = gconf->sc_theta_end;
+  _ph = gconf->in_phi_start;
+  _phstp = gconf->in_phi_step;
+  _phlst = gconf->in_phi_end;
+  _phs = gconf->sc_phi_start;
+  _phsstp = gconf->sc_phi_step;
+  _phslst = gconf->sc_phi_end;
+  double small = 1.0e-3;
+  _nth = 0;
+  _nph = 0;
+  if (_thstp != 0.0) _nth = int((_thlst - _th) / _thstp + small);
+  _nth++;
+  if (_phstp != 0.0) _nph = int((_phlst - _ph) / _phstp + small);
+  _nph++;
+  _nths = 0;
+  _nphs = 0;
+  _thsca = 0.0;
+  if (isam > 1) {
+    _nths = 1;
+    _thsca = _ths - _th;
+  } else { // ISAM <= 1
+    if (_thsstp == 0.0) _nths = 0;
+    else _nths = int ((_thslst - _ths) / _thsstp + small);
+    _nths++;
+  }
+  if (isam >= 1) {
+    _nphs = 1;
+  } else {
+    if (_phsstp == 0.0) _nphs = 0;
+    else _nphs = int((_phslst - _phs) / _phsstp + small);
+    _nphs++;
+  }
+  _nk = nth * nph;
+  _nks = nths * nphs;
+  _nkks = nk * nks;
+}
+
+ScatteringAngles::ScatteringAngles(const ScatteringAngles &rhs) {
+  _th = rhs._th;
+  _thstp = rhs._thstp;
+  _thlst = rhs._thlst;
+  _ths = rhs._ths;
+  _thsstp = rhs._thsstp;
+  _thslst = rhs._thslst;
+  _ph = rhs._ph;
+  _phstp = rhs._phstp;
+  _phlst = rhs._phlst;
+  _phs = rhs._phs;
+  _phsstp = rhs._phsstp;
+  _phslst = rhs._phslst;
+  _nth = rhs._nth;
+  _nph = rhs._nph;
+  _nths = rhs._nths;
+  _nphs = rhs._nphs;
+  _thsca = rhs._thsca;
+  _nk = rhs._nk;
+  _nks = rhs._nks;
+  _nkks = rhs._nkks;
+}

src/libnptm/Parsers.cpp

@@ -19,7 +19,7 @@
 #include "../include/Parsers.h"
 #endif
 
-std::string *load_file(std::string file_name, int *count = 0) {
+std::string *load_file(const std::string& file_name, int *count = 0) {
   std::fstream input_file(file_name.c_str(), std::ios::in);
   List<std::string> *file_lines = new List<std::string>();
   std::string line;

src/libnptm/TransitionMatrix.cpp

@@ -94,7 +94,7 @@ TransitionMatrix::TransitionMatrix(
   }
 }
 
-TransitionMatrix* TransitionMatrix::from_binary(string file_name, string mode) {
+TransitionMatrix* TransitionMatrix::from_binary(const std::string& file_name, const std::string& mode) {
   TransitionMatrix *tm = NULL;
   if (mode.compare("LEGACY") == 0) {
     tm = TransitionMatrix::from_legacy(file_name);
@@ -107,7 +107,7 @@ TransitionMatrix* TransitionMatrix::from_binary(string file_name, string mode) {
   return tm;
 }
 
-TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
+TransitionMatrix* TransitionMatrix::from_hdf5(const std::string& file_name) {
   TransitionMatrix *tm = NULL;
   unsigned int flags = H5F_ACC_RDONLY;
   HDFFile *hdf_file = new HDFFile(file_name, flags);
@@ -160,7 +160,7 @@ TransitionMatrix* TransitionMatrix::from_hdf5(string file_name) {
   return tm;
 }
 
-TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
+TransitionMatrix* TransitionMatrix::from_legacy(const std::string& file_name) {
   fstream ttms;
   TransitionMatrix *tm = NULL;
   ttms.open(file_name, ios::binary | ios::in);
@@ -207,7 +207,7 @@ TransitionMatrix* TransitionMatrix::from_legacy(string file_name) {
   return tm;
 }
 
-void TransitionMatrix::write_binary(string file_name, string mode) {
+void TransitionMatrix::write_binary(const std::string& file_name, const std::string& mode) {
   if (mode.compare("LEGACY") == 0) {
     write_legacy(file_name);
   } else if (mode.compare("HDF5") == 0) {
@@ -218,7 +218,36 @@ void TransitionMatrix::write_binary(string file_name, string mode) {
   }
 }
 
-void TransitionMatrix::write_hdf5(string file_name) {
+void TransitionMatrix::write_binary(
+				   const std::string& file_name, np_int _nlemt, int _lm, double _vk,
+				   double _exri, dcomplex **_am0m, const std::string& mode
+) {
+  if (mode.compare("LEGACY") == 0) {
+    write_legacy(file_name, _nlemt, _lm, _vk, _exri, _am0m);
+  } else if (mode.compare("HDF5") == 0) {
+    write_hdf5(file_name, _nlemt, _lm, _vk, _exri, _am0m);
+  } else {
+    string message = "Unknown format mode: \"" + mode + "\"";
+    throw UnrecognizedFormatException(message);
+  }
+}
+
+void TransitionMatrix::write_binary(
+				    const std::string& file_name, int _lm, double _vk, double _exri,
+				    dcomplex **_rmi, dcomplex **_rei, double _sphere_radius,
+				    const std::string& mode
+) {
+  if (mode.compare("LEGACY") == 0) {
+    write_legacy(file_name, _lm, _vk, _exri, _rmi, _rei, _sphere_radius);
+  } else if (mode.compare("HDF5") == 0) {
+    write_hdf5(file_name, _lm, _vk, _exri, _rmi, _rei, _sphere_radius);
+  } else {
+    string message = "Unknown format mode: \"" + mode + "\"";
+    throw UnrecognizedFormatException(message);
+  }
+}
+
+void TransitionMatrix::write_hdf5(const std::string& file_name) {
   if (is == 1 || is == 1111) {
     List<string> rec_name_list(1);
     List<string> rec_type_list(1);
@@ -237,17 +266,10 @@ void TransitionMatrix::write_hdf5(string file_name) {
     rec_type_list.append("FLOAT64_(1)");
     rec_ptr_list.append(&exri);
     rec_name_list.append("ELEMENTS");
-    str_type = "FLOAT64_(" + to_string(shape[0]) + "," + to_string(2 * shape[1]) + ")";
+    str_type = "COMPLEX128_(" + to_string(shape[0]) + "," + to_string(shape[1]) + ")";
     rec_type_list.append(str_type);
-    // The (N x M) matrix of complex is converted to a (N x 2M) matrix of double,
-    // where REAL(E_i,j) -> E_i,(2 j) and IMAG(E_i,j) -> E_i,(2 j + 1)
-    int num_elements = 2 * shape[0] * shape[1];
-    double *ptr_elements = new double[num_elements]();
-    for (int ei = 0; ei < num_elements / 2; ei++) {
-      ptr_elements[2 * ei] = real(elements[ei]);
-      ptr_elements[2 * ei + 1] = imag(elements[ei]);
-    }
-    rec_ptr_list.append(ptr_elements);
+    dcomplex *p_first = elements;
+    rec_ptr_list.append(p_first);
     if (is == 1111) {
       rec_name_list.append("RADIUS");
       rec_type_list.append("FLOAT64_(1)");
@@ -264,7 +286,7 @@ void TransitionMatrix::write_hdf5(string file_name) {
       hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
     hdf_file->close();
     
-    delete[] ptr_elements;
+    p_first = NULL;
     delete[] rec_names;
     delete[] rec_types;
     delete[] rec_pointers;
@@ -275,7 +297,104 @@ void TransitionMatrix::write_hdf5(string file_name) {
   }
 }
 
-void TransitionMatrix::write_legacy(string file_name) {
+void TransitionMatrix::write_hdf5(
+				  const std::string& file_name, np_int _nlemt, int _lm, double _vk,
+				  double _exri, dcomplex **_am0m
+) {
+  int is = 1;
+  List<string> rec_name_list(1);
+  List<string> rec_type_list(1);
+  List<void *> rec_ptr_list(1);
+  string str_type, str_name;
+  rec_name_list.set(0, "IS");
+  rec_type_list.set(0, "INT32_(1)");
+  rec_ptr_list.set(0, &is);
+  rec_name_list.append("L_MAX");
+  rec_type_list.append("INT32_(1)");
+  rec_ptr_list.append(&_lm);
+  rec_name_list.append("VK");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&_vk);
+  rec_name_list.append("EXRI");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&_exri);
+  rec_name_list.append("ELEMENTS");
+  str_type = "COMPLEX128_(" + to_string(_nlemt) + "," + to_string(_nlemt) + ")";
+  rec_type_list.append(str_type);
+  // The (N x M) matrix of complex is converted to a (N x 2M) matrix of double,
+  // where REAL(E_i,j) -> E_i,(2 j) and IMAG(E_i,j) -> E_i,(2 j + 1)
+  dcomplex *p_first = _am0m[0];
+  rec_ptr_list.append(p_first);
+  
+  string *rec_names = rec_name_list.to_array();
+  string *rec_types = rec_type_list.to_array();
+  void **rec_pointers = rec_ptr_list.to_array();
+  const int rec_num = rec_name_list.length();
+  FileSchema schema(rec_num, rec_types, rec_names);
+  HDFFile *hdf_file = HDFFile::from_schema(schema, file_name, H5F_ACC_TRUNC);
+  for (int ri = 0; ri < rec_num; ri++)
+    hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
+  hdf_file->close();
+
+  p_first = NULL;
+  delete[] rec_names;
+  delete[] rec_types;
+  delete[] rec_pointers;
+  delete hdf_file;
+}
+
+void TransitionMatrix::write_hdf5(
+				  const std::string& file_name, int _lm, double _vk, double _exri,
+				  dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
+) {
+  int is = 1111;
+  List<string> rec_name_list(1);
+  List<string> rec_type_list(1);
+  List<void *> rec_ptr_list(1);
+  string str_type, str_name;
+  rec_name_list.set(0, "IS");
+  rec_type_list.set(0, "INT32_(1)");
+  rec_ptr_list.set(0, &is);
+  rec_name_list.append("L_MAX");
+  rec_type_list.append("INT32_(1)");
+  rec_ptr_list.append(&_lm);
+  rec_name_list.append("VK");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&_vk);
+  rec_name_list.append("EXRI");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&_exri);
+  dcomplex *_elements = new dcomplex[2 * _lm]();
+  for (int ei = 0; ei < _lm; ei++) {
+    _elements[2 * ei] = -1.0 / _rmi[ei][0];
+    _elements[2 * ei + 1] = -1.0 / _rei[ei][0];
+  }
+  rec_name_list.append("ELEMENTS");
+  str_type = "COMPLEX128_(" + to_string(_lm) + "," + to_string(2) + ")";
+  rec_type_list.append(str_type);
+  rec_ptr_list.append(_elements);
+  rec_name_list.append("RADIUS");
+  rec_type_list.append("FLOAT64_(1)");
+  rec_ptr_list.append(&_sphere_radius);
+
+  string *rec_names = rec_name_list.to_array();
+  string *rec_types = rec_type_list.to_array();
+  void **rec_pointers = rec_ptr_list.to_array();
+  const int rec_num = rec_name_list.length();
+  FileSchema schema(rec_num, rec_types, rec_names);
+  HDFFile *hdf_file = HDFFile::from_schema(schema, file_name, H5F_ACC_TRUNC);
+  for (int ri = 0; ri < rec_num; ri++)
+    hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
+  hdf_file->close();
+    
+  delete[] _elements;
+  delete[] rec_names;
+  delete[] rec_types;
+  delete[] rec_pointers;
+  delete hdf_file;
+}
+
+void TransitionMatrix::write_legacy(const std::string& file_name) {
   fstream ttms;
   if (is == 1111 || is == 1) {
     ttms.open(file_name, ios::binary | ios::out);
@@ -308,6 +427,66 @@ void TransitionMatrix::write_legacy(string file_name) {
   }
 }
 
+void TransitionMatrix::write_legacy(
+				    const std::string& file_name, np_int _nlemt, int _lm, double _vk,
+				    double _exri, dcomplex **_am0m
+) {
+  fstream ttms;
+  int is = 1;
+  ttms.open(file_name, ios::binary | ios::out);
+  if (ttms.is_open()) {
+    ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
+    ttms.write(reinterpret_cast<char *>(&_lm), sizeof(int));
+    ttms.write(reinterpret_cast<char *>(&_vk), sizeof(double));
+    ttms.write(reinterpret_cast<char *>(&_exri), sizeof(double));
+    double rval, ival;
+    for (np_int ei = 0; ei < _nlemt; ei++) {
+      for (np_int ej = 0; ej < _nlemt; ej++) {
+	rval = real(_am0m[ei][ej]);
+	ival = imag(_am0m[ei][ej]);
+	ttms.write(reinterpret_cast<char *>(&rval), sizeof(double));
+	ttms.write(reinterpret_cast<char *>(&ival), sizeof(double));
+      }
+    }
+    ttms.close();
+  } else { // Failed to open output file. Should never happen.
+    printf("ERROR: could not open Transition Matrix file for writing.\n");
+  }
+}
+
+void TransitionMatrix::write_legacy(
+				    const std::string& file_name, int _lm, double _vk, double _exri,
+				    dcomplex **_rmi, dcomplex **_rei, double _sphere_radius
+) {
+  fstream ttms;
+  int is = 1111;
+  ttms.open(file_name, ios::binary | ios::out);
+  if (ttms.is_open()) {
+    ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
+    ttms.write(reinterpret_cast<char *>(&_lm), sizeof(int));
+    ttms.write(reinterpret_cast<char *>(&_vk), sizeof(double));
+    ttms.write(reinterpret_cast<char *>(&_exri), sizeof(double));
+    double rval, ival;
+    dcomplex element;
+    for (int ei = 0; ei < _lm; ei++) {
+      element = -1.0 / _rmi[ei][0];
+      rval = real(element);
+      ival = imag(element);
+      ttms.write(reinterpret_cast<char *>(&rval), sizeof(double));
+      ttms.write(reinterpret_cast<char *>(&ival), sizeof(double));
+      element = -1.0 / _rei[ei][0];
+      rval = real(element);
+      ival = imag(element);
+      ttms.write(reinterpret_cast<char *>(&rval), sizeof(double));
+      ttms.write(reinterpret_cast<char *>(&ival), sizeof(double));
+    }
+    ttms.write(reinterpret_cast<char *>(&_sphere_radius), sizeof(double));
+    ttms.close();
+  } else { // Failed to open output file. Should never happen.
+    printf("ERROR: could not open Transition Matrix file for writing.\n");
+  }
+}
+
 bool TransitionMatrix::operator ==(TransitionMatrix &other) {
   if (is != other.is) {
     return false;

src/libnptm/clu_subs.cpp

@@ -9,6 +9,10 @@
 #include "../include/types.h"
 #endif
 
+#ifndef INCLUDE_CONFIGURATION_H_
+#include "../include/Configuration.h"
+#endif
+
 #ifndef INCLUDE_COMMONS_H_
 #include "../include/Commons.h"
 #endif
@@ -1814,7 +1818,7 @@ void scr2(
   } // ipo1 loop
 }
 
-void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
+void str(ScattererConfiguration *sconf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
   dcomplex *ylm;
   const double pi = acos(-1.0);
   c3->gcs = 0.0;
@@ -1826,7 +1830,7 @@ void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
       c1->gcsv[i18 - 1] = gcss;
       int nsh = c1->nshl[i18 - 1];
       for (int j16 = 1; j16 <= nsh; j16++) {
-	c1->rc[i18 - 1][j16 - 1] = rcf[i18 - 1][j16 - 1] * c1->ros[i18 - 1];
+	c1->rc[i18 - 1][j16 - 1] = sconf->get_rcf(i18 - 1, j16 - 1) * c1->ros[i18 - 1];
       } // j16 loop
     }
     c3->gcs += gcss;

src/libnptm/sph_subs.cpp

@@ -8,6 +8,10 @@
 #include "../include/types.h"
 #endif
 
+#ifndef INCLUDE_CONFIGURATION_H_
+#include "../include/Configuration.h"
+#endif
+
 #ifndef INCLUDE_COMMONS_H_
 #include "../include/Commons.h"
 #endif