Merge branch 'sphere' into 'master'

Sphere

See merge request giacomo.mulas/np_tmcode!6

.gitignore

@@ -3,3 +3,4 @@ build/cluster/*
 build/sphere/*
 build/trapping/*
 doc/build/*
+

src/Makefile

 SUBDIRS := cluster sphere trapping
-BUILDDIR=../build
+SRCDIR=$(PWD)
+BUILDDIR=$(SRCDIR)/../build
+DOCSDIR=$(SRCDIR)/../doc
 
 all: $(SUBDIRS)
 
+docs:
+	cd $(DOCSDIR)/src; doxygen config.dox
+
 $(SUBDIRS):
 	$(MAKE) -C $@
 
@@ -15,5 +20,7 @@ wipe:
 	rm -f $(BUILDDIR)/cluster/*
 	rm -f $(BUILDDIR)/sphere/*
 	rm -f $(BUILDDIR)/trapping/*
+	if [ -d $(DOCSDIR)/build/html ]; then rm -r $(DOCSDIR)/build/html; fi
+	if [ -d $(DOCSDIR)/build/latex ]; then rm -r $(DOCSDIR)/build/latex; fi
 
 .PHONY: all $(SUBDIRS)

src/README.md

@@ -7,3 +7,27 @@ This directory collects the source code of the original programs and the develop
 The original code is contained in the folders named `cluster`, `sphere` and `trapping`. Each folder contains a `Makefile` to compile either the whole program set or the single programs. A global `Makefile`, which contains instructions to build all the original source code, is available directly in the `src` folder.
 
 In all cases, build commands executed through `make` will output the object files and the linked binaries in the proper folders under the build directory.
+
+## FORTRAN code setup and execution (requires `gfortran` and GNU `make`)
+
+1. cd to the `src` folder
+2. run `make`
+
+   > make
+
+3. cd to the `build/sphere` folder
+4. run `sph` following the instructions given in `build\README.md`
+
+## C++ code setup and execution (requires `g++` and GNU `make`)
+
+1. cd to the `src` folder
+2. run `make np_sphere`
+
+   > make np_sphere
+
+3. cd to the `build/sphere` folder
+4. run `np_sphere`
+
+   > ./np_sphere
+
+5. check the consistency between the text files named `OSPH` and `c_OSPH`

src/include/Commons.h

+/*! \file Commons.h
+ *
+ * \brief C++ porting of common data structures.
+ *
+ * Many functions of the original FORTRAN code share complex data blocks in
+ * form of COMMON blocks. This poses the limit of freezing the structure of
+ * the data blocks in the code, therefore implying the necessity to modify
+ * the code to adapt it to the input and to recompile before running. C++,
+ * on the contrary, offers the possibility to represent the necessary data
+ * structures as classes that can dinamically instantiate the shared information
+ * in the most convenient format for the current configuration. This approach
+ * adds an abstraction layer that lifts the need to modify and recompile the
+ * code depending on the input.
+ *
+ */
+
+#ifndef SRC_INCLUDE_COMMONS_
+#define SRC_INCLUDE_COMMONS_
+
+#include <complex>
+
+/*! \brief Representation of the FORTRAN C1 common blocks.
+ *
+ * C1 common blocks are used to store vector field expansions and geometric
+ * properties, such as sphere sizes, positions and cross-sections. These are
+ * used by functions such as `aps`, `diel`, `pwma`, `rabas`, `sscr0`, `sscr2`,
+ * `wmamp`, `wmasp` and `dme`. QUESTION: correct?
+ *
+ * Due to the necessity to share the class contents with many functions that
+ * need to read and update various fields, all shared members of C1 are declared
+ * public (i.e., the class is just a structure with more advanced constructor
+ * and destroyer). Further development may go in the direction of creating
+ * a better encapsulation, either by unpacking the contents (recommended) or
+ * by creating public methods to access protected fields (standard way, but not
+ * recommended for HPC applications).
+ */
+class C1 {
+protected:
+	//! \brief Number of spheres.
+	int nsph;
+	//! \brief Maximum order of field expansion.
+	int lm;
+	//! \brief NLMMT. QUESTION: definition?
+	int nlmmt;
+public:
+	//! \brief QUESTION: definition?
+	std::complex<double> **rmi;
+	//! \brief QUESTION: definition?
+	std::complex<double> **rei;
+	//! \brief QUESTION: definition?
+	std::complex<double> **w;
+	//! \brief QUESTION: definition?
+	std::complex<double> *fsas;
+	//! \brief QUESTION: definition?
+	std::complex<double> **vints;
+	//! \brief QUESTION: definition?
+	double *sscs;
+	//! \brief QUESTION: definition?
+	double *sexs;
+	//! \brief QUESTION: definition?
+	double *sabs;
+	//! \brief QUESTION: definition?
+	double *sqscs;
+	//! \brief QUESTION: definition?
+	double *sqexs;
+	//! \brief QUESTION: definition?
+	double *sqabs;
+	//! \brief QUESTION: definition?
+	double *gcsv;
+	//! \brief Vector of sphere X coordinates.
+	double *rxx;
+	//! \brief Vector of sphere X coordinates.
+	double *ryy;
+	//! \brief Vector of sphere X coordinates.
+	double *rzz;
+	//! \brief Vector of sphere radii.
+	double *ros;
+	//! \brief Matrix of spherical layer break radii. QUESTION: correct?
+	double **rc;
+	//! \brief Vector of spherical component identifiers.
+	int *iog;
+	//! \brief Vector of numbers of spherical layers.
+	int *nshl;
+	//! \brief QUESTION: definition?
+	std::complex<double> ***sas;
+
+	/*! \brief C1 instance constructor.
+	 *
+	 * \param ns: `int` Number of spheres.
+	 * \param l_max: `int` Maximum order of field expansion.
+	 */
+	C1(int ns, int l_max);
+
+	//! \brief C1 instance destroyer.
+	~C1();
+};
+
+/*! \brief Representation of the FORTRAN C2 blocks.
+ *
+ */
+class C2 {
+	//! \brief Number of spheres.
+	int nsph;
+	//! \brief Number of required orders.
+	int nhspo;
+public:
+	//! \brief QUESTION: definition?
+	std::complex<double> *ris;
+	//! \brief QUESTION: definition?
+	std::complex<double> *dlri;
+	//! \brief QUESTION: definition?
+	std::complex<double> *dc0;
+	//! \brief QUESTION: definition?
+	std::complex<double> *vkt;
+	//! Vector of scaled sizes. QUESTION: correct?
+	double *vsz;
+
+	/*! \brief C2 instance constructor.
+	 *
+	 * \param ns: `int` Number of spheres.
+	 * \param nl: `int`
+	 * \param npnt: `int`
+	 * \param npntts: `int`
+	 */
+	C2(int ns, int nl, int npnt, int npntts);
+
+	//! \brief C2 instance destroyer.
+	~C2();
+};
+
+#endif

src/include/Configuration.h

+/*! \file Configuration.h
+ */
+
+#ifndef INCLUDE_CONFIGURATION_H_
+#define INCLUDE_CONFIGURATION_H_
+
+#include <complex>
+#include <exception>
+#include <string>
+
+/**
+ * \brief Exception for open file error handlers.
+ */
+class OpenConfigurationFileException: public std::exception {
+protected:
+	//! \brief Name of the file that was accessed.
+	std::string file_name;
+
+public:
+	/**
+	 * \brief Exception instance constructor.
+	 *
+	 * \param name: `string` Name of the file that was accessed.
+	 */
+	OpenConfigurationFileException(std::string name) { file_name = name; }
+
+	/**
+	 * \brief Exception message.
+	 */
+	virtual const char* what() const throw() {
+		std::string message = "Error opening configuration file: " + file_name;
+		return message.c_str();
+	}
+};
+
+/**
+ * \brief Exception for unrecognized configuration data sets.
+ */
+class UnrecognizedConfigurationException: public std::exception {
+protected:
+	//! Description of the problem.
+	std::string message;
+public:
+	/**
+	 * \brief Exception instance constructor.
+	 *
+	 * \param problem: `string` Description of the problem that occurred.
+	 */
+	UnrecognizedConfigurationException(std::string problem) { message = problem; }
+	/**
+	 * \brief Exception message.
+	 */
+	virtual const char* what() const throw() {
+		return message.c_str();
+	}
+};
+
+/**
+ * \brief A class to represent the configuration of the scattering geometry.
+ *
+ * GeometryConfiguration is a class designed to store the necessary configuration
+ * data to describe the scattering geometry, including the distribution of the
+ * particle components, the orientation of the incident and scattered radiation
+ * fields and their polarization properties.
+ */
+class GeometryConfiguration {
+	//! Temporary work-around to allow sphere() peeking in.
+	friend void sphere();
+protected:
+	//! \brief Number of spherical components.
+	int number_of_spheres;
+	//! \brief Maximum expansion order of angular momentum.
+	int l_max;
+	//! \brief Incident field polarization status (0 - linear, 1 - circular).
+	int in_pol;
+	//! \brief Number of transition points. QUESTION: correct?
+	int npnt;
+	//! \brief Transition smoothness. QUESTION: correct?
+	int npntts;
+	//! \brief Type of meridional plane definition.
+	int meridional_type;
+	//! \brief Transition matrix layer ID. QUESTION: correct?
+	int jwtm;
+	//! \brief Incident field initial azimuth.
+	double in_theta_start;
+	//! \brief Incident field azimuth step.
+	double in_theta_step;
+	//! \brief Incident field final azimuth.
+	double in_theta_end;
+	//! \brief Scattered field initial azimuth.
+	double sc_theta_start;
+	//! \brief Scattered field azimuth step.
+	double sc_theta_step;
+	//! \brief Scattered field final azimuth.
+	double sc_theta_end;
+	//! \brief Incident field initial elevation.
+	double in_phi_start;
+	//! \brief Incident field elevation step.
+	double in_phi_step;
+	//! \brief Incident field final elevation.
+	double in_phi_end;
+	//! \brief Scattered field initial elevation.
+	double sc_phi_start;
+	//! \brief Scattered field elevation step.
+	double sc_phi_step;
+	//! \brief Scattered field final elevation.
+	double sc_phi_end;
+	//! \brief Vector of spherical components X coordinates.
+	double *sph_x;
+	//! \brief Vector of spherical components Y coordinates.
+	double *sph_y;
+	//! \brief Vector of spherical components Z coordinates.
+	double *sph_z;
+
+public:
+	/*! \brief Build a scattering geometry configuration structure.
+	 *
+	 * \param nsph: `int` Number of spheres to be used in calculation.
+	 * \param lm: `int` Maximum field angular momentum expansion order.
+	 * \param in_pol: `int` Incident field polarization status
+	 * \param npnt: `int` Number of transition points. QUESTION: correct?
+	 * \param npntts: `int` Transition smoothness. QUESTION: correct?
+	 * \param meridional_type: `int` Type of meridional plane definition (<0
+	 * for incident angles, 0 if determined by incidence and observation, =1
+	 * accross z-axis for incidence and observation, >1 across z-axis as a
+	 * function of incidence angles for fixed scattering).
+	 * \param x: `double*` Vector of spherical components X coordinates.
+	 * \param y: `double*` Vector of spherical components Y coordinates.
+	 * \param z: `double*` Vector of spherical components Z coordinates.
+	 * \param in_th_start: `double` Incident field starting azimuth angle.
+	 * \param in_th_step: `double` Incident field azimuth angle step.
+	 * \param in_th_end: `double` Incident field final azimuth angle.
+	 * \param sc_th_start: `double` Scattered field starting azimuth angle.
+	 * \param sc_th_step: `double` Scattered field azimuth angle step.
+	 * \param sc_th_end: `double` Scattered field final azimuth angle.
+	 * \param in_ph_start: `double` Incident field starting elevation angle.
+	 * \param in_ph_step: `double` Incident field elevation angle step.
+	 * \param in_ph_end: `double` Incident field final elevation angle.
+	 * \param sc_ph_start: `double` Scattered field starting elevation angle.
+	 * \param sc_ph_step: `double` Scattered field elevation angle step.
+	 * \param sc_ph_end: `double` Scattered field final elevation angle.
+	 * \param jwtm: `int` Transition Matrix layer ID. QUESTION: correct?
+	 */
+	GeometryConfiguration(
+			int nsph, int lm, int in_pol, int npnt, int npntts, int meridional_type,
+			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 in_ph_start, double in_ph_step, double in_ph_end,
+			double sc_ph_start, double sc_ph_step, double sc_ph_end,
+			int jwtm
+	);
+
+	/*! \brief Destroy a GeometryConfiguration instance.
+	 */
+	~GeometryConfiguration();
+
+	/*! \brief Build geometry configuration from legacy configuration input file.
+	 *
+	 * To allow for consistency tests and backward compatibility, geometry
+	 * configurations can be built from legacy configuration files. This function
+	 * replicates the approach implemented by the FORTRAN SPH and CLU codes, but
+	 * using a C++ oriented work-flow.
+	 *
+	 * \param file_name: `string` Name of the legacy configuration data file.
+	 * \return config: `GeometryConfiguration*` Pointer to object containing the
+	 * configuration data.
+	 */
+	static GeometryConfiguration *from_legacy(std::string file_name);
+};
+
+/**
+ * \brief A class to represent scatterer configuration objects.
+ *
+ * ScattererConfiguration is a class designed to store the necessary configuration
+ * data to describe the scatterer properties.
+ */
+class ScattererConfiguration {
+	//! Temporary work-around to allow sphere() peeking in.
+	friend void sphere();
+protected:
+	//! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x LAYERS].
+	std::complex<double> ***dc0_matrix;
+	//! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].
+	double *radii_of_spheres;
+	//! \brief Matrix of fractional transition radii with size [N_SPHERES x LAYERS].
+	double **rcf;
+	//! \brief Vector of sphere ID numbers, with size [N_SPHERES].
+	int *iog_vec;
+	//! \brief Vector of layer numbers for every sphere, with size [N_SPHERES].
+	int *nshl_vec;
+	//! \brief Vector of scale parameters, with size [N_SCALES].
+	double *scale_vec;
+	//! \brief Name of the reference variable type (one of XIV, WNS, WLS, PUS, EVS).
+	std::string reference_variable_name;
+	//! \brief Number of spherical components.
+	int number_of_spheres;
+	//! \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).
+	int idfc;
+	//! \brief External medium dielectric constant. QUESTION: correct?
+	double exdc;
+	//! \brief WP. QUESTION: better definition?
+	double wp;
+	//! \brief Peak XI. QUESTION: correct?
+	double xip;
+	//! \brief Flag to control whether to add an external layer.
+	bool use_external_sphere;
+
+public:
+	/*! \brief Build a scatterer configuration structure.
+	 *
+	 * 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: Matrix of reference dielectric constants.
+	 * \param has_external: `bool` Flag to set whether to add an external spherical layer.
+	 * \param exdc: `double` EXDC
+	 * \param wp: `double` wp
+	 * \param xip: `double` xip
+	 */
+	ScattererConfiguration(
+			int nsph,
+			double *scale_vector,
+			int nxi,
+			std::string variable_name,
+			int *iog_vector,
+			double *ros_vector,
+			int *nshl_vector,
+			double **rcf_vector,
+			int dielectric_func_type,
+			std::complex<double> ***dc_matrix,
+			bool has_external,
+			double exdc,
+			double wp,
+			double xip
+			);
+
+	/*! \brief Destroy a scatterer configuration instance.
+	 */
+	~ScattererConfiguration();
+
+	/*! \brief Build configuration from binary configuration input file.
+	 *
+	 * The configuration step can save configuration data as a binary file. The original
+	 * FORTRAN code used this possibility to manage communication between the configuring
+	 * code and the calculation program. This possibility is maintained, in case the
+	 * configuration step needs to be separated from the calculation execution. In this
+	 * case, `from_binary()` is the class method that restores a ScattererConfiguration
+	 * object from a previously saved binary file.
+	 *
+	 * \param file_name: `string` Name of the binary configuration data file.
+	 * \param mode: `string` Binary encoding. Can be one of "LEGACY", ... . Optional
+	 * (default is "LEGACY").
+	 * \return config: `ScattererConfiguration*` Pointer to object containing the
+	 * scatterer configuration data.
+	 */
+	static ScattererConfiguration* from_binary(std::string file_name, std::string mode = "LEGACY");
+
+	/*! \brief Build scatterer configuration from legacy configuration input file.
+	 *
+	 * To allow for consistency tests and backward compatibility, ScattererConfiguration
+	 * objects can be built from legacy configuration files. This function replicates
+	 * the approach implemented by the FORTRAN EDFB code, but using a C++ oriented
+	 * work-flow.
+	 *
+	 * \param file_name: `string` Name of the legacy configuration data file.
+	 * \return config: `ScattererConfiguration*` Pointer to object containing the
+	 * scatterer configuration data.
+	 */
+	static ScattererConfiguration* from_dedfb(std::string file_name);
+
+	/*! \brief Print the contents of the configuration object to terminal.
+	 *
+	 * In case of quick debug testing, `ScattererConfiguration.print()` allows printing
+	 * a formatted summary of the configuration data to terminal.
+	 */
+	void print();
+
+	/*! \brief Write the scatterer configuration data to binary output.
+	 *
+	 * The execution work-flow may be split in a configuration step and one or more
+	 * calculation steps. In case the calculation is not being run all-in-one, it can
+	 * be useful to save the configuration data. `ScattererConfiguration.write_binary()`
+	 * performs the operation of saving the configuration in binary format. This function
+	 * can work in legacy mode, to write backward compatible configuration files, as well
+	 * as by wrapping the data into common scientific formats (NB: this last option still
+	 * needs to be implemented).
+	 *
+	 * \param file_name: `string` Name of the file to be written.
+	 * \param mode: `string` Binary encoding. Can be one of "LEGACY", ... . Optional
+	 * (default is "LEGACY").
+	 */
+	void write_binary(std::string file_name, std::string mode = "LEGACY");
+
+	/*! \brief Write the scatterer configuration data to formatted text output.
+	 *
+	 * Writing configuration to formatted text is an optional operation, which may turn
+	 * out to be useful for consistency checks. As a matter of fact, formatted configuration
+	 * output is not read back by the FORTRAN code work-flow and it can be safely omitted,
+	 * unless there is a specific interest in assessing that the legacy code and the
+	 * updated one are doing the same thing.
+	 *
+	 * \param file_name: `string` Name of the file to be written.
+	 */
+	void write_formatted(std::string file_name);
+};
+
+#endif

src/include/List.h

 /*! \file List.h
  */
 
-#ifndef LIST_OUT_OF_BOUNDS_EXCEPTION
-#define LIST_OUT_OF_BOUNDS_EXCEPTION 1
-#endif
+#ifndef INCLUDE_LIST_H_
+#define INCLUDE_LIST_H_
+
+#include <exception>
+#include <string>
+
+/**
+ * \brief Exception for out of bounds List requests.
+ */
+class ListOutOfBoundsException: public std::exception {
+protected:
+	//! \brief Minimum index defined in the List.
+	int min_index;
+	//! \brief Maximum index defined in the List.
+	int max_index;
+	//! \brief List index requested by user.
+	int requested_index;
+
+public:
+	/**
+	 * \brief Exception instance constructor.
+	 *
+	 * \param requested: `int` The index that was requested.
+	 * \param min: `int` The minimum index allowed by the list.
+	 * \param max: `int` The maximum index allowed by the list.
+	 */
+	ListOutOfBoundsException(int requested, int min, int max) {
+		min_index = min;
+		max_index = max;
+		requested_index = requested;
+	}
+	/**
+	 * \brief Exception message.
+	 */
+	virtual const char* what() const throw() {
+		std::string message = "Error: requested index ";
+		message += requested_index;
+		message += " is out of range [";
+		message += min_index;
+		message += ", ";
+		message += (max_index - 1);
+		message += "]";
+		return message.c_str();
+	}
+};
 
 /**
  * \brief A class to represent dynamic lists.
@@ -16,7 +58,7 @@
  *
  * For this reason, the best use of List objects is to collect all the
  * desired members and then, once the element number is known, to convert
- * the List to C array, by calling List.to_array(). This function returns
+ * the List to C array, by calling `List.to_array()`. This function returns
  * a contiguous array of type T[SIZE] that can be used for indexed access.
  */
 template<class T> class List {
@@ -33,18 +75,19 @@ template<class T> class List {
  public:
   /*! \brief List constructor.
    *
-   * Use the constructor List<T>([int length]) to create a new list with a given
+   * Use the constructor `List<T>([int length])` to create a new list with a given
    * size. If the required size is not known in advance, it is recommended
    * to create a List with SIZE=1 (this is the default behavior) and then
-   * to append the elements dynamically, using List.append(ELEMENT) (where
+   * to append the elements dynamically, using `List.append(ELEMENT)` (where
    * ELEMENT needs to be a value of type T, corresponding to the class
    * template specialization). Note that, due to the default behavior, the
    * following calls are equivalent and they both produce an integer List
    * with size equal to 1:
    *
-   * a = List<int>(1);
-   *
-   * b = List<int>();
+   * \code{.cpp}
+   * List<int> a = List<int>(1);
+   * List<int> b = List<int>();
+   * \endcode
    *
    * \param length: `int` The size of the list to be constructed [OPTIONAL, default=1].
    */
@@ -52,7 +95,7 @@ template<class T> class List {
     size = length;
     first = new element;
     first->p_prev = NULL;
-    element *current = first;
+    current = first;
     element *p_prev = first;
     for (int i = 1; i < size; i++) {
       current = new element;
@@ -62,6 +105,8 @@ template<class T> class List {
     last = current;
   }
 
+  /*! \brief Destroy a List instance.
+   */
   ~List() {
     current = last;
     element *old;
@@ -72,14 +117,14 @@ template<class T> class List {
     }
   }
   
-  /*! \brief Append an element at the end of the list.
+  /*! \brief Append an element at the end of the List.
    *
    * To dynamically create a list whose size is not known in advance,
    * elements can be appended in an iterative way. Note that element
    * manipulation is much more effective in a C array than in a List
    * object. For this reason, after the List has been created, it is
    * strongly advised to convert it to a C array by calling the function
-   * List.to_array().
+   * `List.to_array()`.
    *
    * \param value: `T` The value of the element to be appended.
    */
@@ -99,22 +144,22 @@ template<class T> class List {
    *
    * \param index: `int` The index of the element to be retrieved. 0 for first.
    * \return value `T` The value of the element at the requested position.
-   * \throws LIST_OUT_OF_BOUNDS_EXCEPTION: Raised if the index is out of bounds.
+   * \throws ListOutOfBoundsException: Raised if the index is out of bounds.
    */
   T get(int index) {
     if (index < 0 || index > size - 1) {
-      throw LIST_OUT_OF_BOUNDS_EXCEPTION;
+      throw ListOutOfBoundsException(index, 0, size - 1);
     }
     current = last;
     for (int i = size - 1; i > index; i--) current = current->p_prev;
     return current->value;
   }
 
-  /*! \brief Get the number of elements in the list.
+  /*! \brief Get the number of elements in the List.
    *
-   * Get the number of elements currently stored in the list.
+   * Get the number of elements currently stored in the List.
    *
-   * \return size `int` The size of the list.
+   * \return size `int` The size of the List.
    */
   int length() {
     return size;
@@ -127,11 +172,11 @@ template<class T> class List {
    *
    * \param index: `int` The index of the element to be set. 0 for first.
    * \param value: `int` The value to store in the pointed element.
-   * \throws LIST_OUT_OF_BOUNDS_EXCEPTION: Raised if the index is out of bounds.
+   * \throws ListOutOfBoundsException: Raised if the index is out of bounds.
    */
   void set(int index, T value) {
     if (index < 0 || index > size - 1) {
-      throw LIST_OUT_OF_BOUNDS_EXCEPTION;
+        throw ListOutOfBoundsException(index, 0, size - 1);
     }
     current = last;
     for (int i = size - 1; i > index; i--) current = current->p_prev;
@@ -145,7 +190,7 @@ template<class T> class List {
    * resulting object is not contiguosly stored in memory. As a result,
    * access to specific elements in the middle of the list is not very
    * effective, because the list needs to be walked every time up to
-   * the desired position. In order to avoid this, List.to_array() makes
+   * the desired position. In order to avoid this, `List.to_array()` makes
    * a conversion from List to C array, returning a contiguous object,
    * where indexed access can be used.
    *
@@ -161,3 +206,5 @@ template<class T> class List {
     return array;
   }
 };
+
+#endif

src/include/Parsers.h

+/*! \file Parsers.h
+ */
+
+#ifndef INCLUDE_PARSERS_H_
+#define INCLUDE_PARSERS_H_
+
+#ifndef FILE_NOT_FOUND_ERROR
+//! Error code if a file is not found.
+#define FILE_NOT_FOUND_ERROR 21
+#endif
+
+/*! \brief Load a text file as a sequence of strings in memory.
+ *
+ * The configuration of the field expansion code in FORTRAN uses
+ * shared memory access and file I/O operations managed by different
+ * functions. Although this approach could be theoretically replicated,
+ * it is more convenient to handle input and output to distinct files
+ * using specific functions. load_file() helps in the task of handling
+ * input such as configuration files or text data structures that need
+ * to be loaded entirely. The function performs a line-by line scan of
+ * the input file and returns an array of strings that can be later
+ * parsed and ingested by the concerned code blocks. An optional pointer
+ * to integer allows the function to keep track of the number of file
+ * lines that were read, if needed.
+ *
+ * \param file_name: `string` The path of the file to be read.
+ * \param count: `int*` Pointer to an integer recording the number of
+ * read lines [OPTIONAL, default=NULL].
+ * \return array_lines `string*` An array of strings, one for each input
+ * file line.
+ */
+std::string *load_file(std::string file_name, int *count);
+
+#endif /* INCLUDE_PARSERS_H_ */

src/include/sph_subs.h

+/*! \file sph_subs.h
+ *
+ * \brief C++ porting of SPH functions and subroutines.
+ *
+ * Remember that FORTRAN passes arguments by reference, so, every time we use
+ * a subroutine call, we need to add a referencing layer to the C++ variable.
+ * All the functions defined below need to be properly documented and ported
+ * to C++.
+ *
+ * Currently, only basic documenting information about functions and parameter
+ * types are given, to avoid doxygen warning messages.
+ */
+
+#ifndef INCLUDE_COMMONS_H_
+#include "Commons.h"
+#endif
+
+#ifndef INCLUDE_SPH_SUBS_H_
+#define INCLUDE_SPH_SUBS_H_
+
+#include <complex>
+
+/*! \brief Conjugate of a double precision complex number
+ *
+ * \param z: `std::complex\<double\>` The input complex number.
+ * \return result: `std::complex\<double\>` The conjugate of the input number.
+ */
+std::complex<double> dconjg(std::complex<double> z) {
+	double zreal = z.real();
+	double zimag = z.imag();
+	return std::complex<double>(zreal, -zimag);
+}
+
+/*! \brief C++ porting of CG1
+ *
+ * \param lmpml: `int`
+ * \param mu: `int`
+ * \param l: `int`
+ * \param m: `int`
+ * \return result: `double`
+ */
+double cg1(int lmpml, int mu, int l, int m) {
+	double result = 0.0;
+	double xd, xn;
+	if (lmpml == -1) { // Interpreted as GOTO 30
+		xd = 2.0 * l * (2 * l - 1);
+		if (mu == -1) {
+			xn = 1.0 * (l - 1 - m) * (l - m);
+		} else if (mu == 0) {
+			xn = 2.0 * (l - m) * (l + m);
+		} else if (mu == 1) {
+			xn = 1.0 * (l - 1 + m) * (l + m);
+		} else {
+			throw 111; // Need an exception for unpredicted indices.
+		}
+		result = sqrt(xn / xd);
+	} else if (lmpml == 0) { // Interpreted as GOTO 5
+		bool goto10 = (m != 0) || (mu != 0);
+		if (!goto10) {
+			result = 0.0;
+			return result;
+		}
+		if (mu != 0) {
+			xd = 2.0 * l * (l + 1);
+			if (mu == -1) {
+				xn = 1.0 * (l - m) * (l + m + 1);
+				result = -sqrt(xn / xd);
+			} else if (mu == 1) { // mu > 0
+				xn = 1.0 * (l + m) * (l - m + 1);
+				result = sqrt(xn / xd);
+			} else {
+				throw 111; // Need an exception for unpredicted indices.
+			}
+		} else { // mu = 0
+			xd = 1.0 * l * (l + 1);
+			xn = -1.0 * m;
+			result = xn / sqrt(xd);
+		}
+	} else if (lmpml == 1) { // Interpreted as GOTO 60
+		xd = 2.0 * (l * 2 + 3) * (l + 1);
+		if (mu == -1) {
+			xn = 1.0 * (l + 1 + m) * (l + 2 + m);
+			result = sqrt(xn / xd);
+		} else if (mu == 0) {
+			xn = 2.0 * (l + 1 - m) * (l + 1 + m);
+			result = -sqrt(xn / xd);
+		} else if (mu == 1) {
+			xn = 1.0 * (l + 1 - m) * (l + 2 - m);
+			result = sqrt(xn / xd);
+		} else { // mu was not recognized.
+			throw 111; // Need an exception for unpredicted indices.
+		}
+	} else { // lmpml was not recognized
+		throw 111; // Need an exception for unpredicted indices.
+	}
+	return result;
+}
+
+/*! \brief C++ porting of APS
+ *
+ * \param zpv: `double ****`
+ * \param li: `int`
+ * \param nsph: `int`
+ * \param c1: `C1 *`
+ * \param sqk: `double`
+ * \param gaps: `double *`
+ */
+void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
+	std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
+	std::complex<double> summ, sume, suem, suee, sum;
+	double half_pi = acos(0.0);
+	double cofs = half_pi * 2.0 / sqk;
+	for (int i40 = 1; i40 <= nsph; i40++) {
+		int iogi = c1->iog[i40 - 1];
+		if (iogi >= i40) {
+			sum = cc0;
+			for (int l30 = 1; l30 <= li; l30++) {
+				int ltpo = l30 + l30 + 1;
+				for (int ilmp = 1; ilmp <= 3; ilmp++) {
+					bool goto30 = (l30 == 1 && ilmp == 1) || (l30 == li && ilmp == 3);
+					if (!goto30) {
+						int lmpml = ilmp - 2;
+						int lmp = l30 + lmpml;
+						double cofl = sqrt(1.0 * (ltpo * (lmp + lmp + 1)));
+						summ = zpv[l30 - 1][ilmp - 1][0][0] /
+								(
+										dconjg(c1->rmi[l30 - 1][i40 - 1]) *
+										c1->rmi[lmp - 1][i40 - 1]
+								);
+						sume = zpv[l30 - 1][ilmp - 1][0][1] /
+								(
+										dconjg(c1->rmi[l30 - 1][i40 - 1]) *
+										c1->rei[lmp - 1][i40 - 1]
+								);
+						suem = zpv[l30 - 1][ilmp - 1][1][0] /
+								(
+										dconjg(c1->rei[l30 - 1][i40 - 1]) *
+										c1->rmi[lmp - 1][i40 - 1]
+								);
+						suee = zpv[l30 - 1][ilmp - 1][1][1] /
+								(
+										dconjg(c1->rei[l30 - 1][i40 - 1]) *
+										c1->rei[lmp - 1][i40 - 1]
+								);
+						sum += (cg1(lmpml, 0, l30, -1) * (summ - sume - suem + suee) +
+								cg1(lmpml, 0, l30, 1) * (summ + sume + suem + suee)) * cofl;
+					}
+				}
+			}
+		}
+		gaps[i40 - 1] = sum.real() * cofs;
+		printf("DEBUG: gaps[%d] = %lE\n", i40, gaps[i40 - 1]);
+	}
+}
+
+/*! \brief C++ porting of DIEL
+ *
+ * \param npntmo: `int`
+ * \param ns: `int`
+ * \param i: `int`
+ * \param ic: `int`
+ * \param vk: `double`
+ * \param c1: `C1 *`
+ * \param c2: `C2 *`
+ */
+void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
+	double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1];
+	double half_step = 0.5 * dif / npntmo;
+	double rr = c1->rc[i - 1][ns - 1];
+	std::complex<double> delta = c2->dc0[ic] - c2->dc0[ic - 1];
+	int kpnt = npntmo + npntmo;
+	c2->ris[kpnt] = c2->dc0[ic];
+	c2->dlri[kpnt] = std::complex<double>(0.0, 0.0);
+	for (int np90 = 0; np90 < kpnt; np90++) {
+		double ff = (rr - c1->rc[i - 1][ns - 1]) / dif;
+		c2->ris[np90] = delta * ff * ff * (-2.0 * ff + 3.0) + c2->dc0[ic - 1];
+		c2->dlri[np90] = 3.0 * delta * ff * (1.0 - ff) / (dif * vk * c2->ris[np90]);
+		rr += half_step;
+	}
+}
+
+/*! \brief C++ porting of ENVJ
+ *
+ * \param n: `int`
+ * \param x: `double`
+ * \return result: `double`
+ */
+double envj(int n, double x) {
+	double result = 0.0;
+	double xn;
+	if (n == 0) {
+		xn = 1.0e-100;
+		result = 0.5 * log10(6.28 * xn) - xn * log10(1.36 * x / xn);
+	} else {
+		result = 0.5 * log10(6.28 * n) - n * log10(1.36 * x / n);
+	}
+	return result;
+}
+
+/*! \brief C++ porting of MSTA1
+ *
+ * \param x: `double`
+ * \param mp: `int`
+ * \return result: `int`
+ */
+int msta1(double x, int mp) {
+	int result = 0;
+	double a0 = x;
+	if (a0 < 0.0) a0 *= -1.0;
+	int n0 = (int)(1.1 * a0) + 1;
+	double f0 = envj(n0, a0) - mp;
+	int n1 = n0 + 5;
+	double f1 = envj(n1, a0) - mp;
+	for (int it10 = 0; it10 < 20; it10++) {
+		int nn = n1 - (int)((n1 - n0) / (1.0 - f0 / f1));
+		double f = envj(nn, a0) - mp;
+		int test_n = nn - n1;
+		if (test_n < 0) test_n *= -1;
+		if (test_n < 1) {
+			return nn;
+		}
+		n0 = n1;
+		f0 = f1;
+		n1 = nn;
+		f1 = f;
+		result = nn;
+	}
+	return result;
+}
+
+/*! \brief C++ porting of MSTA2
+ *
+ * \param x: `double`
+ * \param n: `int`
+ * \param mp: `int`
+ * \return result: `int`
+ */
+int msta2(double x, int n, int mp) {
+	int result = 0;
+	double a0 = x;
+	if (a0 < 0) a0 *= -1.0;
+	double half_mp = 0.5 * mp;
+	double ejn = envj(n, a0);
+	double obj;
+	int n0;
+	if (ejn <= half_mp) {
+		obj = 1.0 * mp;
+		n0 = (int)(1.1 * a0) + 1;
+	} else {
+		obj = half_mp + ejn;
+		n0 = n;
+	}
+	double f0 = envj(n0, a0) - obj;
+	int n1 = n0 + 5;
+	double f1 = envj(n1, a0) - obj;
+	for (int it10 = 0; it10 < 20; it10 ++) {
+		int nn = n1 - (int)((n1 - n0) / (1.0 - f0 / f1));
+		double f = envj(nn, a0) - obj;
+		int test_n = nn - n1;
+		if (test_n < 0) test_n *= -1;
+		if (test_n < 1) return (nn + 10);
+		n0 = n1;
+		f0 = f1;
+		n1 = nn;
+		f1 = f;
+		result = nn + 10;
+	}
+	return result;
+}
+
+/*! \brief C++ porting of CBF
+ *
+ * This is the Complex Bessel Function.
+ *
+ * \param n: `int`
+ * \param z: `complex\<double\>`
+ * \param nm: `int &`
+ * \param csj: Vector of complex.
+ */
+void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
+	/*
+	 * FROM CSPHJY OF LIBRARY specfun
+	 *
+	 * ==========================================================
+	 *   Purpose: Compute spherical Bessel functions j
+	 *   Input :  z --- Complex argument of j
+	 *            n --- Order of j ( n = 0,1,2,... )
+	 *   Output:  csj(n+1) --- j
+	 *            nm --- Highest order computed
+	 *   Routines called:
+	 *            msta1 and msta2 for computing the starting
+	 *            point for backward recurrence
+	 * ==========================================================
+	 */
+	double zz = z.real() * z.real() + z.imag() * z.imag();
+	double a0 = sqrt(zz);
+	nm = n;
+	if (a0 < 1.0e-60) {
+		for (int k = 2; k <= n + 1; k++) {
+			csj[k - 1] = 0.0;
+		}
+		csj[0] = std::complex<double>(1.0, 0.0);
+		return;
+	}
+	csj[0] = std::sin(z) / z;
+	if (n == 0) {
+		return;
+	}
+	csj[1] = (csj[0] -std::cos(z)) / z;
+	if (n == 1) {
+		return;
+	}
+	std::complex<double> csa = csj[0];
+	std::complex<double> csb = csj[1];
+	int m = msta1(a0, 200);
+	if (m < n) nm = m;
+	else m = msta2(a0, n, 15);
+	std::complex<double> cf0 = 0.0;
+	std::complex<double> cf1 = 1.0e-100;
+	std::complex<double> cf, cs;
+	for (int k = m; k >= 0; k--) {
+		cf = (2.0 * k + 3.0) * cf1 / z - cf0;
+		if (k <= nm) csj[k] = cf;
+		cf0 = cf1;
+		cf1 = cf;
+	}
+	double abs_csa = sqrt(csa.real() * csa.real() + csa.imag() * csa.imag());
+	double abs_csb = sqrt(csb.real() * csb.real() + csb.imag() * csb.imag());
+	if (abs_csa > abs_csb) cs = csa / cf;
+	else cs = csb / cf0;
+	for (int k = 0; k <= nm; k++) {
+		csj[k] = cs * csj[k];
+	}
+}
+
+/*! \brief C++ porting of MMULC
+ *
+ * \param vint: Vector of complex.
+ * \param cmullr: `double **`
+ * \param cmul: `double **`
+ */
+void mmulc(std::complex<double> *vint, double **cmullr, double **cmul) {
+	double sm2 = vint[0].real();
+	double s24 = vint[1].real();
+	double d24 = vint[1].imag();
+	double sm4 = vint[5].real();
+	double s23 = vint[8].real();
+	double d32 = vint[8].imag();
+    double s34 = vint[9].real();
+    double d34 = vint[9].imag();
+    double sm3 = vint[10].real();
+    double s31 = vint[11].real();
+    double d31 = vint[11].imag();
+    double s21 = vint[12].real();
+    double d12 = vint[12].imag();
+    double s41 = vint[13].real();
+    double d14 = vint[13].imag();
+    double sm1 = vint[15].real();
+    cmullr[0][0] = sm2;
+    cmullr[0][1] = sm3;
+    cmullr[0][2] = -s23;
+    cmullr[0][3] = -d32;
+    cmullr[1][0] = sm4;
+    cmullr[1][1] = sm1;
+    cmullr[1][2] = -s41;
+    cmullr[1][3] = -d14;
+    cmullr[2][0] = -s24 * 2.0;
+    cmullr[2][1] = -s31 * 2.0;
+    cmullr[2][2] = s21 + s34;
+    cmullr[2][3] = d34 + d12;
+    cmullr[3][0] = -d24 * 2.0;
+    cmullr[3][1] = -d31 * 2.0;
+    cmullr[3][2] = d34 - d12;
+    cmullr[3][3] = s21 - s34;
+    cmul[0][0] = (sm2 + sm3 + sm4 + sm1) * 0.5;
+    cmul[0][1] = (sm2 - sm3 + sm4 - sm1) * 0.5;
+    cmul[0][2] = -s23 - s41;
+    cmul[0][3] = -d32 - d14;
+    cmul[1][0] = (sm2 + sm3 - sm4 - sm1) * 0.5;
+    cmul[1][1] = (sm2 - sm3 - sm4 + sm1) * 0.5;
+    cmul[1][2] = -s23 + s41;
+    cmul[1][3] = -d32 + d14;
+    cmul[2][0] = -s24 - s31;
+    cmul[2][1] = -s24 + s31;
+    cmul[2][2] = s21 + s34;
+    cmul[2][3] = d34 + d12;
+    cmul[3][0] = -d24 - d31;
+    cmul[3][1] = -d24 + d31;
+    cmul[3][2] = d34 - d12;
+    cmul[3][3] = s21 - s34;
+}
+
+/*! \brief C++ porting of ORUNVE
+ *
+ * \param u1: `double *`
+ * \param u2: `double *`
+ * \param u3: `double *`
+ * \param iorth: `int`
+ * \param torth: `double`
+ */
+void orunve( double *u1, double *u2, double *u3, int iorth, double torth) {
+	if (iorth <= 0) {
+		double cp = u1[0] * u2[0] + u1[1] * u2[1] + u1[2] * u2[2];
+		double abs_cp = cp;
+		if (abs_cp < 0.0) abs_cp *= -1.0;
+		if (iorth == 0 || abs_cp >= torth) {
+			double fn = 1.0 / sqrt(1.0 - cp * cp);
+			u3[0] = (u1[1] * u2[2] - u1[2] * u2[1]) * fn;
+			u3[1] = (u1[2] * u2[0] - u1[0] * u2[2]) * fn;
+			u3[2] = (u1[0] * u2[1] - u1[1] * u2[0]) * fn;
+			return;
+		}
+	}
+	u3[0] = u1[1] * u2[2] - u1[2] * u2[1];
+    u3[1] = u1[2] * u2[0] - u1[0] * u2[2];
+    u3[2] = u1[0] * u2[1] - u1[1] * u2[0];
+}
+
+/*! \brief C++ porting of PWMA
+ *
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param ylm: Vector of complex
+ * \param inpol: `int`
+ * \param lw: `int`
+ * \param isq: `int`
+ * \param c1: `C1 *`
+ */
+void pwma(
+		double *up, double *un, std::complex<double> *ylm, int inpol, int lw,
+		int isq, C1 *c1
+) {
+	//std::complex<double> cp1, cm1, cc1, cp2, c2, cc2, uim, cl;
+	double four_pi = 8.0 * acos(0.0);
+	int is = isq;
+	if (isq == -1) is = 0;
+	int ispo = is + 1;
+	int ispt = is + 2;
+	int nlwm = lw * (lw + 2);
+	int nlwmt = nlwm + nlwm;
+	double sqrtwi = 1.0 / sqrt(2.0);
+	std::complex<double> uim(0.0, 1.0);
+	std::complex<double> cm1 = 0.5 * std::complex<double>(up[0], up[1]);
+	std::complex<double> cp1 = 0.5 * std::complex<double>(up[0], -up[1]);
+	double cz1 = up[2];
+	std::complex<double> cm2 = 0.5 * std::complex<double>(un[0], un[1]);
+	std::complex<double> cp2 = 0.5 * std::complex<double>(un[0], -un[1]);
+	double cz2 = un[2];
+	//printf("DEBUG: in PWMA YLM(2) = (%lE, %lE)\n", ylm[1].real(), ylm[1].imag());
+	for (int l20 = 1; l20 <= lw; l20++) {
+		//if (ispo == 1) printf("DEBUG: l20 = %d\n", l20);
+		int lf = l20 + 1;
+		int lftl = lf * l20;
+		double x = 1.0 * lftl;
+		std::complex<double> cl = std::complex<double>(four_pi / sqrt(x), 0.0);
+		for (int powi = 1; powi <= l20; powi++) cl *= uim;
+		int mv = l20 + lf;
+		int m = -lf;
+		for (int mf20 = 1; mf20 <= mv; mf20++) {
+			m += 1;
+			int k = lftl + m;
+			x = 1.0 * (lftl - m * (m + 1));
+			double cp = sqrt(x);
+			x = 1.0 * (lftl - m * (m - 1));
+			double cm = sqrt(x);
+			double cz = 1.0 * m;
+			c1->w[k - 1][ispo - 1] = dconjg(
+					cp1 * cp * ylm[k + 1] +
+					cm1 * cm * ylm[k - 1] +
+					cz1 * cz * ylm[k]
+			) * cl;
+			//if (ispo == 1) printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k, ispo, c1->w[k - 1][ispo - 1].real(), c1->w[k - 1][ispo - 1].imag());
+			c1->w[k - 1][ispt - 1] = dconjg(
+					cp2 * cp * ylm[k + 1] +
+					cm2 * cm * ylm[k - 1] +
+					cz2 * cz * ylm[k]
+			) * cl;
+			//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k, ispt, c1->w[k - 1][ispt - 1].real(), c1->w[k - 1][ispt - 1].imag());
+		}
+	}
+	for (int k30 = 1; k30 <= nlwm; k30++) {
+		int i = k30 + nlwm;
+		c1->w[i - 1][ispo - 1] = uim * c1->w[k30 - 1][ispt - 1];
+		//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispo, c1->w[i - 1][ispo - 1].real(), c1->w[i - 1][ispo - 1].imag());
+		c1->w[i - 1][ispt - 1] = -uim * c1->w[k30 - 1][ispo - 1];
+		//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispt, c1->w[i - 1][ispt - 1].real(), c1->w[i - 1][ispt - 1].imag());
+	}
+	if (inpol != 0) {
+		for (int k40 = 1; k40 <= nlwm; k40++) {
+			int i = k40 + nlwm;
+			std::complex<double> cc1 = sqrtwi * (c1->w[k40 - 1][ispo - 1] + uim * c1->w[k40 - 1][ispt - 1]);
+			std::complex<double> cc2 = sqrtwi * (c1->w[k40 - 1][ispo - 1] - uim * c1->w[k40 - 1][ispt - 1]);
+			c1->w[k40 - 1][ispo - 1] = cc2;
+			//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k40, ispo, c1->w[k40 - 1][ispo - 1].real(), c1->w[k40 - 1][ispo - 1].imag());
+			c1->w[i - 1][ispo - 1] = -cc2;
+			//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispo, c1->w[i - 1][ispo - 1].real(), c1->w[i - 1][ispo - 1].imag());
+			c1->w[k40 - 1][ispt - 1] = cc1;
+			//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k40, ispt, c1->w[k40 - 1][ispt - 1].real(), c1->w[k40 - 1][ispt - 1].imag());
+			c1->w[i - 1][ispt - 1] = cc1;
+			//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispt, c1->w[i - 1][ispt - 1].real(), c1->w[i - 1][ispt - 1].imag());
+		}
+	} else {
+		if (isq == 0) {
+			return;
+		}
+	}
+	if (isq != 0) {
+		for (int i50 = 1; i50 <= 2; i50++) {
+			int ipt = i50 + 2;
+			int ipis = i50 + is;
+			for (int k50 = 1; k50 <= nlwmt; k50++) {
+				c1->w[k50 - 1][ipt - 1] = dconjg(c1->w[k50 - 1][ipis - 1]);
+				//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k50, ipt, c1->w[k50 - 1][ipt - 1].real(), c1->w[k50 - 1][ipt - 1].imag());
+			}
+		}
+	}
+	//printf("DEBUG: out PWMA W(1,1) = (%lE, %lE)\n", c1->w[0][0].real(), c1->w[0][0].imag());
+}
+
+/*! \brief C++ porting of RABAS
+ *
+ * \param inpol: `int`
+ * \param li: `int`
+ * \param nsph: `int`
+ * \param c1: `C1 *`
+ * \param tqse: Matrix of complex.
+ * \param tqspe: Matrix of complex.
+ * \param tqss: Matrix of complex.
+ * \param tqsps: Matrix of complex.
+ */
+void rabas(
+			int inpol, int li, int nsph, C1 *c1, double **tqse, std::complex<double> **tqspe,
+			double **tqss, std::complex<double> **tqsps
+) {
+	std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
+	std::complex<double> uim = std::complex<double>(0.0, 1.0);
+	double two_pi = 4.0 * acos(0.0);
+	for (int i80 = 1; i80 <= nsph; i80++) {
+		if(c1->iog[i80 - 1] >= i80) {
+			tqse[0][i80 - 1] = 0.0;
+			tqse[1][i80 - 1] = 0.0;
+			tqspe[0][i80 - 1] = cc0;
+			tqspe[1][i80 - 1] = cc0;
+			tqss[0][i80 - 1] = 0.0;
+			tqss[1][i80 - 1] = 0.0;
+			tqsps[0][i80 - 1] = cc0;
+			tqsps[1][i80 - 1] = cc0;
+			for (int l70 = 1; l70 <= li; l70++) {
+				double fl = 1.0 * l70 + l70 + 1;
+				std::complex<double> rm = 1.0 / c1->rmi[l70 - 1][i80 - 1];
+				double rmm = (rm * dconjg(rm)).real();
+				std::complex<double> re = 1.0 / c1->rei[l70 - 1][i80 - 1];
+				double rem = (re * dconjg(re)).real();
+				if (inpol == 0) {
+					std::complex<double> pce = ((rm + re) * uim) * fl;
+					std::complex<double> pcs = ((rmm + rem) * fl) * uim;
+					tqspe[0][i80 - 1] -= pce;
+					tqspe[1][i80 - 1] += pce;
+					tqsps[0][i80 - 1] -= pcs;
+					tqsps[1][i80 - 1] += pcs;
+				} else {
+					double ce = (rm + re).real() * fl;
+					double cs = (rmm + rem) * fl;
+					tqse[0][i80 - 1] -= ce;
+					tqse[1][i80 - 1] += ce;
+					tqss[0][i80 - 1] -= cs;
+					tqss[1][i80 - 1] += cs;
+				}
+			}
+			if (inpol == 0) {
+				tqspe[0][i80 - 1] *= two_pi;
+				tqspe[1][i80 - 1] *= two_pi;
+				tqsps[0][i80 - 1] *= two_pi;
+				tqsps[1][i80 - 1] *= two_pi;
+				printf("DEBUG: TQSPE(1, %d ) = ( %lE, %lE)\n", i80, tqspe[0][i80 - 1].real(), tqspe[0][i80 - 1].imag());
+				printf("DEBUG: TQSPE(2, %d ) = ( %lE, %lE)\n", i80, tqspe[1][i80 - 1].real(), tqspe[1][i80 - 1].imag());
+				printf("DEBUG: TQSPS(1, %d ) = ( %lE, %lE)\n", i80, tqsps[0][i80 - 1].real(), tqsps[0][i80 - 1].imag());
+				printf("DEBUG: TQSPS(2, %d ) = ( %lE, %lE)\n", i80, tqsps[1][i80 - 1].real(), tqsps[1][i80 - 1].imag());
+			} else {
+				tqse[0][i80 - 1] *= two_pi;
+				tqse[1][i80 - 1] *= two_pi;
+				tqss[0][i80 - 1] *= two_pi;
+				tqss[1][i80 - 1] *= two_pi;
+				printf("DEBUG: TQSE(1, %d ) = %lE)\n", i80, tqse[0][i80 - 1]);
+				printf("DEBUG: TQSE(2, %d ) = %lE)\n", i80, tqse[1][i80 - 1]);
+				printf("DEBUG: TQSS(1, %d ) = %lE)\n", i80, tqss[0][i80 - 1]);
+				printf("DEBUG: TQSS(2, %d ) = %lE)\n", i80, tqss[1][i80 - 1]);
+			}
+		}
+	}
+}
+
+/*! \brief C++ porting of RBF
+ *
+ * This is the Real Bessel Function.
+ *
+ * \param n: `int`
+ * \param x: `double`
+ * \param nm: `int &`
+ * \param sj: `double[]`
+ */
+void rbf(int n, double x, int &nm, double sj[]) {
+	/*
+	 * FROM SPHJ OF LIBRARY specfun
+	 *
+	 * ==========================================================
+	 *   Purpose: Compute spherical Bessel functions j
+	 *   Input :  x --- Argument of j
+	 *            n --- Order of j ( n = 0,1,2,... )
+	 *   Output:  sj(n+1) --- j
+	 *            nm --- Highest order computed
+	 *   Routines called:
+	 *            msta1 and msta2 for computing the starting
+	 *            point for backward recurrence
+	 * ==========================================================
+	 */
+	double a0 = x;
+	if (a0 < 0.0) a0 *= -1.0;
+	nm = n;
+	if (a0 < 1.0e-60) {
+		for (int k = 2; k <= n + 1; k++)
+			sj[k - 1] = 0.0;
+		sj[0] = 1.0;
+		return;
+	}
+	sj[0] = sin(x) / x;
+	if (n == 0) {
+		return;
+	}
+	sj[1] = (sj[0] - cos(x)) / x;
+	if (n == 1) {
+		return;
+	}
+	double sa = sj[0];
+	double sb = sj[1];
+	int m = msta1(a0, 200);
+	if (m < n) nm = m;
+	else m = msta2(a0, n, 15);
+	double f0 = 0.0;
+	double f1 = 1.0e-100;
+	double f;
+	for (int k = m; k >= 0; k--) {
+		f = (2.0 * k +3.0) * f1 / x - f0;
+		if (k <= nm) sj[k] = f;
+		f0 = f1;
+		f1 = f;
+	}
+	double cs;
+	double abs_sa = sa;
+	if (abs_sa < 0.0) abs_sa *= -1.0;
+	double abs_sb = sb;
+	if (abs_sb < 0.0) abs_sb *= -1.0;
+	if (abs_sa > abs_sb) cs = sa / f;
+	else cs = sb / f0;
+	for (int k = 0; k <= nm; k++) {
+		sj[k] = cs * sj[k];
+	}
+}
+
+/*! \brief C++ porting of RKC
+ *
+ * \param npntmo: `int`
+ * \param step: `double`
+ * \param dcc: `complex\<double\>`
+ * \param x: `double &`
+ * \param lpo: `int`
+ * \param y1: `complex\<double\> &`
+ * \param y2: `complex\<double\> &`
+ * \param dy1: `complex\<double\> &`
+ * \param dy2: `complex\<double\> &`
+ */
+void rkc(
+		int npntmo, double step, std::complex<double> dcc, double &x, int lpo,
+		std::complex<double> &y1, std::complex<double> &y2, std::complex<double> &dy1,
+		std::complex<double> &dy2
+) {
+	std::complex<double> cy1, cdy1, c11, cy23, yc2, c12, c13;
+	std::complex<double> cy4, yy, c14, c21, c22, c23, c24;
+	double half_step = 0.5 * step;
+	double cl = 1.0 * lpo * (lpo - 1);
+	for (int ipnt60 = 1; ipnt60 <= npntmo; ipnt60++) {
+		cy1 = cl / (x * x) - dcc;
+		cdy1 = -2.0 / x;
+		c11 = (cy1 * y1 + cdy1 * dy1) * step;
+		double xh = x + half_step;
+		cy23 = cl / (xh * xh) - dcc;
+		double cdy23 = -2.0 / xh;
+		yc2 = y1 + dy1 * half_step;
+		c12 = (cy23 * yc2 + cdy23 * (dy1 + 0.5 * c11)) * step;
+		c13 = (cy23 * (yc2 + 0.25 * c11 * step) + cdy23 * (dy1 + 0.5 * c12)) * step;
+		double xn = x + step;
+		cy4 = cl / (xn * xn) - dcc;
+		double cdy4 = -2.0 / xn;
+		yy = y1 + dy1 * step;
+		c14 = (cy4 * (yy + 0.5 * c12 * step) + cdy4 * (dy1 + c13)) * step;
+		y1 = yy + (c11 + c12 + c13) * step / 6.0;
+		dy1 += (0.5 * c11 + c12 + c13 + 0.5 * c14) / 3.0;
+		c21 = (cy1 * y2 + cdy1 * dy2) * step;
+		yc2 = y2 + dy2 * half_step;
+		c22 = (cy23 * yc2 + cdy23 * (dy2 + 0.5 * c21)) * step;
+		c23= (cy23 * (yc2 + 0.25 * c21 * step) + cdy23 * (dy2 + 0.5 * c22)) * step;
+		yy = y2 + dy2 * step;
+		c24 = (cy4 * (yc2 + 0.5 * c22 * step) + cdy4 * (dy2 + c23)) * step;
+		y2 = yy + (c21 + c22 + c23) * step / 6.0;
+		dy2 += (0.5 * c21 + c22 + c23 + 0.5 * c24) / 3.0;
+		x = xn;
+	}
+}
+
+/*! \brief C++ porting of RKT
+ *
+ * \param npntmo: `int`
+ * \param step: `double`
+ * \param x: `double &`
+ * \param lpo: `int`
+ * \param y1: `complex\<double\> &`
+ * \param y2: `complex\<double\> &`
+ * \param dy1: `complex\<double\> &`
+ * \param dy2: `complex\<double\> &`
+ * \param c2: `C2 *`
+ */
+void rkt(
+		int npntmo, double step, double &x, int lpo, std::complex<double> &y1,
+		std::complex<double> &y2, std::complex<double> &dy1, std::complex<double> &dy2,
+		C2 *c2
+) {
+	std::complex<double> cy1, cdy1, c11, cy23, cdy23, yc2, c12, c13;
+	std::complex<double> cy4, cdy4, yy, c14, c21, c22, c23, c24;
+	double half_step = 0.5 * step;
+	double cl = 1.0 * lpo * (lpo - 1);
+	for (int ipnt60 = 1; ipnt60 <= npntmo; ipnt60++) {
+		int jpnt = ipnt60 + ipnt60 - 1;
+		cy1 = cl / (x * x) - c2->ris[jpnt - 1];
+		cdy1 = -2.0 / x;
+		c11 = (cy1 * y1 + cdy1 * dy1) * step;
+		double xh = x + half_step;
+		int jpntpo = jpnt + 1;
+		cy23 = cl / (xh * xh) - c2->ris[jpntpo - 1];
+		cdy23 = -2.0 / xh;
+		yc2 = y1 + dy1 * half_step;
+		c12 = (cy23 * yc2 + cdy23 * (dy1 + 0.5 * c11)) * step;
+		c13= (cy23 * (yc2 + 0.25 * c11 *step) + cdy23 * (dy1 + 0.5 * c12)) * step;
+		double xn = x + step;
+		int jpntpt = jpnt + 2;
+		cy4 = cl / (xn * xn) - c2->ris[jpntpt - 1];
+		cdy4 = -2.0 / xn;
+		yy = y1 + dy1 * step;
+		c14 = (cy4 * (yy + 0.5 * c12 * step) + cdy4 * (dy1 + c13)) * step;
+		y1= yy + (c11 + c12 + c13) * step / 6.0;
+		dy1 += (0.5 * c11 + c12 + c13 + 0.5 * c14) /3.0;
+		cy1 -= cdy1 * c2->dlri[jpnt - 1];
+		cdy1 += 2.0 * c2->dlri[jpnt - 1];
+		c21 = (cy1 * y2 + cdy1 * dy2) * step;
+		cy23 -= cdy23 * c2->dlri[jpntpo - 1];
+		cdy23 += 2.0 * c2->dlri[jpntpo - 1];
+		yc2 = y2 + dy2 * half_step;
+		c22 = (cy23 * yc2 + cdy23 * (dy2 + 0.5 * c21)) * step;
+		c23 = (cy23 * (yc2 + 0.25 * c21 * step) + cdy23 * (dy2 + 0.5 * c22)) * step;
+		cy4 -= cdy4 * c2->dlri[jpntpt - 1];
+		cdy4 += 2.0 * c2->dlri[jpntpt - 1];
+		yy = y2 + dy2 * step;
+		c24 = (cy4 * (yc2 + 0.5 * c22 * step) + cdy4 * (dy2 + c23)) * step;
+		y2 = yy + (c21 + c22 + c23) * step / 6.0;
+		dy2 += (0.5 * c21 + c22 + c23 + 0.5 * c24) / 3.0;
+		x = xn;
+	}
+}
+
+/*! \brief C++ porting of RNF
+ *
+ * This is a real spherical Bessel function.
+ *
+ * \param n: `int`
+ * \param x: `double`
+ * \param nm: `int &`
+ * \param sy: `double[]`
+ */
+void rnf(int n, double x, int &nm, double sy[]) {
+	/*
+	 * FROM SPHJY OF LIBRARY specfun
+	 *
+	 * ==========================================================
+	 *   Purpose: Compute spherical Bessel functions y
+	 *   Input :  x --- Argument of y ( x > 0 )
+	 *            n --- Order of y ( n = 0,1,2,... )
+	 *   Output:  sy(n+1) --- y
+	 *            nm --- Highest order computed
+	 * ==========================================================
+	 */
+	if (x < 1.0e-60) {
+		for (int k = 0; k <= n; k++)
+			sy[k] = -1.0e300;
+		return;
+	}
+	sy[0] = -1.0 * cos(x) / x;
+	if (n == 0) {
+		return;
+	}
+	sy[1] = (sy[0] - sin(x)) / x;
+	if (n == 1) {
+		return;
+	}
+	double f0 = sy[0];
+	double f1 = sy[1];
+	double f;
+	for (int k = 2; k <= n; k++) {
+		f = (2.0 * k - 1.0) * f1 / x - f0;
+		sy[k] = f;
+		double abs_f = f;
+		if (abs_f < 0.0) abs_f *= -1.0;
+		if (abs_f >= 1.0e300) {
+			nm = k;
+			break;
+		}
+		f0 = f1;
+		f1 = f;
+		nm = k;
+	}
+	return;
+}
+
+/*! \brief C++ porting of SSCR0
+ *
+ * \param tfsas: `complex<double> &`
+ * \param nsph: `int`
+ * \param lm: `int`
+ * \param vk: `double`
+ * \param exri: `double`
+ * \param c1: `C1 *`
+ */
+void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
+	std::complex<double> sum21, rm, re, csam;
+	std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
+	double exdc = exri * exri;
+	double ccs = 4.0 * acos(0.0) / (vk * vk);
+	double cccs = ccs / exdc;
+	csam = -(ccs / (exri * vk)) * std::complex<double>(0.0, 0.5);
+	tfsas = cc0;
+	for (int i12 = 1; i12 <= nsph; i12++) {
+		int iogi = c1->iog[i12 - 1];
+		if (iogi >= i12) {
+			double sums = 0.0;
+			std::complex<double> sum21 = cc0;
+			for (int l10 = 1; l10 <= lm; l10++) {
+				double fl = 1.0 + l10 + l10;
+				rm = 1.0 / c1->rmi[l10 - 1][i12 - 1];
+				re = 1.0 / c1->rei[l10 - 1][i12 - 1];
+				std::complex<double> rm_cnjg = std::complex<double>(rm.real(), -rm.imag());
+				std::complex<double> re_cnjg = std::complex<double>(re.real(), -re.imag());
+				sums += (rm_cnjg * rm + re_cnjg * re).real() * fl;
+				sum21 += (rm + re) * fl;
+			}
+			sum21 *= -1.0;
+			double scasec = cccs * sums;
+			double extsec = -cccs * sum21.real();
+			double abssec = extsec - scasec;
+			c1->sscs[i12 - 1] = scasec;
+			c1->sexs[i12 - 1] = extsec;
+			c1->sabs[i12 - 1] = abssec;
+			double gcss = c1-> gcsv[i12 - 1];
+			c1->sqscs[i12 - 1] = scasec / gcss;
+			c1->sqexs[i12 - 1] = extsec / gcss;
+			c1->sqabs[i12 - 1] = abssec / gcss;
+			c1->fsas[i12 - 1] = sum21 * csam;
+			//printf("DEBUG: FSAS( %d ) = (%lE,%lE)\n", i12, c1->fsas[i12 - 1].real(), c1->fsas[i12 - 1].imag());
+		}
+		tfsas += c1->fsas[iogi - 1];
+	}
+}
+
+/*! \brief C++ porting of SSCR2
+ *
+ * \param nsph: `int`
+ * \param lm: `int`
+ * \param vk: `double`
+ * \param exri: `double`
+ * \param c1: `C1 *`
+ */
+void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
+	std::complex<double> s11, s21, s12, s22, rm, re, csam, cc;
+	std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
+	double ccs = 1.0 / (vk * vk);
+	csam = -(ccs / (exri * vk)) * std::complex<double>(0.0, 0.5);
+	double pigfsq = 64.0 * acos(0.0) * acos(0.0);
+	double cfsq = 4.0 / (pigfsq * ccs * ccs);
+	int nlmm = lm * (lm + 2);
+	//printf("DEBUG: in SSCR2 W(1,1) = (%lE,%lE)\n", c1->w[0][0].real(), c1->w[0][0].imag());
+	for (int i14 = 1; i14 <= nsph; i14++) {
+		int iogi = c1->iog[i14 - 1];
+		if (iogi >= i14) {
+			int k = 0;
+			s11 = cc0;
+			s21 = cc0;
+			s12 = cc0;
+			s22 = cc0;
+			for (int l10 = 1; l10 <= lm; l10++) {
+				rm = 1.0 / c1->rmi[l10 - 1][i14 - 1];
+				re = 1.0 / c1->rei[l10 - 1][i14 - 1];
+				//printf("DEBUG: rm = (%lE,%lE)\n", rm.real(), rm.imag());
+				//printf("DEBUG: re = (%lE,%lE)\n", re.real(), re.imag());
+				int ltpo = l10 + l10 + 1;
+				for (int im10 = 1; im10 <= ltpo; im10++) {
+					k += 1;
+					int ke = k + nlmm;
+					//printf("DEBUG: W( %d, 3) = (%lE,%lE)\n", k, c1->w[k - 1][2].real(), c1->w[k - 1][2].imag());
+					//printf("DEBUG: W( %d, 1) = (%lE,%lE)\n", k, c1->w[k - 1][0].real(), c1->w[k - 1][0].imag());
+					//printf("DEBUG: W( %d, 3) = (%lE,%lE)\n", ke, c1->w[ke - 1][2].real(), c1->w[ke - 1][2].imag());
+					//printf("DEBUG: W( %d, 1) = (%lE,%lE)\n", ke, c1->w[ke - 1][0].real(), c1->w[ke - 1][0].imag());
+					s11 = s11 - c1->w[k - 1][2] * c1->w[k - 1][0] * rm - c1->w[ke - 1][2] * c1->w[ke - 1][0] * re;
+					//printf("DEBUG: W( %d, 4) = (%lE,%lE)\n", k, c1->w[k - 1][3].real(), c1->w[k - 1][3].imag());
+					//printf("DEBUG: W( %d, 4) = (%lE,%lE)\n", ke, c1->w[ke - 1][3].real(), c1->w[ke - 1][3].imag());
+					s21 = s21 - c1->w[k - 1][3] * c1->w[k - 1][0] * rm - c1->w[ke - 1][3] * c1->w[ke - 1][0] * re;
+					//printf("DEBUG: W( %d, 2) = (%lE,%lE)\n", k, c1->w[k - 1][1].real(), c1->w[k - 1][1].imag());
+					//printf("DEBUG: W( %d, 2) = (%lE,%lE)\n", ke, c1->w[ke - 1][1].real(), c1->w[ke - 1][1].imag());
+					s12 = s12 - c1->w[k - 1][2] * c1->w[k - 1][1] * rm - c1->w[ke - 1][2] * c1->w[ke - 1][1] * re;
+					s22 = s22 - c1->w[k - 1][3] * c1->w[k - 1][1] * rm - c1->w[ke - 1][3] * c1->w[ke - 1][1] * re;
+				}
+			}
+			c1->sas[i14 - 1][0][0] = s11 * csam;
+			c1->sas[i14 - 1][1][0] = s21 * csam;
+			c1->sas[i14 - 1][0][1] = s12 * csam;
+			c1->sas[i14 - 1][1][1] = s22 * csam;
+		}
+	} // loop i14
+	for (int i24 = 1; i24 <= nsph; i24++) {
+		int iogi = c1->iog[i24 - 1];
+		if (iogi >= i24) {
+			int j = 0;
+			for (int ipo1 = 0; ipo1 < 2; ipo1++) {
+				for (int jpo1 = 0; jpo1 < 2; jpo1++) {
+					std::complex<double> cc = dconjg(c1->sas[i24 - 1][jpo1][ipo1]);
+					for (int ipo2 = 0; ipo2 < 2; ipo2++) {
+						for (int jpo2 = 0; jpo2 < 2; jpo2++) {
+							j += 1;
+							c1->vints[i24 - 1][j - 1] = c1->sas[i24 - 1][jpo2][ipo2] * cc * cfsq;
+						}
+					}
+				}
+			}
+		}
+	}
+}
+
+/*! \brief C++ porting of SPHAR
+ *
+ * \param cosrth: `double`
+ * \param sinrth: `double`
+ * \param cosrph: `double`
+ * \param sinrph: `double`
+ * \param ll: `int`
+ * \param ylm: Vector of complex.
+ */
+void sphar(
+		double cosrth, double sinrth, double cosrph, double sinrph,
+		int ll, std::complex<double> *ylm
+) {
+	double sinrmp[40], cosrmp[40], plegn[861];
+	double four_pi = 8.0 * acos(0.0);
+	double pi4irs = 1.0 / sqrt(four_pi);
+	double x = cosrth;
+	double y = sinrth;
+	//printf("DEBUG: X = %lE\n", x);
+	if (y < 0.0) y *= -1.0;
+	//printf("DEBUG: Y = %lE\n", y);
+	double cllmo = 3.0;
+	double cll = 1.5;
+	double ytol = y;
+	plegn[0] = 1.0;
+	//printf("DEBUG: PLEGN( %d ) = %lE\n", 1, plegn[0]);
+	plegn[1] = x * sqrt(cllmo);
+	//printf("DEBUG: PLEGN( %d ) = %lE\n", 2, plegn[1]);
+	plegn[2] = ytol * sqrt(cll);
+	//printf("DEBUG: PLEGN( %d ) = %lE\n", 3, plegn[2]);
+	sinrmp[0] = sinrph;
+	cosrmp[0] = cosrph;
+	if (ll >= 2) {
+		int k = 3;
+		for (int l20 = 2; l20 <= ll; l20++) {
+			int lmo = l20 - 1;
+			int ltpo = l20 + l20 + 1;
+			int ltmo = ltpo - 2;
+			int lts = ltpo * ltmo;
+			double cn = 1.0 * lts;
+			for (int mpo10 = 1; mpo10 <= lmo; mpo10++) {
+				int m = mpo10 - 1;
+				int mpopk = mpo10 + k;
+				int ls = (l20 + m) * (l20 - m);
+				double cd = 1.0 * ls;
+				double cnm = 1.0 * ltpo * (lmo + m) * (l20 - mpo10);
+				double cdm = 1.0 * ls * (ltmo - 2);
+				plegn[mpopk - 1] = plegn[mpopk - l20 - 1] * x * sqrt(cn / cd) -
+						plegn[mpopk - ltmo - 1] * sqrt(cnm / cdm);
+				//printf("DEBUG: PLEGN( %d ) = %lE\n", mpopk, plegn[mpopk - 1]);
+			}
+			int lpk = l20 + k;
+			double cltpo = 1.0 * ltpo;
+			plegn[lpk - 1] = plegn[k - 1] * x * sqrt(cltpo);
+			//printf("DEBUG: PLEGN( %d ) = %lE\n", lpk, plegn[lpk - 1]);
+			k = lpk + 1;
+			double clt = 1.0 * (ltpo - 1);
+			cll *= (cltpo / clt);
+			ytol *= y;
+			plegn[k - 1] = ytol * sqrt(cll);
+			//printf("DEBUG: PLEGN( %d ) = %lE\n", k, plegn[k - 1]);
+			sinrmp[l20 - 1] = sinrph * cosrmp[lmo - 1] + cosrph * sinrmp[lmo - 1];
+			cosrmp[l20 - 1] = cosrph * cosrmp[lmo - 1] - sinrph * sinrmp[lmo - 1];
+		} // end l20 loop
+	}
+	// label 30
+	int l = 0;
+	int m, k, l0y, l0p, lmy, lmp;
+	double save;
+label40:
+	m = 0;
+	k = l * (l + 1);
+	l0y = k + 1;
+	l0p = k / 2 + 1;
+	ylm[l0y - 1] = pi4irs * plegn[l0p - 1];
+	//printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", l0y, ylm[l0y - 1].real(), ylm[l0y - 1].imag());
+	goto label45;
+label44:
+	lmp = l0p + m;
+	save = pi4irs * plegn[lmp - 1];
+	lmy = l0y + m;
+	ylm[lmy - 1] = save * std::complex<double>(cosrmp[m - 1], sinrmp[m - 1]);
+	if (m % 2 != 0) ylm[lmy - 1] *= -1.0;
+	//printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", lmy, ylm[lmy - 1].real(), ylm[lmy - 1].imag());
+	lmy = l0y - m;
+	ylm[lmy - 1] = save * std::complex<double>(cosrmp[m - 1], -sinrmp[m - 1]);
+	//printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", lmy, ylm[lmy - 1].real(), ylm[lmy - 1].imag());
+label45:
+	if (m >= l) goto label47;
+	m += 1;
+	goto label44;
+label47:
+	if (l >= ll) return;
+	l += 1;
+	goto label40;
+}
+
+/*! \brief C++ porting of THDPS
+ *
+ * \param lm: `int`
+ * \param zpv: `double ****`
+ */
+void thdps(int lm, double ****zpv) {
+	// WARNING: unclear nested loop in THDPS
+	// The optimized interpretation was adopted here.
+	for (int l10 = 1; l10 <= lm; l10++) {
+		for (int ilmp = 1; ilmp <= 3; ilmp++) {
+			zpv[l10 - 1][ilmp - 1][0][0] = 0.0;
+			zpv[l10 - 1][ilmp - 1][0][1] = 0.0;
+			zpv[l10 - 1][ilmp - 1][1][0] = 0.0;
+			zpv[l10 - 1][ilmp - 1][1][1] = 0.0;
+		}
+	}
+	for (int l15 = 1; l15 <= lm; l15++) {
+		double xd = 1.0 * l15 * (l15 + 1);
+		double zp = -1.0 / sqrt(xd);
+		zpv[l15 - 1][1][0][1] = zp;
+		zpv[l15 - 1][1][1][0] = zp;
+	}
+	if (lm != 1) {
+		for (int l20 = 2; l20 <= lm; l20++) {
+			double xn = 1.0 * (l20 - 1) * (l20 + 1);
+			double xd = 1.0 * l20 * (l20 + l20 + 1);
+			double zp = sqrt(xn / xd);
+			zpv[l20 - 1][0][0][0] = zp;
+			zpv[l20 - 1][0][1][1] = zp;
+		}
+		int lmmo = lm - 1;
+		for (int l25 = 1; l25 <= lmmo; l25++) {
+			double xn = 1.0 * l25 * (l25 + 2);
+			double xd = (l25 + 1) * (l25 + l25 + 1);
+			double zp = -1.0 * sqrt(xn / xd);
+			zpv[l25 - 1][2][0][0] = zp;
+			zpv[l25 - 1][2][1][1] = zp;
+		}
+	}
+}
+
+/*! \brief C++ porting of UPVMP
+ *
+ * \param thd: `double`
+ * \param phd: `double`
+ * \param icspnv: `int`
+ * \param cost: `double`
+ * \param sint: `double`
+ * \param cosp: `double`
+ * \param sinp: `double`
+ * \param u: `double *`
+ * \param up: `double *`
+ * \param un: `double *`
+ */
+void upvmp(
+		double thd, double phd, int icspnv, double &cost, double &sint,
+		double &cosp, double &sinp, double *u, double *up, double *un
+) {
+	double half_pi = acos(0.0);
+	double rdr = half_pi / 90.0;
+	double th = thd * rdr;
+	double ph = phd * rdr;
+	cost = cos(th);
+	sint = sin(th);
+	cosp = cos(ph);
+	sinp = sin(ph);
+	//printf("DEBUG: cost = %lE\n", cost);
+	//printf("DEBUG: sint = %lE\n", sint);
+	//printf("DEBUG: cosp = %lE\n", cosp);
+	//printf("DEBUG: sinp = %lE\n", sinp);
+	u[0] = cosp * sint;
+	u[1] = sinp * sint;
+	u[2] = cost;
+	up[0] = cosp * cost;
+	up[1] = sinp * cost;
+	up[2] = -sint;
+	un[0] = -sinp;
+	un[1] = cosp;
+	un[2] = 0.0;
+	if (icspnv != 0) {
+		up[0] *= -1.0;
+		up[1] *= -1.0;
+		up[2] *= -1.0;
+		un[0] *= -1.0;
+		un[1] *= -1.0;
+	}
+}
+
+/*! \brief C++ porting of UPVSP
+ *
+ * \param u: `double *`
+ * \param upmp: `double *`
+ * \param unmp: `double *`
+ * \param us: `double *`
+ * \param upsmp: `double *`
+ * \param unsmp: `double *`
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param ups: `double *`
+ * \param uns: `double *`
+ * \param duk: `double *`
+ * \param isq: `int &`
+ * \param ibf: `int &`
+ * \param scand: `double &`
+ * \param cfmp: `double &`
+ * \param sfmp: `double &`
+ * \param cfsp: `double &`
+ * \param sfsp: `double &`
+ */
+void upvsp(
+		double *u, double *upmp, double *unmp, double *us, double *upsmp, double *unsmp,
+		double *up, double *un, double *ups, double *uns, double *duk, int &isq,
+		int &ibf, double &scand, double &cfmp, double &sfmp, double &cfsp, double &sfsp
+) {
+	double rdr = acos(0.0) / 90.0;
+	double small = 1.0e-6;
+	isq = 0;
+	scand = u[0] * us[0] + u[1] * us[1] + u[2] * us[2];
+	double abs_scand = scand - 1.0;
+	if (abs_scand < 0.0) abs_scand *= -1.0;
+	if (abs_scand >= small) {
+		abs_scand = scand + 1.0;
+		if (abs_scand < 0.0) abs_scand *= -1.0;
+		if (abs_scand >= small) {
+			scand = acos(scand) / rdr;
+			duk[0] = u[0] - us[0];
+			duk[1] = u[1] - us[1];
+			duk[2] = u[2] - us[2];
+			ibf = 0;
+		} else {
+			// label 15
+			scand = 180.0;
+			duk[0] = 2.0 * u[0];
+			duk[1] = 2.0 * u[1];
+			duk[2] = 2.0 * u[2];
+			ibf = 1;
+			ups[0] = -upsmp[0];
+			ups[1] = -upsmp[1];
+			ups[2] = -upsmp[2];
+			uns[0] = -unsmp[0];
+			uns[1] = -unsmp[1];
+			uns[2] = -unsmp[2];
+		}
+	} else {
+		// label 10
+		scand = 0.0;
+		duk[0] = 0.0;
+		duk[1] = 0.0;
+		duk[2] = 0.0;
+		ibf = -1;
+		isq = -1;
+		ups[0] = upsmp[0];
+		ups[1] = upsmp[1];
+		ups[2] = upsmp[2];
+		uns[0] = unsmp[0];
+		uns[1] = unsmp[1];
+		uns[2] = unsmp[2];
+	}
+	if (ibf == -1 || ibf == 1) { // label 20
+		up[0] = upmp[0];
+		up[1] = upmp[1];
+		up[2] = upmp[2];
+		un[0] = unmp[0];
+		un[1] = unmp[1];
+		un[2] = unmp[2];
+	} else { // label 25
+		orunve(u, us, un, -1, small);
+		uns[0] = un[0];
+		uns[1] = un[1];
+		uns[2] = un[2];
+		orunve(un, u, up, 1, small);
+		orunve(uns, us, ups, 1, small);
+	}
+	// label 85
+	cfmp = upmp[0] * up[0] + upmp[1] * up[1] + upmp[2] * up[2];
+	sfmp = unmp[0] * up[0] + unmp[1] * up[1] + unmp[2] * up[2];
+    cfsp = ups[0] * upsmp[0] + ups[1] * upsmp[1] + ups[2] * upsmp[2];
+    sfsp = uns[0] * upsmp[0] + uns[1] * upsmp[1] + uns[2] * upsmp[2];
+}
+
+/*! \brief C++ porting of WMAMP
+ *
+ * \param iis: `int`
+ * \param cost: `double`
+ * \param sint: `double`
+ * \param cosp: `double`
+ * \param sinp: `double`
+ * \param inpol: `int`
+ * \param lm: `int`
+ * \param idot: `int`
+ * \param nsph: `int`
+ * \param arg: `double *`
+ * \param u: `double *`
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param c1: `C1 *`
+ */
+void wmamp(
+		int iis, double cost, double sint, double cosp, double sinp, int inpol,
+		int lm, int idot, int nsph, double *arg, double *u, double *up,
+		double *un, C1 *c1
+) {
+	std::complex<double> *ylm = new std::complex<double>[1682];
+	int nlmp = lm * (lm + 2) + 2;
+	ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
+	if (idot != 0) {
+		if (idot != 1) {
+			for (int n40 = 0; n40 < nsph; n40++) {
+				arg[n40] = u[0] * c1->rxx[n40] + u[1] * c1->ryy[n40] + u[2] * c1->rzz[n40];
+			}
+		} else {
+			for (int n50 = 0; n50 < nsph; n50++) {
+				arg[n50] = c1->rzz[n50];
+			}
+		}
+		if (iis == 2) {
+			for (int n60 = 0; n60 < nsph; n60++) arg[n60] *= -1;
+		}
+	}
+	sphar(cost, sint, cosp, sinp, lm, ylm);
+	//printf("DEBUG: in WMAMP and calling PWMA with lm = %d and iis = %d\n", lm, iis);
+	pwma(up, un, ylm, inpol, lm, iis, c1);
+	delete[] ylm;
+}
+
+/*! \brief C++ porting of WMASP
+ *
+ * \param cost: `double`
+ * \param sint: `double`
+ * \param cosp: `double`
+ * \param sinp: `double`
+ * \param costs: `double`
+ * \param sints: `double`
+ * \param cosps: `double`
+ * \param sinps: `double`
+ * \param u: `double *`
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param us: `double *`
+ * \param ups: `double *`
+ * \param uns: `double *`
+ * \param isq: `int`
+ * \param ibf: `int`
+ * \param inpol: `int`
+ * \param lm: `int`
+ * \param idot: `int`
+ * \param nsph: `int`
+ * \param argi: `double *`
+ * \param args: `double *`
+ * \param c1: `C1 *`
+ */
+void wmasp(
+		double cost, double sint, double cosp, double sinp, double costs, double sints,
+		double cosps, double sinps, double *u, double *up, double *un, double *us,
+		double *ups, double *uns, int isq, int ibf, int inpol, int lm, int idot,
+		int nsph, double *argi, double *args, C1 *c1
+) {
+	std::complex<double> *ylm = new std::complex<double>[1682];
+	int nlmp = lm * (lm + 2) + 2;
+	ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
+	if (idot != 0) {
+		if (idot != 1) {
+			for (int n40 = 1; n40 <= nsph; n40++) {
+				argi[n40 - 1] = u[0] * c1->rxx[n40 - 1] + u[1] * c1->ryy[n40 - 1] + u[2] * c1->rzz[n40 - 1];
+				if (ibf != 0) {
+					args[n40 - 1] = argi[n40 - 1] * ibf;
+				} else {
+					args[n40 - 1] = -1.0 * (us[0] * c1->rxx[n40 - 1] + us[1] * c1->ryy[n40 - 1] + us[2] * c1->rzz[n40 - 1]);
+				}
+			}
+		} else { // label 50
+			for (int n60 = 1; n60 <= nsph; n60++) {
+				argi[n60 - 1] = cost * c1->rzz[n60 - 1];
+				if (ibf != 0) {
+					args[n60 - 1] = argi[n60 - 1] * ibf;
+				} else {
+					args[n60 - 1] = -costs * c1->rzz[n60 - 1];
+				}
+			}
+		}
+	}
+	sphar(cost, sint, cosp, sinp, lm, ylm);
+	//printf("DEBUG: in WMASP and calling PWMA with isq = %d\n", isq);
+	pwma(up, un, ylm, inpol, lm, isq, c1);
+	if (ibf >= 0) {
+		sphar(costs, sints, cosps, sinps, lm, ylm);
+		//printf("DEBUG: in WMASP and calling PWMA with isq = 2 and ibf = %d\n", ibf);
+		pwma(ups, uns, ylm, inpol, lm, 2, c1);
+	}
+	delete[] ylm;
+}
+
+
+/*! \brief C++ porting of DME
+ *
+ * \param li: `int`
+ * \param i: `int`
+ * \param npnt: `int`
+ * \param npntts: `int`
+ * \param vk: `double`
+ * \param exdc: `double`
+ * \param exri: `double`
+ * \param c1: `C1 *`
+ * \param c2: `C2 *`
+ * \param jer: `int &`
+ * \param lcalc: `int &`
+ * \param arg: `complex<double> &`.
+ */
+void dme(
+		int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
+		C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg) {
+	double *rfj = new double[42];
+	double *rfn = new double[42];
+	std::complex<double> cfj[42], fbi[42], fb[42], fn[42];
+	std::complex<double> rmf[40], drmf[40], ref[40], dref[40];
+	std::complex<double> dfbi, dfb, dfn, ccna, ccnb, ccnc, ccnd;
+	std::complex<double> y1, dy1, y2, dy2, arin, cri, uim;
+	jer = 0;
+	uim = std::complex<double>(0.0, 1.0);
+	int nstp = npnt - 1;
+	int nstpts = npntts - 1;
+	int lipo = li + 1;
+	int lipt = li + 2;
+	double sz = vk * c1->ros[i - 1];
+	c2->vsz[i - 1] = sz;
+	double vkr1 = vk * c1->rc[i - 1][0];
+	int nsh = c1->nshl[i - 1];
+	c2->vkt[i - 1] = std::sqrt(c2->dc0[0]);
+	arg = vkr1 * c2->vkt[i - 1];
+	arin = arg;
+	bool goto32 = false;
+	if (arg.imag() != 0.0) {
+		cbf(lipo, arg, lcalc, cfj);
+		if (lcalc < lipo) {
+			jer = 5;
+			delete[] rfj;
+			delete[] rfn;
+			return;
+		}
+		for (int j24 = 1; j24 <= lipt; j24++) fbi[j24 - 1] = cfj[j24 - 1];
+		goto32 = true;
+	}
+	if (!goto32) {
+		rbf(lipo, arg.real(), lcalc, rfj);
+		if (lcalc < lipo) {
+			jer = 5;
+			delete[] rfj;
+			delete[] rfn;
+			return;
+		}
+		for (int j30 = 1; j30 <= lipt; j30++) fbi[j30 - 1] = rfj[j30 - 1];
+	}
+	double arex = sz * exri;
+	arg = arex;
+	rbf(lipo, arex, lcalc, rfj);
+	if (lcalc < lipo) {
+		jer = 7;
+		delete[] rfj;
+		delete[] rfn;
+		return;
+	}
+	rnf(lipo, arex, lcalc, rfn);
+	if (lcalc < lipo) {
+		jer = 8;
+		delete[] rfj;
+		delete[] rfn;
+		return;
+	}
+	for (int j43 = 1; j43 <= lipt; j43++) {
+		fb[j43 - 1] = rfj[j43 - 1];
+		//printf("DEBUG: fb[%d] = (%lE,%lE)\n", j43, fb[j43 - 1].real(), fb[j43 - 1].imag());
+		fn[j43 - 1] = rfn[j43 - 1];
+		//printf("DEBUG: fn[%d] = (%lE,%lE)\n", j43, fn[j43 - 1].real(), fn[j43 - 1].imag());
+	}
+	if (nsh <= 1) {
+		cri = c2->dc0[0] / exdc;
+		for (int l60 = 1; l60 <= li; l60++) {
+			int lpo = l60 + 1;
+			int ltpo = lpo + l60;
+			int lpt = lpo + 1;
+			dfbi = ((1.0 * l60) * fbi[l60 - 1] - (1.0 * lpo) * fbi[lpt - 1]) * arin + fbi[lpo - 1] * (1.0 * ltpo);
+			dfb = ((1.0 * l60) * fb[l60 - 1] - (1.0 * lpo) * fb[lpt - 1]) * arex + fb[lpo - 1] * (1.0 * ltpo);
+			dfn = ((1.0 * l60) * fn[l60 - 1] - (1.0 * lpo) * fn[lpt - 1]) * arex + fn[lpo - 1] * (1.0 * ltpo);
+			ccna = fbi[lpo - 1] * dfn;
+			ccnb = fn[lpo - 1] * dfbi;
+			ccnc = fbi[lpo - 1] * dfb;
+			ccnd = fb[lpo - 1] * dfbi;
+			c1->rmi[l60 - 1][i - 1] = 1.0 + uim * (ccna - ccnb) / (ccnc - ccnd);
+			c1->rei[l60 - 1][i - 1] = 1.0 + uim * (cri * ccna - ccnb) / (cri * ccnc - ccnd);
+		}
+	} else { // nsh > 1
+		int ic = 1;
+		for (int l80 = 1; l80 <= li; l80++) {
+			int lpo = l80 + 1;
+			int ltpo = lpo + l80;
+			int lpt = lpo + 1;
+			int dltpo = ltpo;
+			y1 = fbi[lpo - 1];
+			dy1 = ((1.0 * l80) * fbi[l80 - 1] - (1.0 * lpo) * fbi[lpt - 1]) * c2->vkt[i - 1] / (1.0 * dltpo);
+			y2 = y1;
+			dy2 = dy1;
+			ic = 1;
+			for (int ns76 = 2; ns76 <= nsh; ns76++) {
+				int nsmo = ns76 - 1;
+				double vkr = vk * c1->rc[i - 1][nsmo - 1];
+				if (ns76 % 2 != 0) {
+					ic += 1;
+					double step = 1.0 * nstp;
+					step = vk * (c1->rc[i - 1][ns76 - 1] - c1->rc[i - 1][nsmo - 1]) / step;
+					arg = c2->dc0[ic - 1];
+					rkc(nstp, step, arg, vkr, lpo, y1, y2, dy1, dy2);
+				} else {
+					diel(nstpts, nsmo, i, ic, vk, c1, c2);
+					double stepts = 1.0 * nstpts;
+					stepts = vk * (c1->rc[i - 1][ns76 - 1] - c1->rc[i - 1][nsmo - 1]) / stepts;
+					rkt(nstpts, stepts, vkr, lpo, y1, y2, dy1, dy2, c2);
+				}
+			}
+			rmf[l80 - 1] = y1 * sz;
+			drmf[l80 - 1] = dy1 * sz + y1;
+			ref[l80 - 1] = y2 * sz;
+			dref[l80 - 1] = dy2 * sz + y2;
+		}
+		cri = 1.0 + uim * 0.0;
+		if (nsh % 2 != 0) cri = c2->dc0[ic - 1] / exdc;
+		for (int l90 = 1; l90 <= li; l90++) {
+			int lpo = l90 + 1;
+			int ltpo = lpo + l90;
+			int lpt = lpo + 1;
+			dfb = ((1.0 * l90) * fb[l90 - 1] - (1.0 * lpo) * fb[lpt - 1]) * arex + fb[lpo - 1] * (1.0 * ltpo);
+			dfn = ((1.0 * l90) * fn[l90 - 1] - (1.0 * lpo) * fn[lpt - 1]) * arex + fn[lpo - 1] * (1.0 * ltpo);
+			ccna = rmf[l90 - 1] * dfn;
+			ccnb = drmf[l90 - 1] * fn[lpo - 1] * (1.0 * sz * ltpo);
+			ccnc = rmf[l90 - 1] * dfb;
+			ccnd = drmf[l90 - 1] * fb[lpo -1] * (1.0 * sz * ltpo);
+			c1->rmi[l90 - 1][i - 1] = 1.0 + uim *(ccna - ccnb) / (ccnc - ccnd);
+			//printf("DEBUG: gone 90, rmi[%d][%d] = (%lE,%lE)\n", l90, i, c1->rmi[l90 - 1][i - 1].real(), c1->rmi[l90 - 1][i - 1].imag());
+			ccna = ref[l90 - 1] * dfn;
+			ccnb = dref[l90 - 1] * fn[lpo - 1] * (1.0 * sz * ltpo);
+			ccnc = ref[l90 - 1] * dfb;
+			ccnd = dref[l90 - 1] *fb[lpo - 1] * (1.0 * sz * ltpo);
+			c1->rei[l90 - 1][i - 1] = 1.0 + uim * (cri * ccna - ccnb) / (cri * ccnc - ccnd);
+			//printf("DEBUG: gone 90, rei[%d][%d] = (%lE,%lE)\n", l90, i, c1->rei[l90 - 1][i - 1].real(), c1->rei[l90 - 1][i - 1].imag());
+		}
+	} // nsh <= 1 ?
+	delete[] rfj;
+	delete[] rfn;
+	return;
+}
+
+#endif /* SRC_INCLUDE_SPH_SUBS_H_ */

src/libnptm/Commons.cpp

+/*! \file Commons.cpp
+ *
+ */
+
+#include "../include/Commons.h"
+
+using namespace std;
+
+C1::C1(int ns, int l_max) {
+	nlmmt = 2 * (l_max * (l_max + 2));
+	nsph = ns;
+	lm = l_max;
+
+	rmi = new complex<double>*[lm];
+	rei = new complex<double>*[lm];
+	for (int ri = 0; ri < lm; ri++) {
+		rmi[ri] = new complex<double>[nsph];
+		rei[ri] = new complex<double>[nsph];
+	}
+	w = new complex<double>*[nlmmt];
+	for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4];
+	vints = new complex<double>*[nsph];
+	rc = new double*[nsph];
+	for (int vi = 0; vi < nsph; vi++) {
+		rc[vi] = new double[lm];
+		vints[vi] = new complex<double>[16];
+	}
+	fsas = new complex<double>[nsph];
+	sscs = new double[nsph];
+	sexs = new double[nsph];
+	sabs = new double[nsph];
+	sqscs = new double[nsph];
+	sqexs = new double[nsph];
+	sqabs = new double[nsph];
+	gcsv = new double[nsph];;
+	rxx = new double[nsph];
+	ryy = new double[nsph];
+	rzz = new double[nsph];
+	ros = new double[nsph];
+	iog = new int[nsph];
+	nshl = new int[nsph];
+
+	sas = new complex<double>**[nsph];
+	for (int si = 0; si < nsph; si++) {
+		sas[si] = new complex<double>*[2];
+		sas[si][0] = new complex<double>[2];
+		sas[si][1] = new complex<double>[2];
+	}
+}
+
+C1::~C1() {
+	delete[] rmi;
+	delete[] rei;
+	for (int wi = 1; wi <= nlmmt; wi++) delete[] w[wi];
+	for (int vi = 1; vi <= nsph; vi++) {
+		delete[] rc[nsph - vi];
+		delete[] vints[nsph - vi];
+	}
+	for (int si = 1; si <= nsph; si++) {
+		delete[] sas[nsph - si][1];
+		delete[] sas[nsph - si][0];
+	}
+	delete[] fsas;
+	delete[] sscs;
+	delete[] sexs;
+	delete[] sabs;
+	delete[] sqscs;
+	delete[] sqexs;
+	delete[] sqabs;
+	delete[] gcsv;
+	delete[] rxx;
+	delete[] ryy;
+	delete[] rzz;
+	delete[] ros;
+	delete[] iog;
+	delete[] nshl;
+}
+
+C2::C2(int ns, int nl, int npnt, int npntts) {
+	nsph = ns;
+	int max_n = (npnt > npntts) ? npnt : npntts;
+	nhspo = 2 * max_n - 1;
+	ris = new complex<double>[nhspo];
+	dlri = new complex<double>[nhspo];
+	vkt = new complex<double>[nsph];
+	dc0 = new complex<double>[nl];
+	vsz = new double[nsph];
+}
+
+C2::~C2() {
+	delete[] ris;
+	delete[] dlri;
+	delete[] vkt;
+	delete[] dc0;
+	delete[] vsz;
+}

src/libnptm/Configuration.cpp

+/*! \file Configuration.cpp
+*/
+
+#include <cmath>
+#include <cstdio>
+#include <fstream>
+#include <string>
+#include "../include/List.h"
+#include "../include/Parsers.h"
+#include "../include/Configuration.h"
+
+using namespace std;
+
+GeometryConfiguration::GeometryConfiguration(
+		int nsph, int lm, int _in_pol, int _npnt, int _npntts, int isam,
+		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 in_ph_start, double in_ph_step, double in_ph_end,
+		double sc_ph_start, double sc_ph_step, double sc_ph_end,
+		int _jwtm
+	) {
+	number_of_spheres = nsph;
+	l_max = lm;
+	in_pol = _in_pol;
+	npnt = _npnt;
+	npntts = _npntts;
+	meridional_type = isam;
+	in_theta_start = in_th_start;
+	in_theta_step = in_th_step;
+	in_theta_end = in_th_end;
+	in_phi_start = in_ph_start;
+	in_phi_step = in_ph_step;
+	in_phi_end = in_ph_end;
+	sc_theta_start = sc_th_start;
+	sc_theta_step = sc_th_step;
+	sc_theta_end = sc_th_end;
+	sc_phi_start = sc_ph_start;
+	sc_phi_step = sc_ph_step;
+	sc_phi_end = sc_ph_end;
+	jwtm = _jwtm;
+	sph_x = x;
+	sph_y = y;
+	sph_z = z;
+}
+
+GeometryConfiguration::~GeometryConfiguration() {
+	delete[] sph_x;
+	delete[] sph_y;
+	delete[] sph_z;
+}
+
+GeometryConfiguration* GeometryConfiguration::from_legacy(string file_name) {
+	int num_lines = 0;
+	int last_read_line = 0;
+	string *file_lines;
+	try {
+		file_lines = load_file(file_name, &num_lines);
+	} catch (exception &ex) {
+		throw OpenConfigurationFileException(file_name);
+	}
+	int nsph, lm, _in_pol, _npnt, _npntts, isam;
+	sscanf(
+			file_lines[last_read_line++].c_str(),
+			" %d %d %d %d %d %d",
+			&nsph, &lm, &_in_pol, &_npnt, &_npntts, &isam
+	);
+	double *x, *y, *z;
+	x = new double[nsph];
+	y = new double[nsph];
+	z = new double[nsph];
+	if (nsph == 1) {
+		x[0] = 0.0;
+		y[0] = 0.0;
+		z[0] = 0.0;
+	} else {
+		for (int i = 0; i < nsph; i++) {
+			double sph_x, sph_y, sph_z;
+			int sph_x_exp, sph_y_exp, sph_z_exp;
+			sscanf(
+					file_lines[last_read_line++].c_str(),
+					" %lf D%d %lf D%d %lf D%d",
+					&sph_x, &sph_x_exp, &sph_y, &sph_y_exp, &sph_z, &sph_z_exp
+			);
+			x[i] = sph_x * pow(10.0, 1.0 * sph_x_exp);
+			y[i] = sph_y * pow(10.0, 1.0 * sph_y_exp);
+			z[i] = sph_z * pow(10.0, 1.0 * sph_z_exp);
+		}
+	}
+	double in_th_start, in_th_end, in_th_step, sc_th_start, sc_th_end, sc_th_step;
+	int in_th_start_exp, in_th_end_exp, in_th_step_exp, sc_th_start_exp, sc_th_end_exp, sc_th_step_exp;
+	sscanf(
+			file_lines[last_read_line++].c_str(),
+			" %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d",
+			&in_th_start, &in_th_start_exp,
+			&in_th_step, &in_th_step_exp,
+			&in_th_end, &in_th_end_exp,
+			&sc_th_start, &sc_th_start_exp,
+			&sc_th_step, &sc_th_step_exp,
+			&sc_th_end, &sc_th_end_exp
+	);
+	in_th_start *= pow(10.0, 1.0 * in_th_start_exp);
+	in_th_step *= pow(10.0, 1.0 * in_th_step_exp);
+	in_th_end *= pow(10.0, 1.0 * in_th_end_exp);
+	sc_th_start *= pow(10.0, 1.0 * sc_th_start_exp);
+	sc_th_step *= pow(10.0, 1.0 * sc_th_step_exp);
+	sc_th_end *= pow(10.0, 1.0 * sc_th_end_exp);
+	double in_ph_start, in_ph_end, in_ph_step, sc_ph_start, sc_ph_end, sc_ph_step;
+	int in_ph_start_exp, in_ph_end_exp, in_ph_step_exp, sc_ph_start_exp, sc_ph_end_exp, sc_ph_step_exp;
+	sscanf(
+			file_lines[last_read_line++].c_str(),
+			" %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d %lf D%d",
+			&in_ph_start, &in_ph_start_exp,
+			&in_ph_step, &in_ph_step_exp,
+			&in_ph_end, &in_ph_end_exp,
+			&sc_ph_start, &sc_ph_start_exp,
+			&sc_ph_step, &sc_ph_step_exp,
+			&sc_ph_end, &sc_ph_end_exp
+	);
+	in_ph_start *= pow(10.0, 1.0 * in_ph_start_exp);
+	in_ph_step *= pow(10.0, 1.0 * in_ph_step_exp);
+	in_ph_end *= pow(10.0, 1.0 * in_ph_end_exp);
+	sc_ph_start *= pow(10.0, 1.0 * sc_ph_start_exp);
+	sc_ph_step *= pow(10.0, 1.0 * sc_ph_step_exp);
+	sc_ph_end *= pow(10.0, 1.0 * sc_ph_end_exp);
+	int _jwtm;
+	sscanf(file_lines[last_read_line++].c_str(), " %d", &_jwtm);
+	GeometryConfiguration *conf = new GeometryConfiguration(
+			nsph, lm, _in_pol, _npnt, _npntts, isam,
+			x, y, z,
+			in_th_start, in_th_step, in_th_end,
+			sc_th_start, sc_th_step, sc_th_end,
+			in_ph_start, in_ph_step, in_ph_end,
+			sc_ph_start, sc_ph_step, sc_ph_end,
+			_jwtm
+	);
+	return conf;
+}
+
+ScattererConfiguration::ScattererConfiguration(
+		int nsph,
+		double *scale_vector,
+		int nxi,
+		string variable_name,
+		int *iog_vector,
+		double *ros_vector,
+		int *nshl_vector,
+		double **rcf_vector,
+		int dielectric_func_type,
+		complex<double> ***dc_matrix,
+		bool is_external,
+		double ex,
+		double w,
+		double x
+		) {
+	number_of_spheres = nsph;
+	scale_vec = scale_vector;
+	number_of_scales = nxi;
+	reference_variable_name = variable_name;
+	iog_vec = iog_vector;
+	radii_of_spheres = ros_vector;
+	nshl_vec = nshl_vector;
+	rcf = rcf_vector;
+	idfc = dielectric_func_type,
+	dc0_matrix = dc_matrix;
+	use_external_sphere = is_external;
+	exdc = ex;
+	wp = w;
+	xip = x;
+}
+
+ScattererConfiguration::~ScattererConfiguration() {
+	int max_ici = 0;
+	for (int i = 1; i <= number_of_spheres; i++) {
+		if (iog_vec[i - 1] < i) continue;
+		int nsh = nshl_vec[i - 1];
+		if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
+	}
+	for (int i = 0; i < max_ici; i++) {
+		for (int j = 0; j < number_of_spheres; j++) {
+			delete[] dc0_matrix[i][j];
+		}
+	}
+	for (int i = 0; i < number_of_spheres; i++) {
+		delete[] rcf[i];
+	}
+	delete[] nshl_vec;
+	delete[] radii_of_spheres;
+	delete[] iog_vec;
+	delete[] scale_vec;
+}
+
+ScattererConfiguration* ScattererConfiguration::from_binary(string file_name, string mode) {
+	int nsph;
+	int *iog;
+	double _exdc, _wp, _xip;
+	int _idfc, nxi;
+	int *nshl_vector;
+	double *xi_vec;
+	double *ros_vector;
+	double **rcf_vector;
+	complex<double> ***dc0m;
+	if (mode.compare("LEGACY") == 0) { // Legacy mode was chosen.
+		int max_ici = 0;
+		fstream input;
+		input.open(file_name.c_str(), ios::in | ios::binary);
+		if (input.is_open()) {
+			input.read(reinterpret_cast<char *>(&nsph), sizeof(int));
+			iog = new int[nsph];
+			for (int i = 0; i < nsph; i++) {
+				input.read(reinterpret_cast<char *>(&(iog[i])), sizeof(int));
+			}
+			input.read(reinterpret_cast<char *>(&_exdc), sizeof(double));
+			input.read(reinterpret_cast<char *>(&_wp), sizeof(double));
+			input.read(reinterpret_cast<char *>(&_xip), sizeof(double));
+			input.read(reinterpret_cast<char *>(&_idfc), sizeof(int));
+			input.read(reinterpret_cast<char *>(&nxi), sizeof(int));
+			try {
+				xi_vec = new double[nxi];
+			} catch (bad_alloc &ex) {
+				throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of scales " + nxi);
+			}
+			for (int i = 0; i < nxi; i++) {
+				input.read(reinterpret_cast<char *>(&(xi_vec[i])), sizeof(double));
+			}
+			nshl_vector = new int[nsph];
+			ros_vector = new double[nsph];
+			rcf_vector = new double*[nsph];
+			for (int i115 = 1; i115 <= nsph; i115++) {
+				if (iog[i115 - 1] < i115) continue;
+				input.read(reinterpret_cast<char *>(&(nshl_vector[i115 - 1])), sizeof(int));
+				input.read(reinterpret_cast<char *>(&(ros_vector[i115 - 1])), sizeof(double));
+				int nsh = nshl_vector[i115 - 1];
+				if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
+				try {
+					rcf_vector[i115 - 1] = new double[nsh];
+				} catch (bad_alloc &ex) {
+					throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of layers " + nsh);
+				}
+				for (int nsi = 0; nsi < nsh; nsi++) {
+					input.read(reinterpret_cast<char *>(&(rcf_vector[i115 - 1][nsi])), sizeof(double));
+				}
+			}
+			dc0m = new complex<double>**[max_ici];
+			for (int dim1 = 0; dim1 < max_ici; dim1++) {
+				dc0m[dim1] = new complex<double>*[nsph];
+			    for (int dim2 = 0; dim2 < nsph; dim2++) {
+			    	dc0m[dim1][dim2] = new complex<double>[nxi];
+			    }
+			}
+			for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
+				if (_idfc != 0 && jxi468 > 1) continue;
+			    for (int i162 = 1; i162 <= nsph; i162++) {
+			    	if (iog[i162 - 1] < i162) continue;
+			    	int nsh = nshl_vector[i162 - 1];
+			    	int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
+			    	for (int i157 = 0; i157 < ici; i157++) {
+			    		double dc0_real, dc0_img;
+			    		input.read(reinterpret_cast<char *>(&dc0_real), sizeof(double));
+			    		input.read(reinterpret_cast<char *>(&dc0_img), sizeof(double));
+			    		dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;
+			    	}
+			    }
+			}
+			input.close();
+		} else { // Opening of the input file did not succeed.
+			throw OpenConfigurationFileException(file_name);
+		}
+	} else { // A different binary format was chosen.
+		//TODO: this part is not yet implemented.
+		//      Functions to write optimized file formats may be invoked here.
+	}
+	ScattererConfiguration *conf = new ScattererConfiguration(
+			nsph,
+			xi_vec,
+			nxi,
+			"XIV",
+			iog,
+			ros_vector,
+			nshl_vector,
+			rcf_vector,
+			_idfc,
+			dc0m,
+			false,
+			_exdc,
+			_wp,
+			_xip
+	);
+	return conf;
+}
+
+ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_name) {
+	int num_lines = 0;
+	int last_read_line = 0;
+	string *file_lines;
+	try {
+		file_lines = load_file(dedfb_file_name, &num_lines);
+	} catch (exception &ex) {
+		throw OpenConfigurationFileException(dedfb_file_name);
+	}
+	int nsph, ies;
+	int max_ici = 0;
+	sscanf(file_lines[last_read_line].c_str(), " %d %d", &nsph, &ies);
+	if (ies != 0) ies = 1;
+	double _exdc, _wp, _xip;
+	int exdc_exp, wp_exp, xip_exp;
+	int _idfc, nxi, instpc, insn;
+	sscanf(
+			file_lines[++last_read_line].c_str(),
+			" %lf D%d %lf D%d %lf D%d %d %d %d %d",
+			&_exdc, &exdc_exp,
+			&_wp, &wp_exp,
+			&_xip, &xip_exp,
+			&_idfc, &nxi, &instpc, &insn
+	);
+	_exdc *= pow(10.0, 1.0 * 1.0 * exdc_exp);
+	_wp *= pow(10.0, 1.0 * wp_exp);
+	_xip *= pow(10.0, 1.0 * xip_exp);
+
+	double *variable_vector;
+	string variable_name;
+
+	if (_idfc < 0) { // Diel. functions at XIP value and XI is scale factor
+		variable_name = "XIV";
+		if (instpc < 1) { // The variable vector is explicitly defined.
+			double xi;
+			int xi_exp;
+			List<double> xi_vector;
+			sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
+			xi *= pow(10.0, 1.0 * xi_exp);
+			xi_vector.set(0, xi);
+			for (int jxi310 = 1; jxi310 < nxi; jxi310++) {
+				sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d", &xi, &xi_exp);
+				xi *= pow(10.0, 1.0 * xi_exp);
+				xi_vector.append(xi);
+			}
+			variable_vector = xi_vector.to_array();
+		} else { // instpc >= 1: the variable vector is defined in steps
+			double xi, xi_step;
+			int xi_exp, xi_step_exp;
+			sscanf(file_lines[++last_read_line].c_str(), " %9lE D%d %9lE D%d", &xi, &xi_exp, &xi_step, &xi_step_exp);
+			xi *= pow(10.0, 1.0 * xi_exp);
+			xi_step *= pow(10.0, 1.0 * xi_step_exp);
+			variable_vector = new double[nxi];
+			for (int jxi320 = 0; jxi320 < nxi; jxi320++) {
+				variable_vector[jxi320] = xi;
+				xi += xi_step;
+			}
+		}
+	} else { // idfc >= 0
+		variable_vector = new double[nxi];
+		if (instpc == 0) { // The variable vector is explicitly defined
+			double vs;
+			int vs_exp;
+			for (int jxi_r = 0; jxi_r < nxi; jxi_r++) {
+				sscanf(file_lines[++last_read_line].c_str(), " %lf D%d", &vs, &vs_exp);
+				vs *= pow(10.0, 1.0 * vs_exp);
+				variable_vector[jxi_r] = vs;
+			}
+			switch (insn) {
+			case 1: //xi vector definition
+				variable_name = "XIV";
+				break;
+			case 2: //wave number vector definition
+				variable_name = "WNS";
+		    	break;
+			case 3: //wavelength vector definition
+				variable_name = "WLS";
+				break;
+			case 4: //pu vector definition
+				variable_name = "PUS";
+				break;
+			case 5: //eV vector definition
+				variable_name = "EVS";
+				break;
+		    }
+		} else { // The variable vector needs to be computed in steps
+			double vs, vs_step;
+			int vs_exp, vs_step_exp;
+			sscanf(file_lines[++last_read_line].c_str(), " %lf D%d %lf D%d", &vs, &vs_exp, &vs_step, &vs_step_exp);
+			vs *= pow(10.0, 1.0 * vs_exp);
+			vs_step *= pow(10.0, 1.0 * vs_step_exp);
+			for (int jxi110w = 0; jxi110w < nxi; jxi110w++) {
+				variable_vector[jxi110w] = vs;
+				vs += vs_step;
+			}
+			switch (insn) {
+			case 1: //xi vector definition
+				variable_name = "XIV";
+				break;
+			case 2: //wave number vector definition
+				variable_name = "WNS";
+				break;
+			case 3: //wavelength vector definition
+				variable_name = "WLS";
+				break;
+			case 4: //pu vector definition
+				variable_name = "PUS";
+				break;
+			case 5: //eV vector definition
+				variable_name = "EVS";
+				break;
+			}
+		}
+	}
+	last_read_line++;
+	int *iog_vector = new int[nsph];
+	double *ros_vector = new double[nsph];
+	double **rcf_vector = new double*[nsph];
+	int *nshl_vector = new int[nsph];
+	for (int i = 0; i < nsph; i++) {
+		string read_format = "";
+		for (int j = 0; j < i; j++) read_format += " %*d";
+		read_format += " %d";
+		sscanf(file_lines[last_read_line].c_str(), read_format.c_str(), (iog_vector + i));
+	}
+	for (int i113 = 1; i113 <= nsph; i113++) {
+		int i_val, nsh;
+	    double ros_val;
+	    int ros_val_exp;
+	    if (iog_vector[i113 - 1] < i113) continue;
+	    sscanf(file_lines[++last_read_line].c_str(), " %d %lf D%d", &i_val, &ros_val, &ros_val_exp);
+	    nshl_vector[i113 - 1] = i_val;
+	    if (max_ici < (i_val + 1) / 2) max_ici = (i_val + 1) / 2;
+	    ros_vector[i113 - 1] = ros_val * pow(10.0, 1.0 * ros_val_exp);
+	    nsh = nshl_vector[i113 - 1];
+	    if (i113 == 1) nsh += ies;
+	    rcf_vector[i113 - 1] = new double[nsh];
+	    for (int ns = 0; ns < nsh; ns++) {
+	    	double ns_rcf;
+	    	int ns_rcf_exp;
+	    	sscanf(file_lines[++last_read_line].c_str(), " %lf D%d", &ns_rcf, &ns_rcf_exp);
+	    	rcf_vector[i113 -1][ns] = ns_rcf * pow(10.0, 1.0 * ns_rcf_exp);
+	    }
+	}
+	complex<double> ***dc0m = new complex<double>**[max_ici];
+	for (int dim1 = 0; dim1 < max_ici; dim1++) {
+		dc0m[dim1] = new complex<double>*[nsph];
+	    for (int dim2 = 0; dim2 < nsph; dim2++) {
+	    	dc0m[dim1][dim2] = new complex<double>[nxi];
+	    }
+	}
+	for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
+		if (_idfc != 0 && jxi468 > 1) continue;
+	    for (int i162 = 1; i162 <= nsph; i162++) {
+	    	if (iog_vector[i162 - 1] < i162) continue;
+	    	int nsh = nshl_vector[i162 - 1];
+	    	int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
+	    	if (i162 == 1) ici = ici + ies;
+	    	for (int i157 = 0; i157 < ici; i157++) {
+	    		double dc0_real, dc0_img;
+	    		int dc0_real_exp, dc0_img_exp;
+	    		sscanf(file_lines[++last_read_line].c_str(), " (%lf D%d, %lf D%d)", &dc0_real, &dc0_real_exp, &dc0_img, &dc0_img_exp);
+	    		dc0_real *= pow(10.0, 1.0 * dc0_real_exp);
+	    		dc0_img *= pow(10.0, 1.0 * dc0_img_exp);
+	    		dc0m[i157][i162 - 1][jxi468 - 1] = dc0_real + 1i * dc0_img;
+	    	}
+	    }
+	}
+
+	ScattererConfiguration *config = new ScattererConfiguration(
+			nsph,
+			variable_vector,
+			nxi,
+			variable_name,
+			iog_vector,
+			ros_vector,
+			nshl_vector,
+			rcf_vector,
+			_idfc,
+			dc0m,
+			(ies > 0),
+			_exdc,
+			_wp,
+			_xip
+	);
+	return config;
+}
+
+void ScattererConfiguration::print() {
+	int ies = (use_external_sphere)? 1 : 0;
+	int max_ici = 0;
+	printf("### CONFIGURATION DATA ###\n");
+	printf("NSPH  = %d\n", number_of_spheres);
+	printf("ROS   = [");
+	for (int i = 0; i < number_of_spheres; i++) printf("\t%lg", radii_of_spheres[i]);
+	printf("\t]\n");
+	printf("IOG   = [");
+	for (int i = 0; i < number_of_spheres; i++) printf("\t%d", iog_vec[i]);
+	printf("\t]\n");
+	printf("NSHL  = [");
+	for (int i = 0; i < number_of_spheres; i++) printf("\t%d", nshl_vec[i]);
+	printf("\t]\n");
+	printf("RCF   = [\n");
+	for (int i = 1; i <= number_of_spheres; i++) {
+		int nsh = nshl_vec[i - 1];
+	    if (i == 1) nsh += ies;
+	    if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
+		printf("         [");
+	    for (int ns = 0; ns < nsh; ns++) {
+	    	printf("\t%lg", rcf[i - 1][ns]);
+	    }
+	    printf("\t]\n");
+	}
+	printf("        ]\n");
+	printf("SCALE = %s\n", reference_variable_name.c_str());
+	printf("NXI   = %d\n", number_of_scales);
+	printf("VEC   = [");
+	for (int i = 0; i < number_of_scales; i++) printf("\t%lg", scale_vec[i]);
+	printf("\t]\n");
+	printf("DC0M  = [\n");
+	for (int i = 0; i < max_ici; i++) {
+		printf("         [\n");
+	    for (int j = 0; j < number_of_spheres; j++) {
+			printf("          [");
+			for (int k = 0; k < number_of_scales; k++) {
+				if (idfc != 0 and k > 0) continue;
+		    	printf("\t%lg + i(%lg)", dc0_matrix[i][j][k].real(), dc0_matrix[i][j][k].imag());
+			}
+			printf("\t]\n");
+	    }
+	    printf("         ]\n");
+	}
+	printf("        ]\n");
+}
+
+void ScattererConfiguration::write_binary(string file_name, string mode) {
+	const double two_pi = acos(0.0) * 4.0;
+	const double evc = 6.5821188e-16;
+	int max_ici = 0;
+	if (mode.compare("LEGACY") == 0) { // Legacy mode was chosen.
+		fstream output;
+		int ies = (use_external_sphere)? 1 : 0;
+		double *xi_vec;
+		if (reference_variable_name.compare("XIV") == 0) xi_vec = scale_vec;
+		else {
+			xi_vec = new double[number_of_scales];
+			if (reference_variable_name.compare("WNS") == 0) {
+				for (int i = 0; i < number_of_scales; i++)
+					xi_vec[i] = 3.0e8 * scale_vec[i] / wp;
+			} else if (reference_variable_name.compare("WLS") == 0) {
+				for (int i = 0; i < number_of_scales; i++) {
+					double wn = two_pi / scale_vec[i];
+					xi_vec[i] = 3.0e8 * wn / wp;
+				}
+			} else if (reference_variable_name.compare("PUS") == 0) {
+				for (int i = 0; i < number_of_scales; i++)
+					xi_vec[i] = scale_vec[i] / wp;
+			} else if (reference_variable_name.compare("EVS") == 0) {
+				for (int i = 0; i < number_of_scales; i++) {
+					double pu = scale_vec[i] / evc;
+					xi_vec[i] = pu / wp;
+				}
+			} else {
+				throw UnrecognizedConfigurationException(
+						"Wrong parameter set: unrecognized scale type "
+						+ reference_variable_name
+				);
+			}
+		}
+		output.open(file_name.c_str(), ios::out | ios::binary);
+		output.write(reinterpret_cast<char *>(&number_of_spheres), sizeof(int));
+		for (int i = 0; i < number_of_spheres; i++)
+			output.write(reinterpret_cast<char *>(&(iog_vec[i])), sizeof(int));
+		output.write(reinterpret_cast<char *>(&exdc), sizeof(double));
+		output.write(reinterpret_cast<char *>(&wp), sizeof(double));
+		output.write(reinterpret_cast<char *>(&xip), sizeof(double));
+		output.write(reinterpret_cast<char *>(&idfc), sizeof(int));
+		output.write(reinterpret_cast<char *>(&number_of_scales), sizeof(int));
+		for (int i = 0; i < number_of_scales; i++)
+			output.write(reinterpret_cast<char *>(&(xi_vec[i])), sizeof(double));
+		for (int i115 = 1; i115 <= number_of_spheres; i115++) {
+			if (iog_vec[i115 - 1] < i115) continue;
+			output.write(reinterpret_cast<char *>(&(nshl_vec[i115 - 1])), sizeof(int));
+			output.write(reinterpret_cast<char *>(&(radii_of_spheres[i115 - 1])), sizeof(double));
+			int nsh = nshl_vec[i115 - 1];
+			if (i115 == 1) nsh += ies;
+			if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
+			for (int nsi = 0; nsi < nsh; nsi++)
+				output.write(reinterpret_cast<char *>(&(rcf[i115 - 1][nsi])), sizeof(double));
+		}
+		for (int jxi468 = 1; jxi468 <= number_of_scales; jxi468++) {
+			if (idfc != 0 && jxi468 > 1) continue;
+		    for (int i162 = 1; i162 <= number_of_spheres; i162++) {
+		    	if (iog_vec[i162 - 1] < i162) continue;
+		    	int nsh = nshl_vec[i162 - 1];
+		    	int ici = (nsh + 1) / 2; // QUESTION: is integer division really intended here?
+		    	if (i162 == 1) ici = ici + ies;
+		    	for (int i157 = 0; i157 < ici; i157++) {
+		    		double dc0_real, dc0_img;
+		    		dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
+		    		dc0_img = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
+		    		// The FORTRAN code writes the complex numbers as a 16-byte long binary stream.
+		    		// Here we assume that the 16 bytes are equally split in 8 bytes to represent the
+		    		// real part and 8 bytes to represent the imaginary one.
+		    		output.write(reinterpret_cast<char *>(&dc0_real), sizeof(double));
+		    		output.write(reinterpret_cast<char *>(&dc0_img), sizeof(double));
+		    	}
+		    }
+		}
+		output.close();
+	}
+}
+
+void ScattererConfiguration::write_formatted(string file_name) {
+	const double evc = 6.5821188e-16;
+	const double two_pi = acos(0.0) * 4.0;
+	double *xi_vec;
+	int ici;
+	int ies = (use_external_sphere)? 1: 0;
+	FILE *output = fopen(file_name.c_str(), "w");
+	int scale_type = -1;
+	if (reference_variable_name.compare("XIV") == 0) scale_type = 0;
+	else if (reference_variable_name.compare("WNS") == 0) scale_type = 1;
+	else if (reference_variable_name.compare("WLS") == 0) scale_type = 2;
+	else if (reference_variable_name.compare("PUS") == 0) scale_type = 3;
+	else if (reference_variable_name.compare("EVS") == 0) scale_type = 4;
+	if (idfc >= 0) { // Dielectric functions are constant or they depend on XI
+		double  *pu_vec, *ev_vec, *wn_vec, *wl_vec;
+	    xi_vec = new double[number_of_scales];
+	    pu_vec = new double[number_of_scales];
+	    ev_vec = new double[number_of_scales];
+	    wn_vec = new double[number_of_scales];
+	    wl_vec = new double[number_of_scales];
+		switch (scale_type) {
+		case 0:
+			fprintf(output, "  JXI     XIV          WNS          WLS          PUS          EVS\n");
+			for (int i = 0; i < number_of_scales; i++) {
+				xi_vec[i] = scale_vec[i];
+				pu_vec[i] = xi_vec[i] * wp;
+				ev_vec[i] = pu_vec[i] * evc;
+				wn_vec[i] = pu_vec[i] / 3.0e8;
+				wl_vec[i] = two_pi / wn_vec[i];
+				fprintf(
+						output,
+						"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+						(i + 1),
+						xi_vec[i],
+						wn_vec[i],
+						wl_vec[i],
+						pu_vec[i],
+						ev_vec[i]
+				);
+			}
+			break;
+		case 1:
+			fprintf(output, "  JXI     WNS          WLS          PUS          EVS          XIV\n");
+			for (int i = 0; i < number_of_scales; i++) {
+				wn_vec[i] = scale_vec[i];
+				wl_vec[i] = two_pi / wn_vec[i];
+				xi_vec[i] = 3.0e8 * wn_vec[i] / wp;
+				pu_vec[i] = xi_vec[i] * wp;
+				ev_vec[i] = pu_vec[i] * evc;
+				fprintf(
+						output,
+						"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+						(i + 1),
+						wn_vec[i],
+						wl_vec[i],
+						pu_vec[i],
+						ev_vec[i],
+						xi_vec[i]
+				);
+			}
+			break;
+		case 2:
+			fprintf(output, "  JXI     WLS          WNS          PUS          EVS          XIV\n");
+			for (int i = 0; i < number_of_scales; i++) {
+				wl_vec[i] = scale_vec[i];
+				wn_vec[i] = two_pi / wl_vec[i];
+				xi_vec[i] = 3.0e8 * wn_vec[i] / wp;
+				pu_vec[i] = xi_vec[i] * wp;
+				ev_vec[i] = pu_vec[i] * evc;
+				fprintf(
+						output,
+						"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+						(i + 1),
+						wl_vec[i],
+						wn_vec[i],
+						pu_vec[i],
+						ev_vec[i],
+						xi_vec[i]
+				);
+			}
+			break;
+		case 3:
+			fprintf(output, "  JXI     PUS          WNS          WLS          EVS          XIV\n");
+			for (int i = 0; i < number_of_scales; i++) {
+				pu_vec[i] = scale_vec[i];
+				xi_vec[i] = pu_vec[i] / wp;
+				wn_vec[i] = pu_vec[i] / 3.0e8;
+				wl_vec[i] = two_pi / wn_vec[i];
+				ev_vec[i] = pu_vec[i] * evc;
+				fprintf(
+						output,
+						"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+						(i + 1),
+						pu_vec[i],
+						wn_vec[i],
+						wl_vec[i],
+						ev_vec[i],
+						xi_vec[i]
+				);
+			}
+			break;
+		case 4:
+			fprintf(output, "  JXI     EVS          WNS          WLS          PUS          XIV\n");
+			for (int i = 0; i < number_of_scales; i++) {
+				ev_vec[i] = scale_vec[i];
+				pu_vec[i] = ev_vec[i] / evc;
+				xi_vec[i] = pu_vec[i] / wp;
+				wn_vec[i] = pu_vec[i] / 3.0e8;
+				wl_vec[i] = two_pi / wn_vec[i];
+				fprintf(
+						output,
+						"%5d %13.4lE %13.4lE %13.4lE %13.4lE %13.4lE\n",
+						(i + 1),
+						ev_vec[i],
+						wn_vec[i],
+						wl_vec[i],
+						pu_vec[i],
+						xi_vec[i]
+				);
+			}
+			break;
+		default:
+			throw UnrecognizedConfigurationException(
+					"Wonrg parameter set: unrecognized scale variable type "
+					+ reference_variable_name
+			);
+			break;
+		}
+	} else { // idfc < 0, Dielectric functions are at XIP and XI is scale for dimensions
+		double pu, wn;
+		xi_vec = scale_vec;
+		pu = xip * wp;
+		wn = pu / 3.0e8;
+		fprintf(output, "          XIP          WN           WL           PU           EV\n");
+		fprintf(output, "     %13.4lE", xip);
+		fprintf(output, "%13.4lE", wn);
+		fprintf(output, "%13.4lE", two_pi / wn);
+		fprintf(output, "%13.4lE", pu);
+		fprintf(output, "%13.4lE\n", pu * evc);
+		fprintf(output, "  SCALE FACTORS XI\n");
+		for (int i = 0; i < number_of_scales; i++)
+			fprintf(output, "%5d%13.4lE\n", (i + 1), xi_vec[i]);
+	}
+	if (idfc != 0) {
+		fprintf(output, "  DIELECTRIC CONSTANTS\n");
+	    for (int i473 = 1; i473 <= number_of_spheres; i473++) {
+	    	if (iog_vec[i473 - 1] != i473) continue;
+	    	ici = (nshl_vec[i473 - 1] + 1) / 2;
+	    	if (i473 == 1) ici += ies;
+	    	fprintf(output, " SPHERE N. %4d\n", i473);
+	    	for (int ic472 = 0; ic472 < ici; ic472++) {
+	    		double dc0_real = dc0_matrix[ic472][i473 - 1][0].real(), dc0_img = dc0_matrix[ic472][i473 - 1][0].imag();
+	    		fprintf(output, "%5d %12.4lE%12.4lE\n", (ic472 + 1), dc0_real, dc0_img);
+	    	}
+	    }
+	  } else {
+		  fprintf(output, "  DIELECTRIC FUNCTIONS\n");
+		  for (int i478 = 1; i478 <= number_of_spheres; i478++) {
+			  if (iog_vec[i478 - 1] != i478) continue;
+			  ici = (nshl_vec[i478 - 1] + 1) / 2;
+			  if (i478 == 1) ici += ies;
+			  fprintf(output, " SPHERE N. %4d\n", i478);
+			  for (int ic477 = 1; ic477 <= ici; ic477++) {
+				  fprintf(output, " NONTRANSITION LAYER N. %2d , SCALE =  %3s\n", ic477, reference_variable_name.c_str());
+				  for (int jxi476 = 0; jxi476 < number_of_scales; jxi476++) {
+					  double dc0_real = dc0_matrix[ic477 - 1][i478 - 1][jxi476].real();
+					  double dc0_img = dc0_matrix[ic477 - 1][i478 - 1][jxi476].imag();
+					  fprintf(output, "%5d (%12.4lE,%12.4lE)\n", (jxi476 + 1), dc0_real, dc0_img);
+				  }
+			  }
+		  }
+	  }
+	fclose(output);
+}

src/libnptm/Parsers.cpp

+/*! \file Parsers.cpp
+ */
+
+#include <fstream>
+#include <string>
+#include "../include/List.h"
+#include "../include/Parsers.h"
+
+std::string *load_file(std::string file_name, int *count = 0) {
+  std::fstream input_file(file_name.c_str(), std::ios::in);
+  List<std::string> file_lines = List<std::string>();
+  std::string line;
+  if (input_file.is_open()) {
+    getline(input_file, line);
+    file_lines.set(0, line);
+    while (getline(input_file, line)) {
+      file_lines.append(line);
+    }
+    input_file.close();
+  } else {
+	  throw FILE_NOT_FOUND_ERROR;
+  }
+  std::string *array_lines = file_lines.to_array();
+  if (count != 0) *count = file_lines.length();
+  return array_lines;
+}

src/np_sphere.cpp

+/*! \file np_sphere.cpp
+ */
+
+#include <cstdio>
+#include <string>
+#include "include/Configuration.h"
+
+using namespace std;
+
+extern void sphere();
+
+/*! \brief Main program entry point.
+ *
+ * This is the starting point of the execution flow. Here we may choose
+ * how to configure the code, e.g. by loading a legacy configuration file
+ * or some otherwise formatted configuration data set. The code can be
+ * linked to a luncher script or to a GUI oriented application that performs
+ * the configuration and runs the main program.
+ */
+int main(int argc, char **argv) {
+	sphere();
+	return 0;
+}

src/sphere/Makefile

@@ -2,8 +2,11 @@ BUILDDIR=../../build/sphere
 FC=gfortran
 FCFLAGS=-std=legacy -O3
 LFLAGS=
+CXX=g++
+CXXFLAGS=-O0 -ggdb -pg -coverage
+CXXLFLAGS=
 
-all: edfb sph
+all: edfb sph np_sphere
 
 edfb: edfb.o
 	$(FC) $(FCFLAGS) -o $(BUILDDIR)/edfb $(BUILDDIR)/edfb.o $(LFLAGS)
@@ -11,6 +14,24 @@ edfb: edfb.o
 sph: sph.o
 	$(FC) $(FCFLAGS) -o $(BUILDDIR)/sph $(BUILDDIR)/sph.o $(LFLAGS)
 
+np_sphere: $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o
+	$(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_sphere $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o
+
+$(BUILDDIR)/np_sphere.o:
+	$(CXX) $(CXXFLAGS) -c ../np_sphere.cpp -o $(BUILDDIR)/np_sphere.o
+
+$(BUILDDIR)/Commons.o:
+	$(CXX) $(CXXFLAGS) -c ../libnptm/Commons.cpp -o $(BUILDDIR)/Commons.o
+
+$(BUILDDIR)/Configuration.o:
+	$(CXX) $(CXXFLAGS) -c ../libnptm/Configuration.cpp -o $(BUILDDIR)/Configuration.o
+
+$(BUILDDIR)/Parsers.o:
+	$(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o
+
+$(BUILDDIR)/sphere.o:
+	$(CXX) $(CXXFLAGS) -c sphere.cpp -o $(BUILDDIR)/sphere.o
+
 clean:
 	rm -f $(BUILDDIR)/*.o
 

src/sphere/sph.f

@@ -83,9 +83,9 @@ CCC   DIMENSION DC0M(NSPH,NSHL-NTL),XIV(NXI)
  6991 FORMAT('========== JXI =',I3,' ====================')
  6992 FORMAT('********** JTH =',I3,', JPH =',I3,
      1', JTHS =',I3,', JPHS =',I3,' ********************')
-      IR=5
-      IW=6
-      IT=7
+      IR=25
+      IW=26
+      IT=27
       ITIN=17
 CCC
 CCC   OTHER COMMENTS IN FILE GLOBALCOMS
@@ -270,6 +270,7 @@ C
 
       CS0=0.25D0*VK*VK*VK/PIGH
       CALL SSCR0(TFSAS,NSPH,LM,VK,EXRI)
+      PRINT *,"DEBUG: TFSAS =", TFSAS
       SQK=VK*VK*EXDC
       CALL APS(ZPV,LM,NSPH,IOG,RMI,REI,SQK,GAPS)
       CALL RABAS(INPOL,LM,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS)
@@ -455,11 +456,19 @@ CCC  1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
       TQSPE(2,I)=TQSPE(2,I)*PIG2
       TQSPS(1,I)=TQSPS(1,I)*PIG2
       TQSPS(2,I)=TQSPS(2,I)*PIG2
+      PRINT *,"DEBUG: TQSPE(1,",I,") =",TQSPE(1,I)
+      PRINT *,"DEBUG: TQSPE(2,",I,") =",TQSPE(2,I)
+      PRINT *,"DEBUG: TQSPS(1,",I,") =",TQSPS(1,I)
+      PRINT *,"DEBUG: TQSPS(2,",I,") =",TQSPS(2,I)
       GO TO 80
    75 TQSE(1,I)=TQSE(1,I)*PIG2
       TQSE(2,I)=TQSE(2,I)*PIG2
       TQSS(1,I)=TQSS(1,I)*PIG2
       TQSS(2,I)=TQSS(2,I)*PIG2
+      PRINT *,"DEBUG: TQSE(1,",I,") =",TQSE(1,I)
+      PRINT *,"DEBUG: TQSE(2,",I,") =",TQSE(2,I)
+      PRINT *,"DEBUG: TQSS(1,",I,") =",TQSS(1,I)
+      PRINT *,"DEBUG: TQSS(2,",I,") =",TQSS(2,I)
    80 CONTINUE
       RETURN
       END
@@ -575,6 +584,7 @@ CCC  1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
       PIGFSQ=6.4D+1*DACOS(C0)**2
       CFSQ=4.0D0/(PIGFSQ*CCS*CCS)
       NLMM=LM*(LM+2)
+C     PRINT *,"DEBUG: in SSCR2 W(1,1) =",W(1,1)
       DO 14 I=1,NSPH
       IOGI=IOG(I)
       IF(IOGI.LT.I)GO TO 14
@@ -586,14 +596,30 @@ CCC  1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
       DO 10 L=1,LM
       RM=1.0D0/RMI(L,I)
       RE=1.0D0/REI(L,I)
+C     PRINT *,"DEBUG: RM =",RM
+C     PRINT *,"DEBUG: RE =",RE
       LTPO=L+L+1
       DO 10 IM=1,LTPO
       K=K+1
       KE=K+NLMM
+C     PRINT *,"DEBUG: W(",K,",3) =",W(K,3)
+C     PRINT *,"DEBUG: W(",K,",1) =",W(K,1)
+C     PRINT *,"DEBUG: W(",KE,",3) =",W(KE,3)
+C     PRINT *,"DEBUG: W(",KE,",1) =",W(KE,1)
       S11=S11-W(K,3)*W(K,1)*RM-W(KE,3)*W(KE,1)*RE
+C     PRINT *,"DEBUG: S11 =",S11
+C     PRINT *,"DEBUG: W(",K,",4) =",W(K,4)
+C     PRINT *,"DEBUG: W(",KE,",4) =",W(KE,4)
       S21=S21-W(K,4)*W(K,1)*RM-W(KE,4)*W(KE,1)*RE
+C     PRINT *,"DEBUG: W(",K,",2) =",W(K,2)
+C     PRINT *,"DEBUG: W(",KE,",2) =",W(KE,2)
       S12=S12-W(K,3)*W(K,2)*RM-W(KE,3)*W(KE,2)*RE
    10 S22=S22-W(K,4)*W(K,2)*RM-W(KE,4)*W(KE,2)*RE
+C     PRINT *,"DEBUG: CSAM =",CSAM
+C     PRINT *,"DEBUG: S11 =",S11
+C     PRINT *,"DEBUG: S21 =",S21
+C     PRINT *,"DEBUG: S12 =",S12
+C     PRINT *,"DEBUG: S22 =",S22
       SAS(I,1,1)=S11*CSAM
       SAS(I,2,1)=S21*CSAM
       SAS(I,1,2)=S12*CSAM
@@ -644,6 +670,7 @@ CCC  1RMI(LI,NSPH),REI(LI,NSPH),GAPS(NSPH)
      1CG1(LMPML,0,L,1)*(SUMM+SUME+SUEM+SUEE))*COFL
    30 CONTINUE
       GAPS(I)=DREAL(SUM)*COFS
+      PRINT *,"DEBUG: GAPS(",I,") =",GAPS(I)
    40 CONTINUE
       RETURN
       END
@@ -664,6 +691,11 @@ CCC   DIMENSION ZPV(LM,3,2,2)
       ZP=-1.0D0/DSQRT(XD)
       ZPV(L,2,1,2)=ZP
       ZPV(L,2,2,1)=ZP
+      IF(ZP.EQ.C0) GOTO 15
+C      PRINT 9010,L-1,ZP
+C      PRINT 9011,L-1,ZP
+ 9010 FORMAT("DEBUG: zpv[",I1,"][1][0][1] = ",E15.3)
+ 9011 FORMAT("DEBUG: zpv[",I1,"][1][1][0] = ",E15.3)
    15 CONTINUE
       IF(LM.EQ.1)GO TO 30
       DO 20 L=2,LM
@@ -672,6 +704,11 @@ CCC   DIMENSION ZPV(LM,3,2,2)
       ZP=DSQRT(XN/XD)
       ZPV(L,1,1,1)=ZP
       ZPV(L,1,2,2)=ZP
+      IF(ZP.EQ.C0) GOTO 20
+C      PRINT 9012,L-1,ZP
+C      PRINT 9013,L-1,ZP
+ 9012 FORMAT("DEBUG: zpv[",I1,"][0][0][0] = ",E15.3)
+ 9013 FORMAT("DEBUG: zpv[",I1,"][0][1][1] = ",E15.3)
    20 CONTINUE
       LMMO=LM-1
       DO 25 L=1,LMMO
@@ -680,6 +717,11 @@ CCC   DIMENSION ZPV(LM,3,2,2)
       ZP=-DSQRT(XN/XD)
       ZPV(L,3,1,1)=ZP
       ZPV(L,3,2,2)=ZP
+      IF(ZP.EQ.C0) GOTO 25
+C      PRINT 9014,L-1,ZP
+C      PRINT 9015,L-1,ZP
+ 9014 FORMAT("DEBUG: zpv[",I1,"][2][0][0] = ",E15.3)
+ 9015 FORMAT("DEBUG: zpv[",I1,"][2][1][1] = ",E15.3)
    25 CONTINUE
    30 CONTINUE
       RETURN
@@ -739,6 +781,10 @@ CCC   CG1(LMPML,MU,L,M)=CLGO(1,LMP,L;MU,M-MU,M)
       SINT=DSIN(TH)
       COSP=DCOS(PH)
       SINP=DSIN(PH)
+C     PRINT *,"DEBUG: cost =",COST
+C     PRINT *,"DEBUG: sint =",SINT
+C     PRINT *,"DEBUG: cosp =",COSP
+C     PRINT *,"DEBUG: sinp =",SINP
       U(1)=COSP*SINT
       U(2)=SINP*SINT
       U(3)=COST
@@ -782,6 +828,8 @@ CCC   NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
       DO 60 N=1,NSPH
    60 ARG(N)=-ARG(N)
  65   CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
+C     PRINT *,"DEBUG: in WMAMP and calling PWMA with lm =",LM,
+C    1"and iis =",IIS
       CALL PWMA(UP,UN,YLM,INPOL,LM,IIS)
       RETURN
       END
@@ -880,9 +928,11 @@ CCC   NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
    55 ARGS(N)=-COSTS*RZZ(N)
    60 CONTINUE
  75   CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
+C     PRINT *,"DEBUG: in WMASP and calling PWMA"
       CALL PWMA(UP,UN,YLM,INPOL,LM,ISQ)
       IF(IBF.LT.0)RETURN
       CALL SPHAR(COSTS,SINTS,COSPS,SINPS,LM,YLM)
+C     PRINT *,"DEBUG: in WMASP and calling PWMA with IBF =",IBF
       CALL PWMA(UPS,UNS,YLM,INPOL,LM,2)
       RETURN
       END
@@ -913,6 +963,7 @@ CCC   DIMENSION YLM(NLWM+2)
       CM2=.5D0*DCMPLX(UN(1),UN(2))
       CP2=.5D0*DCMPLX(UN(1),-UN(2))
       CZ2=UN(3)
+C     PRINT *,"DEBUG: in PWMA YLM(2) =",YLM(2)
       DO 20 L=1,LW
       LF=L+1
       LFTL=LF*L
@@ -929,26 +980,35 @@ CCC   DIMENSION YLM(NLWM+2)
       CM=DSQRT(X)
       CZ=M
       W(K,ISPO)=DCONJG(CP1*CP*YLM(K+2)+CM1*CM*YLM(K)+CZ1*CZ*YLM(K+1))*CL
+C     IF(ISPO.EQ.1)PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO)
  20   W(K,ISPT)=DCONJG(CP2*CP*YLM(K+2)+CM2*CM*YLM(K)+CZ2*CZ*YLM(K+1))*CL
+C 20  PRINT *,"DEBUG: W(",K,",",ISPT,") =",W(K,ISPT)
       DO 30 K=1,NLWM
       I=K+NLWM
       W(I,ISPO)=UIM*W(K,ISPT)
+C     PRINT *,"DEBUG: W(",I,",",ISPO,") =",W(I,ISPO)
  30   W(I,ISPT)=-UIM*W(K,ISPO)
+C 30  PRINT *,"DEBUG: W(",I,",",ISPT,") =",W(I,ISPT)
       IF(INPOL.EQ.0)GO TO 42
       DO 40 K=1,NLWM
       I=K+NLWM
       C1=(W(K,ISPO)+UIM*W(K,ISPT))*SQRTWI
       C2=(W(K,ISPO)-UIM*W(K,ISPT))*SQRTWI
       W(K,ISPO)=C2
+C     PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO)
       W(I,ISPO)=-C2
+C     PRINT *,"DEBUG: W(",I,",",ISPO,") =",W(I,ISPO)
       W(K,ISPT)=C1
+C     PRINT *,"DEBUG: W(",K,",",ISPT,") =",W(K,ISPT)
  40   W(I,ISPT)=C1
+C 40  PRINT *,"DEBUG: W(",I,",",ISPT,") =",W(I,ISPT)
    42 IF(ISQ.EQ.0)RETURN
       DO 50 I=1,2
       IPT=I+2
       IPIS=I+IS
       DO 50 K=1,NLWMT
  50   W(K,IPT)=DCONJG(W(K,IPIS))
+C 50  PRINT *,"DEBUG: W(",K,",",IPT,") =",W(K,IPT)
       RETURN
       END
       SUBROUTINE ORUNVE(U1,U2,U3,IORTH,TORTH)
@@ -1041,7 +1101,9 @@ CCC  2RMF(LI),DRMF(LI),REF(LI),DREF(LI)
       CCNC=FBI(LPO)*DFB
       CCND=FB(LPO)*DFBI
       RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND)
+C     PRINT *,"DEBUG: gone 60, RMI(",L,",",I,") =",RMI(L,I)
  60   REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
+C 60  PRINT *,"DEBUG: gone 60, REI(",L,",",I,") =",REI(L,I)
       RETURN
    65 DO 80 L=1,LI
       LPO=L+1
@@ -1076,6 +1138,7 @@ cccccccccccccccccccccc
    80 DREF(L)=DY2*SZ+Y2
       CRI=(1.0D0,0.0D0)
       IF(MOD(NSH,2).NE.0)CRI=DC0(IC)/EXDC
+C     PRINT *,"DEBUG: going 90, AREX =",AREX
       DO 90 L=1,LI
       LPO=L+1
       LTPO=LPO+L
@@ -1087,11 +1150,13 @@ cccccccccccccccccccccc
       CCNC=RMF(L)*DFB
       CCND=DRMF(L)*FB(LPO)*SZ*LTPO
       RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND)
+C      PRINT *,"DEBUG: gone 90, RMI(",L,",",I,") =",RMI(L,I)
       CCNA=REF(L)*DFN
       CCNB=DREF(L)*FN(LPO)*SZ*LTPO
       CCNC=REF(L)*DFB
       CCND=DREF(L)*FB(LPO)*SZ*LTPO
  90   REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
+C 90  PRINT *,"DEBUG: gone 90, REI(",L,",",I,") =",REI(L,I)
       RETURN
       END
       SUBROUTINE RKT(NPNTMO,STEP,X,LPO,Y1,Y2,DY1,DY2)
@@ -1433,13 +1498,18 @@ CCC   LL=LM
       PI4=DACOS(0.0D0)*8.0D0
       PI4IRS=1.0D0/DSQRT(PI4)
       X=COSRTH
+C     PRINT *,"DEBUG: X =",X
       Y=DABS(SINRTH)
+C     PRINT *,"DEBUG: Y =",Y
       CLLMO=3.0D0
       CLL=1.5D0
       YTOL=Y
       PLEGN(1)=1.0D0
+C     PRINT *,"DEBUG: PLEGN( 1 ) =",PLEGN(1)
       PLEGN(2)=X*DSQRT(CLLMO)
+C     PRINT *,"DEBUG: PLEGN( 2 ) =",PLEGN(2)
       PLEGN(3)=YTOL*DSQRT(CLL)
+C     PRINT *,"DEBUG: PLEGN( 3 ) =",PLEGN(3)
       SINRMP(1)=SINRPH
       COSRMP(1)=COSRPH
       IF(LL.LT.2)GO TO 30
@@ -1459,14 +1529,17 @@ CCC   LL=LM
       CDM=(LTMO-2)*LS
  10   PLEGN(MPOPK)=PLEGN(MPOPK-L)*X*DSQRT(CN/CD)-
      1PLEGN(MPOPK-LTMO)*DSQRT(CNM/CDM)
+C10   PRINT *,"DEBUG: PLEGN(",MPOPK,") =",PLEGN(MPOPK)
       LPK=L+K
       CLTPO=LTPO
       PLEGN(LPK)=PLEGN(K)*X*DSQRT(CLTPO)
+C     PRINT *,"DEBUG: PLEGN(",LPK,") =",PLEGN(LPK)
       K=LPK+1
       CLT=LTPO-1
       CLL=CLL*(CLTPO/CLT)
       YTOL=YTOL*Y
       PLEGN(K)=YTOL*DSQRT(CLL)
+C     PRINT *,"DEBUG: PLEGN(",K,") =",PLEGN(K)
       SINRMP(L)=SINRPH*COSRMP(LMO)+COSRPH*SINRMP(LMO)
    20 COSRMP(L)=COSRPH*COSRMP(LMO)-SINRPH*SINRMP(LMO)
    30 L=0
@@ -1475,14 +1548,17 @@ CCC   LL=LM
       L0Y=K+1
       L0P=K/2+1
       YLM(L0Y)=PI4IRS*PLEGN(L0P)
+C     PRINT *, "DEBUG: YLM(",L0Y,") =",YLM(L0Y)
       GO TO 45
    44 LMP=L0P+M
       SAVE=PI4IRS*PLEGN(LMP)
       LMY=L0Y+M
       YLM(LMY)=SAVE*DCMPLX(COSRMP(M),SINRMP(M))
       IF(MOD(M,2).NE.0)YLM(LMY)=-YLM(LMY)
+C     PRINT *, "DEBUG: YLM(",LMY,") =",YLM(LMY)
       LMY=L0Y-M
       YLM(LMY)=SAVE*DCMPLX(COSRMP(M),-SINRMP(M))
+C     PRINT *, "DEBUG: YLM(",LMY,") =",YLM(LMY)
  45   IF(M.GE.L)GO TO 47
       M=M+1
       GO TO 44

src/sphere/sphere.cpp

+#include <cstdio>
+#include <fstream>
+#include <string>
+#include <complex>
+#include "../include/Configuration.h"
+#include "../include/sph_subs.h"
+
+using namespace std;
+
+//! \brief C++ implementation of SPH
+void sphere() {
+	complex<double> arg, s0, tfsas;
+	complex<double> **tqspe, **tqsps;
+	double **tqse, **tqss;
+	double *argi, *args, *gaps;
+	double th, ph;
+	printf("INFO: making legacy configuration ...\n");
+	ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/sphere/DEDFB");
+	conf->write_formatted("c_OEDFB");
+	conf->write_binary("c_TEDF");
+	delete conf;
+	printf("INFO: reading binary configuration ...\n");
+	ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF");
+	GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/sphere/DSPH");
+	if (sconf->number_of_spheres == gconf->number_of_spheres) {
+		int isq, ibf;
+		double cost, sint, cosp, sinp;
+		double costs, sints, cosps, sinps;
+		double scan;
+		double *duk = new double[3];
+		double *u = new double[3];
+		double *us = new double[3];
+		double *un = new double[3];
+		double *uns = new double[3];
+		double *up = new double[3];
+		double *ups = new double[3];
+		double *upmp = new double[3];
+		double *upsmp = new double[3];
+		double *unmp = new double[3];
+		double *unsmp = new double[3];
+		double **cmul = new double*[4];
+		double **cmullr = new double*[4];
+		for (int i = 0; i < 4; i++) {
+			cmul[i] = new double[4];
+			cmullr[i] = new double[4];
+		}
+		double frx = 0.0, fry = 0.0, frz = 0.0;
+		double cfmp, cfsp, sfmp, sfsp;
+		complex<double> *vint = new complex<double>[16];
+		int jw;
+		int nsph = gconf->number_of_spheres;
+		C1 *c1 = new C1(nsph, gconf->l_max);
+		for (int i = 0; i < nsph; i++) {
+			c1->iog[i] = sconf->iog_vec[i];
+			c1->nshl[i] = sconf->nshl_vec[i];
+			c1->ros[i] = sconf->radii_of_spheres[i];
+		}
+		for (int i = 0; i < gconf->l_max; i++) {
+			c1->rmi[i][0] = complex<double>(0.0, 0.0);
+			c1->rmi[i][1] = complex<double>(0.0, 0.0);
+			c1->rei[i][0] = complex<double>(0.0, 0.0);
+			c1->rei[i][1] = complex<double>(0.0, 0.0);
+		}
+		C2 *c2 = new C2(nsph, 5, gconf->npnt, gconf->npntts);
+		argi = new double[1];
+		args = new double[1];
+		gaps = new double[2];
+		FILE *output = fopen("c_OSPH", "w");
+		fprintf(output, " READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM\n");
+		fprintf(
+				output,
+				" %5d%5d%5d%5d%5d%5d\n",
+				gconf->number_of_spheres,
+				gconf->l_max,
+				gconf->in_pol,
+				gconf->npnt,
+				gconf->npntts,
+				gconf->meridional_type
+		);
+		fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
+		fprintf(
+				output,
+				"  %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",
+				gconf->in_theta_start,
+				gconf->in_theta_step,
+				gconf->in_theta_end,
+				gconf->sc_theta_start,
+				gconf->sc_theta_step,
+				gconf->sc_theta_end
+		);
+		fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
+		fprintf(
+				output,
+				"  %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",
+				gconf->in_phi_start,
+				gconf->in_phi_step,
+				gconf->in_phi_end,
+				gconf->sc_phi_start,
+				gconf->sc_phi_step,
+				gconf->sc_phi_end
+		);
+		fprintf(output, " READ(IR,*)JWTM\n");
+		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");
+		fprintf(output, "  READ(ITIN)(XIV(I),I=1,NXI)\n");
+		fprintf(output, "  READ(ITIN)NSHL(I),ROS(I)\n");
+		fprintf(output, "  READ(ITIN)(RCF(I,NS),NS=1,NSH)\n \n");
+		double sml = 1.0e-3;
+		int nth = 0, nph = 0;
+		if (gconf->in_theta_step != 0.0)
+			nth = int((gconf->in_theta_end - gconf->in_theta_start) / gconf->in_theta_step + sml);
+		nth += 1;
+		if (gconf->in_phi_step != 0.0)
+			nph = int((gconf->in_phi_end - gconf->in_phi_start) / gconf->in_phi_step + sml);
+		nph += 1;
+		int nths = 0, nphs = 0;
+		double thsca;
+		if (gconf->meridional_type > 1) { // ISAM > 1, fixed scattering angle
+			nths = 1;
+			thsca = gconf->sc_theta_start - gconf->in_theta_start;
+		} else { //ISAM <= 1
+			if (gconf->in_theta_step != 0.0)
+				nths = int((gconf->sc_theta_end - gconf->sc_theta_start) / gconf->sc_theta_step + sml);
+			nths += 1;
+			if (gconf->meridional_type == 1) { // ISAM = 1
+				nphs = 1;
+			} else { // ISAM < 1
+				if (gconf->sc_phi_step != 0.0)
+					nphs = int((gconf->sc_phi_end - gconf->sc_phi_start) / gconf->sc_phi_step + sml);
+				nphs += 1;
+			}
+		}
+		int nk = nth * nph;
+		int nks = nths * nphs;
+		int nkks = nk * nks;
+		double th1 = gconf->in_theta_start;
+		double ph1 = gconf->in_phi_start;
+		double ths1 = gconf->sc_theta_start;
+		double phs1 = gconf->sc_phi_start;
+		const double half_pi = acos(0.0);
+		const double pi = 2.0 * half_pi;
+		double gcs = 0.0;
+		for (int i116 = 1; i116 <= nsph; i116++) {
+			int iogi = c1->iog[i116 - 1];
+			if (iogi >= i116) {
+				double gcss = pi * c1->ros[i116 - 1] * c1->ros[i116 - 1];
+				c1->gcsv[i116 - 1] = gcss;
+				int nsh = c1->nshl[i116 - 1];
+				for (int j115 = 1; j115 <= nsh; j115++) {
+					c1->rc[i116 - 1][j115 - 1] = sconf->rcf[i116 - 1][j115 - 1] * c1->ros[i116 - 1];
+				}
+			}
+			gcs += c1->gcsv[iogi - 1];
+		}
+		double ****zpv = new double***[gconf->l_max]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
+		for (int zi = 0; zi < gconf->l_max; 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];
+				}
+			}
+		}
+		thdps(gconf->l_max, zpv);
+		double exri = sqrt(sconf->exdc);
+		fprintf(output, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
+		fstream tppoan;
+		tppoan.open("c_TPPOAN", ios::binary|ios::out);
+		if (tppoan.is_open()) {
+			int imode = 10;
+			tppoan.write(reinterpret_cast<char *>(&imode), sizeof(int));
+			tppoan.write(reinterpret_cast<char *>(&(gconf->meridional_type)), sizeof(int));
+			tppoan.write(reinterpret_cast<char *>(&(gconf->in_pol)), sizeof(int));
+			tppoan.write(reinterpret_cast<char *>(&(sconf->number_of_scales)), sizeof(int));
+			tppoan.write(reinterpret_cast<char *>(&nth), sizeof(int));
+			tppoan.write(reinterpret_cast<char *>(&nph), sizeof(int));
+			tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
+			tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
+
+			for (int nsi = 0; nsi < nsph; nsi++)
+				tppoan.write(reinterpret_cast<char *>(&(sconf->iog_vec[nsi])), sizeof(int));
+			if (gconf->in_pol == 0) fprintf(output, "   LIN\n");
+			else fprintf(output, "  CIRC\n");
+			fprintf(output, " \n");
+			double wn = sconf->wp / 3.0e8;
+			double sqsfi = 1.0;
+			double vk, vkarg;
+			if (sconf->idfc < 0) {
+				vk = sconf->xip * wn;
+				fprintf(output, "  VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
+				fprintf(output, " \n");
+			}
+			for (int jxi488 = 1; jxi488 <= sconf->number_of_scales; jxi488++) {
+				fprintf(output, "========== JXI =%3d ====================\n", jxi488);
+				double xi = sconf->scale_vec[jxi488 - 1];
+				if (sconf->idfc >= 0) {
+					vk = xi * wn;
+					vkarg = vk;
+					fprintf(output, "  VK=%15.7lE, XI=%15.7lE\n", xi, vk);
+				} else { // IDFC < 0
+					vkarg = xi * vk;
+					sqsfi = 1.0 / (xi * xi);
+					fprintf(output, "  XI=%15.7lE\n", xi);
+				}
+				tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
+				for (int i132 = 1; i132 <= nsph; i132++) {
+					int iogi = sconf->iog_vec[i132 - 1];
+					if (iogi != i132) {
+						for (int l123 = 1; l123 <= gconf->l_max; l123++) {
+							c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];
+							c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];
+						}
+						continue; // i132
+					}
+					int nsh = sconf->nshl_vec[i132 - 1];
+					int ici = (nsh + 1) / 2;
+					if (sconf->idfc == 0) {
+						for (int ic = 0; ic < ici; ic++)
+							c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][0]; // WARNING: IDFC=0 is not tested!
+					} else { // IDFC != 0
+						if (jxi488 == 1) {
+							for (int ic = 0; ic < ici; ic++) {
+								c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][jxi488 - 1];
+							}
+						}
+					}
+					if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
+					int jer = 0;
+					int lcalc = 0;
+					dme(
+							gconf->l_max, i132, gconf->npnt, gconf->npntts, vkarg,
+							sconf->exdc, exri, c1, c2, jer, lcalc, arg
+					);
+					if (jer != 0) {
+						fprintf(output, "  STOP IN DME\n");
+						fprintf(output, "  AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n", jer, lcalc, arg.real(), arg.imag());
+						tppoan.close();
+						fclose(output);
+						return;
+					}
+				} // i132
+				if (sconf->idfc >= 0 and nsph == 1 and jxi488 == gconf->jwtm) {
+					// This is the condition that writes the transition matrix to output.
+					int is = 1111;
+					fstream ttms;
+					ttms.open("c_TTMS", ios::binary | ios::out);
+					if (ttms.is_open()) {
+						ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
+						ttms.write(reinterpret_cast<char *>(&(gconf->l_max)), sizeof(int));
+						ttms.write(reinterpret_cast<char *>(&vk), sizeof(double));
+						ttms.write(reinterpret_cast<char *>(&exri), sizeof(double));
+						for (int lmi = 0; lmi < gconf->l_max; lmi++) {
+							complex<double> element1 = -1.0 / c1->rmi[0][lmi];
+							complex<double> element2 = -1.0 / c1->rei[0][lmi];
+							ttms.write(reinterpret_cast<char *>(&element1), sizeof(complex<double>));
+							ttms.write(reinterpret_cast<char *>(&element2), sizeof(complex<double>));
+						}
+						ttms.write(reinterpret_cast<char *>(&(sconf->radii_of_spheres[0])), sizeof(double));
+						ttms.close();
+					} else { // Failed to open output file. Should never happen.
+						printf("ERROR: could not open TTMS file.\n");
+						tppoan.close();
+						fclose(output);
+					}
+				}
+				double cs0 = 0.25 * vk * vk * vk / half_pi;
+				//printf("DEBUG: cs0 = %lE\n", cs0);
+				sscr0(tfsas, nsph, gconf->l_max, vk, exri, c1);
+				printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag());
+				double sqk = vk * vk * sconf->exdc;
+				aps(zpv, gconf->l_max, nsph, c1, sqk, gaps);
+				tqse = new double*[2];
+				tqss = new double*[2];
+				tqspe = new std::complex<double>*[2];
+				tqsps = new std::complex<double>*[2];
+				for (int ti = 0; ti < 2; ti++) {
+					tqse[ti] = new double[2];
+					tqss[ti] = new double[2];
+					tqspe[ti] = new std::complex<double>[2];
+					tqsps[ti] = new std::complex<double>[2];
+				}
+				rabas(gconf->in_pol, gconf->l_max, nsph, c1, tqse, tqspe, tqss, tqsps);
+				for (int i170 = 1; i170 <= nsph; i170++) {
+					if (c1->iog[i170 - 1] >= i170) {
+						double albeds = c1->sscs[i170 - 1] / c1->sexs[i170 - 1];
+						c1->sqscs[i170 - 1] *= sqsfi;
+						c1->sqabs[i170 - 1] *= sqsfi;
+						c1->sqexs[i170 - 1] *= sqsfi;
+						fprintf(output, "     SPHERE %2d\n", i170);
+						if (c1->nshl[i170 - 1] != 1) {
+							fprintf(output, "  SIZE=%15.7lE\n", c2->vsz[i170 - 1]);
+						} else {
+							fprintf(
+									output,
+									"  SIZE=%15.7lE, REFRACTIVE INDEX=(%15.7lE,%15.7lE)\n",
+									c2->vsz[i170 -1],
+									c2->vkt[i170 - 1].real(),
+									c2->vkt[i170 - 1].imag()
+							);
+						}
+						fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+						fprintf(
+								output,
+								" %14.7lE%15.7lE%15.7lE%15.7lE\n",
+								c1->sscs[i170 - 1], c1->sabs[i170 - 1],
+								c1->sexs[i170 - 1], albeds
+						);
+						fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
+						fprintf(
+								output,
+								" %14.7lE%15.7lE%15.7lE\n",
+								c1->sqscs[i170 - 1], c1->sqabs[i170 - 1],
+								c1->sqexs[i170 - 1]
+						);
+						fprintf(output, "  FSAS=%15.7lE%15.7lE\n", c1->fsas[i170 - 1].real(), c1->fsas[i170 - 1].imag());
+						double csch = 2.0 * vk * sqsfi / c1->gcsv[i170 -1];
+						//printf("DEBUG: csch = %lE\n", csch);
+						s0 = c1->fsas[i170 - 1] * exri;
+						//printf("DEBUG: s0 = (%lE,%lE)\n", s0.real(), s0.imag());
+						double qschu = csch * s0.imag();
+						//printf("DEBUG: qschu = %lE\n", qschu);
+						double pschu = csch * s0.real();
+						//printf("DEBUG: pschu = %lE\n", pschu);
+						double s0mag = cs0 * sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag()));
+						//printf("DEBUG: s0mag = %lE\n", s0mag);
+						fprintf(
+								output,
+								"  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+								qschu, pschu, s0mag
+						);
+						double rapr = c1->sexs[i170 - 1] - gaps[i170 - 1];
+						double cosav = gaps[i170 - 1] / c1->sscs[i170 - 1];
+						fprintf(output, "  COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
+						fprintf(output, "  IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 1, tqse[0][i170 - 1], tqss[0][i170 - 1]);
+						fprintf(output, "  IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170 - 1], tqss[1][i170 - 1]);
+						tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170 - 1])), sizeof(double));
+						tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170 - 1])), sizeof(double));
+						double val = tqspe[0][i170 - 1].real();
+						tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+						val = tqspe[0][i170 - 1].imag();
+						tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+						val = tqsps[0][i170 - 1].real();
+						tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+						val = tqsps[0][i170 - 1].imag();
+						tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+						tppoan.write(reinterpret_cast<char *>(&(tqse[1][i170 - 1])), sizeof(double));
+						tppoan.write(reinterpret_cast<char *>(&(tqss[1][i170 - 1])), sizeof(double));
+						val = tqspe[1][i170 - 1].real();
+						tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+						val = tqspe[1][i170 - 1].imag();
+						tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+						val = tqsps[1][i170 - 1].real();
+						tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+						val = tqsps[1][i170 - 1].imag();
+						tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+					} // End if iog[i170 - 1] >= i170
+				} // i170 loop
+				if (nsph != 1) {
+					fprintf(output, "  FSAT=(%15.7lE,%15.7lE)\n", tfsas.real(), tfsas.imag());
+					double csch = 2.0 * vk * sqsfi / gcs;
+					s0 = tfsas * exri;
+					double qschu = csch * s0.imag();
+					double pschu = csch * s0.real();
+					double s0mag = cs0 * sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag()));
+					fprintf(
+							output,
+							"  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+							qschu, pschu, s0mag
+					);
+				}
+				th = th1;
+				for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP parallelizable section
+					ph = ph1;
+					for (int jph484 = 1; jph484 <= nph; jph484++) {
+						bool goto182 = (nk == 1) && (jxi488 > 1);
+						if (!goto182) {
+							upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
+						}
+						if (gconf->meridional_type >= 0) {
+							wmamp(0, cost, sint, cosp, sinp, gconf->in_pol, gconf->l_max, 0, nsph, argi, u, upmp, unmp, c1);
+							for (int i183 = 0; i183 < nsph; i183++) {
+								double rapr = c1->sexs[i183] - gaps[i183];
+								frx = rapr * u[0];
+								fry = rapr * u[1];
+								frz = rapr * u[2];
+							}
+							jw = 1;
+						}
+						double thsl = ths1;
+						double phsph = 0.0;
+						for (int jths482 = 1; jths482 <= nths; jths482++) {
+							double ths = thsl;
+							int icspnv = 0;
+							if (gconf->meridional_type > 1) ths = th + thsca;
+							if (gconf->meridional_type >= 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;
+							}
+							double phs = phs1;
+							for (int jphs480 = 1; jphs480 <= nphs; jphs480++) {
+								if (gconf->meridional_type >= 1) {
+									phs = ph + phsph;
+									if (phs >= 360.0) phs -= 360.0;
+								}
+								bool goto190 = (nks == 1) && ((jxi488 > 1) || (jth486 > 1) || (jph484 > 1));
+								if (!goto190) {
+									upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
+									if (gconf->meridional_type >= 0) {
+										wmamp(2, costs, sints, cosps, sinps, gconf->in_pol, gconf->l_max, 0, nsph, args, us, upsmp, unsmp, c1);
+									}
+								}
+								if (nkks != 0 || jxi488 == 1) {
+									upvsp(u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns, duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp);
+									if (gconf->meridional_type < 0) {
+										wmasp(
+												cost, sint, cosp, sinp, costs, sints, cosps, sinps,
+												u, up, un, us, ups, uns, isq, ibf, gconf->in_pol,
+												gconf->l_max, 0, nsph, argi, args, c1
+										);
+									}
+									for (int i193 = 1; i193 <= 3; i193++) {
+										un[i193 - 1] = unmp[i193 - 1];
+										uns[i193 - 1] = unsmp[i193 - 1];
+									}
+								}
+								if (gconf->meridional_type < 0) jw = 1;
+								tppoan.write(reinterpret_cast<char *>(&th), sizeof(double));
+								tppoan.write(reinterpret_cast<char *>(&ph), sizeof(double));
+								tppoan.write(reinterpret_cast<char *>(&ths), sizeof(double));
+								tppoan.write(reinterpret_cast<char *>(&phs), sizeof(double));
+								tppoan.write(reinterpret_cast<char *>(&scan), sizeof(double));
+								if (jw != 0) {
+									jw = 0;
+									tppoan.write(reinterpret_cast<char *>(&(u[0])), sizeof(double));
+									tppoan.write(reinterpret_cast<char *>(&(u[1])), sizeof(double));
+									tppoan.write(reinterpret_cast<char *>(&(u[2])), sizeof(double));
+								}
+								fprintf(
+										output,
+										"********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n",
+										jth486, jph484, jths482, jphs480
+								);
+								fprintf(
+										output,
+										"  TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n",
+										th, ph, ths, phs
+								);
+								fprintf(output, "  SCAND=%10.3lE\n", scan);
+								fprintf(output, "  CFMP=%15.7lE, SFMP=%15.7lE\n", cfmp, sfmp);
+								fprintf(output, "  CFSP=%15.7lE, SFSP=%15.7lE\n", cfsp, sfsp);
+								if (gconf->meridional_type >= 0) {
+									fprintf(output, "  UNI=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
+									fprintf(output, "  UNS=(%12.5lE,%12.5lE,%12.5lE)\n", uns[0], uns[1], uns[2]);
+								} else {
+									fprintf(output, "  UN=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
+								}
+								sscr2(nsph, gconf->l_max, vk, exri, c1);
+								for (int ns226 = 1; ns226 <= nsph; ns226++) {
+									fprintf(output, "     SPHERE %2d\n", ns226);
+									fprintf(
+											output, "  SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
+											c1->sas[ns226 - 1][0][0].real(), c1->sas[ns226 - 1][0][0].imag(),
+											c1->sas[ns226 - 1][1][0].real(), c1->sas[ns226 - 1][1][0].imag()
+									);
+									fprintf(
+											output, "  SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
+											c1->sas[ns226 - 1][0][1].real(), c1->sas[ns226 - 1][0][1].imag(),
+											c1->sas[ns226 - 1][1][1].real(), c1->sas[ns226 - 1][1][1].imag()
+									);
+									if (jths482 == 1 && jphs480 == 1)
+										fprintf(
+												output, "  Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
+												frx, fry, frz
+										);
+									for (int i225 = 0; i225 < 16; i225++)
+										vint[i225] = c1->vints[ns226 - 1][i225];
+									mmulc(vint, cmullr, cmul);
+									fprintf(output, "  MULS\n        ");
+									for (int imul = 0; imul < 4; imul++) {
+										for (int jmul = 0; jmul < 4; jmul++) {
+											fprintf(output, "%15.7lE", cmul[imul][jmul]);
+										}
+										if (imul < 3) fprintf(output, "\n        ");
+										else fprintf(output, "\n");
+									}
+									fprintf(output, "  MULSLR\n        ");
+									for (int imul = 0; imul < 4; imul++) {
+										for (int jmul = 0; jmul < 4; jmul++) {
+											fprintf(output, "%15.7lE", cmullr[imul][jmul]);
+										}
+										if (imul < 3) fprintf(output, "\n        ");
+										else fprintf(output, "\n");
+									}
+									for (int vi = 0; vi < 16; vi++) {
+										double value = vint[vi].real();
+										tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+										value = vint[vi].imag();
+										tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+									}
+									for (int imul = 0; imul < 4; imul++) {
+										for (int jmul = 0; jmul < 4; jmul++) {
+											tppoan.write(reinterpret_cast<char *>(&(cmul[imul][jmul])), sizeof(double));
+										}
+									}
+								} // ns226 loop
+								if (gconf->meridional_type < 1) phs += gconf->sc_phi_step;
+							} // jphs480 loop
+							if (gconf->meridional_type <= 1) thsl += gconf->sc_theta_step;
+						} // jths482 loop
+						ph += gconf->in_phi_step;
+					} // jph484 loop on elevation
+					th += gconf->in_theta_step;
+				} // jth486 loop on azimuth
+			} //jxi488 loop on scales
+			tppoan.close();
+		} else { // In case TPPOAN could not be opened. Should never happen.
+			printf("ERROR: failed to open TPPOAN file.\n");
+		}
+		fclose(output);
+		delete c1;
+		delete c2;
+		delete[] duk;
+		delete[] u;
+		delete[] us;
+		delete[] un;
+		delete[] uns;
+		delete[] up;
+		delete[] ups;
+		delete[] upmp;
+		delete[] upsmp;
+		delete[] unmp;
+		delete[] unsmp;
+		delete[] vint;
+		delete[] argi;
+		delete[] args;
+		delete[] gaps;
+		for (int i = 0; i < 4; i++) {
+			delete[] cmul[i];
+			delete[] cmullr[i];
+		}
+		delete[] cmul;
+		delete[] cmullr;
+	} else { // NSPH mismatch between geometry and scatterer configurations.
+		throw UnrecognizedConfigurationException(
+				"Inconsistent geometry and scatterer configurations."
+		);
+	}
+}