diff --git a/.gitignore b/.gitignore
index ecb83700cb4f01a8442484d33b4b02010c8592a6..db0984b6c6b0f612e3c708ac3d67bd154628173a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,4 @@ build/cluster/*
build/sphere/*
build/trapping/*
doc/build/*
+
diff --git a/src/Makefile b/src/Makefile
index ac867ef74084ce78f2e17cff9b63ca1b60585e48..544534f8a4ef39c6e646749269d7ad3b732cc37c 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,8 +1,13 @@
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)
diff --git a/src/README.md b/src/README.md
index 29f340e1f259c9566529d369e6e9dc3923058a8d..83c9bcbc7082cbb0bd38dd4c20536ab6b7e5364d 100644
--- a/src/README.md
+++ b/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`
diff --git a/src/include/Commons.h b/src/include/Commons.h
new file mode 100644
index 0000000000000000000000000000000000000000..088dcc74771f967511b579964e8278000db1f3e3
--- /dev/null
+++ b/src/include/Commons.h
@@ -0,0 +1,131 @@
+/*! \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
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
new file mode 100644
index 0000000000000000000000000000000000000000..fdc36184be5f58f321b221b8752113d210583349
--- /dev/null
+++ b/src/include/Configuration.h
@@ -0,0 +1,321 @@
+/*! \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
diff --git a/src/include/List.h b/src/include/List.h
index 047616aaa574c6b9159600b043bdc2c46fc75065..c0f859922f6238be2b193074b1c82b15046e5973 100644
--- a/src/include/List.h
+++ b/src/include/List.h
@@ -1,9 +1,51 @@
/*! \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
diff --git a/src/include/Parsers.h b/src/include/Parsers.h
new file mode 100644
index 0000000000000000000000000000000000000000..a4523de6d03c78d8efe2467e5e37b067ddcbe1e3
--- /dev/null
+++ b/src/include/Parsers.h
@@ -0,0 +1,34 @@
+/*! \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_ */
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
new file mode 100644
index 0000000000000000000000000000000000000000..3c562fb1f586d34b37481c1317daa491f0ed4ca5
--- /dev/null
+++ b/src/include/sph_subs.h
@@ -0,0 +1,1491 @@
+/*! \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_ */
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..809337ea7adb2475af3d0325d3cbd1b139c025d6
--- /dev/null
+++ b/src/libnptm/Commons.cpp
@@ -0,0 +1,96 @@
+/*! \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;
+}
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6bb02bddd092ed30b02aebf2ef88b0e3f524645a
--- /dev/null
+++ b/src/libnptm/Configuration.cpp
@@ -0,0 +1,778 @@
+/*! \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);
+}
diff --git a/src/libnptm/Parsers.cpp b/src/libnptm/Parsers.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a327df810c95b58e3c0e94cfe207f104d482a5bd
--- /dev/null
+++ b/src/libnptm/Parsers.cpp
@@ -0,0 +1,26 @@
+/*! \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;
+}
diff --git a/src/np_sphere.cpp b/src/np_sphere.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ac985ac4e8b292acd06534aad948586a61d39880
--- /dev/null
+++ b/src/np_sphere.cpp
@@ -0,0 +1,23 @@
+/*! \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;
+}
diff --git a/src/sphere/Makefile b/src/sphere/Makefile
index 41d1d3342ed431507d63447b539fe3813fb2d827..c42d9b04f701cfb823a0c443581c96c23207d920 100644
--- a/src/sphere/Makefile
+++ b/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
diff --git a/src/sphere/sph.f b/src/sphere/sph.f
index 5eb62c2fa02e53689cd6a79063d48c98c026a1de..bebacd746a8dac1c960f0c6eb409db89425a3b75 100644
--- a/src/sphere/sph.f
+++ b/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
@@ -781,7 +827,9 @@ CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
55 IF(IIS.NE.2)GO TO 65
DO 60 N=1,NSPH
60 ARG(N)=-ARG(N)
- 65 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
+ 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
@@ -879,10 +927,12 @@ CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
GO TO 60
55 ARGS(N)=-COSTS*RZZ(N)
60 CONTINUE
- 75 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
+ 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
- 20 W(K,ISPT)=DCONJG(CP2*CP*YLM(K+2)+CM2*CM*YLM(K)+CZ2*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)
- 30 W(I,ISPT)=-UIM*W(K,ISPO)
+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
- 40 W(I,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))
+ 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)
@@ -1026,7 +1086,7 @@ CCC 2RMF(LI),DRMF(LI),REF(LI),DREF(LI)
RETURN
42 DO 43 J=1,LIPT
FB(J)=RFJ(J)
- 43 FN(J)=RFN(J)
+ 43 FN(J)=RFN(J)
IF(NSH.GT.1)GO TO 65
CRI=DC0(1)/EXDC
DO 60 L=1,LI
@@ -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)
- 60 REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*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)
+ 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
@@ -1457,16 +1527,19 @@ CCC LL=LM
CD=LS
CNM=LTPO*(LMO+M)*(L-MPO)
CDM=(LTMO-2)*LS
- 10 PLEGN(MPOPK)=PLEGN(MPOPK-L)*X*DSQRT(CN/CD)-
+ 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
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5d92b19255106177f5dc30201704280ecc4a3e90
--- /dev/null
+++ b/src/sphere/sphere.cpp
@@ -0,0 +1,558 @@
+#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."
+ );
+ }
+}