diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 4c08ba9e21ff7cdb82e15175d06a9eaf97fbabf6..ba2c38a4d3d01f47dadc5656ddc230db74729fc4 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -24,10 +24,25 @@ using namespace std;
*/
void cluster(string config_file, string data_file, string output_path) {
printf("INFO: making legacy configuration...");
- ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(config_file);
- sconf->write_formatted(output_path + "/c_OEDFB_clu");
- sconf->write_binary(output_path + "/c_TEDF_clu");
- GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(data_file);
+ ScattererConfiguration *sconf = NULL;
+ try {
+ sconf = ScattererConfiguration::from_dedfb(config_file);
+ } catch(const OpenConfigurationFileException &ex) {
+ printf("\nERROR: failed to open scatterer configuration file.\n");
+ printf("FILE: %s\n", ex.what());
+ exit(1);
+ }
+ sconf->write_formatted(output_path + "/c_OEDFB");
+ sconf->write_binary(output_path + "/c_TEDF");
+ GeometryConfiguration *gconf = NULL;
+ try {
+ gconf = GeometryConfiguration::from_legacy(data_file);
+ } catch (const OpenConfigurationFileException &ex) {
+ printf("\nERROR: failed to open geometry configuration file.\n");
+ printf("FILE: %s\n", ex.what());
+ if (sconf) delete sconf;
+ exit(1);
+ }
printf(" done.\n");
if (sconf->number_of_spheres == gconf->number_of_spheres) {
// Shortcuts to variables stored in configuration objects
diff --git a/src/include/Configuration.h b/src/include/Configuration.h
index e784fa38bf0bd7b2f681102f0cbba5bbb7628028..c23ac6924cf3595fef13f4e72f76c5e85aaddd87 100644
--- a/src/include/Configuration.h
+++ b/src/include/Configuration.h
@@ -13,24 +13,23 @@
*/
class OpenConfigurationFileException: public std::exception {
protected:
- //! \brief Name of the file that was accessed.
- std::string file_name;
+ //! \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 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 message.
+ */
+ virtual const char* what() const throw() {
+ return file_name.c_str();
+ }
};
/**
@@ -38,21 +37,21 @@ public:
*/
class UnrecognizedConfigurationException: public std::exception {
protected:
- //! Description of the problem.
- std::string message;
+ //! 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 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();
+ }
};
/**
@@ -64,97 +63,97 @@ public:
* fields and their polarization properties.
*/
class GeometryConfiguration {
- //! Temporary work-around to allow cluster() and sphere() peeking in.
- friend void cluster(std::string, std::string, std::string);
- friend void sphere(std::string, std::string, std::string);
+ //! Temporary work-around to allow cluster() and sphere() peeking in.
+ friend void cluster(std::string, std::string, std::string);
+ friend void sphere(std::string, std::string, std::string);
protected:
- //! \brief Number of spherical components.
- int number_of_spheres;
- //! \brief Maximum expansion order of angular momentum.
- int l_max;
- //! \brief QUESTION: definition?
- int li;
- //! \brief QUESTION: definition?
- int le;
- //! \brief QUESTION: definition?
- int mxndm;
- //! \brief QUESTION: definition?
- int iavm;
- //! \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;
+ //! \brief Number of spherical components.
+ int number_of_spheres;
+ //! \brief Maximum expansion order of angular momentum.
+ int l_max;
+ //! \brief QUESTION: definition?
+ int li;
+ //! \brief QUESTION: definition?
+ int le;
+ //! \brief QUESTION: definition?
+ int mxndm;
+ //! \brief QUESTION: definition?
+ int iavm;
+ //! \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 li: `int`
- * \param le: `int`
- * \param mxndm: `int`
- * \param iavm: `int`
- * \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(
+ /*! \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 li: `int`
+ * \param le: `int`
+ * \param mxndm: `int`
+ * \param iavm: `int`
+ * \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,
int li, int le, int mxndm, int iavm,
double *x, double *y, double *z,
@@ -163,24 +162,24 @@ public:
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 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 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);
};
/**
@@ -190,147 +189,147 @@ public:
* data to describe the scatterer properties.
*/
class ScattererConfiguration {
- //! Temporary work-around to allow cluster() and sphere() peeking in.
- friend void cluster(std::string, std::string, std::string);
- friend void sphere(std::string, std::string, std::string);
+ //! Temporary work-around to allow cluster() and sphere() peeking in.
+ friend void cluster(std::string, std::string, std::string);
+ friend void sphere(std::string, std::string, std::string);
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;
+ //! \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 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 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 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 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 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 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);
+ /*! \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/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index 4f52c8f78dd66ae848465795de14178e56b4297f..96d055a11f555c217044cec6802461dd465effd7 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -64,8 +64,9 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(string file_name) {
smatch m;
try {
file_lines = load_file(file_name, &num_lines);
- } catch (exception &ex) {
- throw OpenConfigurationFileException(file_name);
+ } catch (...) {
+ OpenConfigurationFileException ex(file_name);
+ throw ex;
}
int _nsph = 0, _lm = 0, _in_pol = 0, _npnt = 0, _npntts = 0, _isam = 0;
int _li = 0, _le = 0, _mxndm = 0, _iavm = 0;
@@ -263,7 +264,7 @@ ScattererConfiguration* ScattererConfiguration::from_binary(string file_name, st
input.read(reinterpret_cast<char *>(&nxi), sizeof(int));
try {
xi_vec = new double[nxi]();
- } catch (bad_alloc &ex) {
+ } catch (const bad_alloc &ex) {
throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of scales " + nxi);
}
for (int i = 0; i < nxi; i++) {
@@ -283,7 +284,7 @@ ScattererConfiguration* ScattererConfiguration::from_binary(string file_name, st
if (max_ici < (nsh + 1) / 2) max_ici = (nsh + 1) / 2;
try {
rcf_vector[i115 - 1] = new double[nsh]();
- } catch (bad_alloc &ex) {
+ } catch (const bad_alloc &ex) {
throw UnrecognizedConfigurationException("Wrong parameter set: invalid number of layers " + nsh);
}
for (int nsi = 0; nsi < nsh; nsi++) {
@@ -313,7 +314,8 @@ ScattererConfiguration* ScattererConfiguration::from_binary(string file_name, st
}
input.close();
} else { // Opening of the input file did not succeed.
- throw OpenConfigurationFileException(file_name);
+ OpenConfigurationFileException ex(file_name);
+ throw ex;
}
} else { // A different binary format was chosen.
//TODO: this part is not yet implemented.
@@ -346,8 +348,9 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
smatch m;
try {
file_lines = load_file(dedfb_file_name, &num_lines);
- } catch (exception &ex) {
- throw OpenConfigurationFileException(dedfb_file_name);
+ } catch (...) {
+ OpenConfigurationFileException ex(dedfb_file_name);
+ throw ex;
}
int nsph, ies;
int max_ici = 0;
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 2cbafe99a37c1e0db0ec952b387c873745288a26..4c4892f9a79528b43ee409313805f30e97a616f7 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -18,566 +18,583 @@ using namespace std;
* \param output_path: `string` Directory to write the output files in.
*/
void sphere(string config_file, string data_file, string output_path) {
- complex<double> arg, s0, tfsas;
- double th, ph;
- printf("INFO: making legacy configuration...\n");
- ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(config_file);
- GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(data_file);
- if (sconf->number_of_spheres == gconf->number_of_spheres) {
- int isq, ibf;
- double *argi, *args, *gaps;
- 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];
+ complex<double> arg, s0, tfsas;
+ double th, ph;
+ printf("INFO: making legacy configuration...\n");
+ ScattererConfiguration *sconf = NULL;
+ try {
+ sconf = ScattererConfiguration::from_dedfb(config_file);
+ } catch(const OpenConfigurationFileException &ex) {
+ printf("\nERROR: failed to open scatterer configuration file.\n");
+ printf("FILE: %s\n", ex.what());
+ exit(1);
+ }
+ sconf->write_formatted(output_path + "c_OEDFB");
+ sconf->write_binary(output_path + "c_TEDF");
+ GeometryConfiguration *gconf = NULL;
+ try {
+ gconf = GeometryConfiguration::from_legacy(data_file);
+ } catch(const OpenConfigurationFileException &ex) {
+ printf("\nERROR: failed to open geometry configuration file.\n");
+ printf("FILE: %s\n", ex.what());
+ if (sconf) delete sconf;
+ exit(1);
+ }
+ if (sconf->number_of_spheres == gconf->number_of_spheres) {
+ int isq, ibf;
+ double *argi, *args, *gaps;
+ 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];
+ }
+ complex<double> **tqspe, **tqsps;
+ double **tqse, **tqss;
+ 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]();
+ }
+ 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, sconf->nshl_vec, sconf->iog_vec);
+ for (int i = 0; i < nsph; i++) {
+ c1->ros[i] = sconf->radii_of_spheres[i];
+ }
+ 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((output_path + "/c_OSPH").c_str(), "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 = 0; i116 < nsph; i116++) {
+ int i = i116 + 1;
+ int iogi = c1->iog[i116];
+ if (iogi >= i) {
+ double gcss = pi * c1->ros[i116] * c1->ros[i116];
+ c1->gcsv[i116] = gcss;
+ int nsh = c1->nshl[i116];
+ for (int j115 = 0; j115 < nsh; j115++) {
+ c1->rc[i116][j115] = sconf->rcf[i116][j115] * c1->ros[i116];
+ }
+ }
+ gcs += c1->gcsv[iogi];
+ }
+ 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;
+ string tppoan_name = output_path + "/c_TPPOAN_sph";
+ tppoan.open(tppoan_name.c_str(), 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 = 0; jxi488 < sconf->number_of_scales; jxi488++) {
+ printf("INFO: running scale iteration...\n");
+ int jxi = jxi488 + 1;
+ fprintf(output, "========== JXI =%3d ====================\n", jxi);
+ double xi = sconf->scale_vec[jxi488];
+ 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 = 0; i132 < nsph; i132++) {
+ int i = i132 + 1;
+ int iogi = sconf->iog_vec[i132];
+ if (iogi != i) {
+ for (int l123 = 0; l123 < gconf->l_max; l123++) {
+ c1->rmi[l123][i132] = c1->rmi[l123][iogi - 1];
+ c1->rei[l123][i132] = c1->rei[l123][iogi - 1];
+ }
+ continue; // i132
+ }
+ int nsh = sconf->nshl_vec[i132];
+ int ici = (nsh + 1) / 2;
+ if (sconf->idfc == 0) {
+ for (int ic = 0; ic < ici; ic++)
+ c2->dc0[ic] = sconf->dc0_matrix[ic][i132][0]; // WARNING: IDFC=0 is not tested!
+ } else { // IDFC != 0
+ if (jxi == 1) {
+ for (int ic = 0; ic < ici; ic++) {
+ c2->dc0[ic] = sconf->dc0_matrix[ic][i132][jxi488];
+ }
+ }
+ }
+ if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
+ int jer = 0;
+ int lcalc = 0;
+ dme(
+ gconf->l_max, i, 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 jxi == gconf->jwtm) {
+ // This is the condition that writes the transition matrix to output.
+ int is = 1111;
+ fstream ttms;
+ string ttms_name = output_path + "/c_TTMS_sph";
+ ttms.open(ttms_name.c_str(), 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();
+ }
+ }
+ 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);
+ rabas(gconf->in_pol, gconf->l_max, nsph, c1, tqse, tqspe, tqss, tqsps);
+ for (int i170 = 0; i170 < nsph; i170++) {
+ int i = i170 + 1;
+ if (c1->iog[i170] >= i) {
+ double albeds = c1->sscs[i170] / c1->sexs[i170];
+ c1->sqscs[i170] *= sqsfi;
+ c1->sqabs[i170] *= sqsfi;
+ c1->sqexs[i170] *= sqsfi;
+ fprintf(output, " SPHERE %2d\n", i);
+ if (c1->nshl[i170] != 1) {
+ fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i170]);
+ } else {
+ fprintf(
+ output,
+ " SIZE=%15.7lE, REFRACTIVE INDEX=(%15.7lE,%15.7lE)\n",
+ c2->vsz[i170],
+ c2->vkt[i170].real(),
+ c2->vkt[i170].imag()
+ );
+ }
+ fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+ fprintf(
+ output,
+ " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+ c1->sscs[i170], c1->sabs[i170],
+ c1->sexs[i170], albeds
+ );
+ fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
+ fprintf(
+ output,
+ " %14.7lE%15.7lE%15.7lE\n",
+ c1->sqscs[i170], c1->sqabs[i170],
+ c1->sqexs[i170]
+ );
+ fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i170].real(), c1->fsas[i170].imag());
+ double csch = 2.0 * vk * sqsfi / c1->gcsv[i170];
+ s0 = c1->fsas[i170] * exri;
+ double qschu = csch * s0.imag();
+ double pschu = csch * s0.real();
+ double s0mag = cs0 * abs(s0);
+ fprintf(
+ output,
+ " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ qschu, pschu, s0mag
+ );
+ double rapr = c1->sexs[i170] - gaps[i170];
+ double cosav = gaps[i170] / c1->sscs[i170];
+ 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], tqss[0][i170]);
+ fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170], tqss[1][i170]);
+ tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170])), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170])), sizeof(double));
+ double val = tqspe[0][i170].real();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqspe[0][i170].imag();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqsps[0][i170].real();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqsps[0][i170].imag();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&(tqse[1][i170])), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&(tqss[1][i170])), sizeof(double));
+ val = tqspe[1][i170].real();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqspe[1][i170].imag();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqsps[1][i170].real();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqsps[1][i170].imag();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ } // End if iog[i170] >= i
+ } // 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 * abs(s0);
+ fprintf(
+ output,
+ " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ qschu, pschu, s0mag
+ );
+ }
+ th = th1;
+ for (int jth486 = 0; jth486 < nth; jth486++) { // OpenMP parallelizable section
+ int jth = jth486 + 1;
+ ph = ph1;
+ for (int jph484 = 0; jph484 < nph; jph484++) {
+ int jph = jph484 + 1;
+ bool goto182 = (nk == 1) && (jxi > 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 = 0; jths482 < nths; jths482++) {
+ int jths = jths482 + 1;
+ 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 = 0; jphs480 < nphs; jphs480++) {
+ int jphs = jphs480 + 1;
+ if (gconf->meridional_type >= 1) {
+ phs = ph + phsph;
+ if (phs >= 360.0) phs -= 360.0;
}
- complex<double> **tqspe, **tqsps;
- double **tqse, **tqss;
- 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]();
+ bool goto190 = (nks == 1) && ((jxi > 1) || (jth > 1) || (jph > 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);
+ }
}
- 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, sconf->nshl_vec, sconf->iog_vec);
- for (int i = 0; i < nsph; i++) {
- c1->ros[i] = sconf->radii_of_spheres[i];
+ if (nkks != 0 || jxi == 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 = 0; i193 < 3; i193++) {
+ un[i193] = unmp[i193];
+ uns[i193] = unsmp[i193];
+ }
+ }
+ 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));
}
- 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((output_path + "/c_OSPH").c_str(), "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");
+ output,
+ "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n",
+ jth, jph, jths, jphs
+ );
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 = 0; i116 < nsph; i116++) {
- int i = i116 + 1;
- int iogi = c1->iog[i116];
- if (iogi >= i) {
- double gcss = pi * c1->ros[i116] * c1->ros[i116];
- c1->gcsv[i116] = gcss;
- int nsh = c1->nshl[i116];
- for (int j115 = 0; j115 < nsh; j115++) {
- c1->rc[i116][j115] = sconf->rcf[i116][j115] * c1->ros[i116];
- }
- }
- gcs += c1->gcsv[iogi];
- }
- 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;
- string tppoan_name = output_path + "/c_TPPOAN_sph";
- tppoan.open(tppoan_name.c_str(), 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 = 0; jxi488 < sconf->number_of_scales; jxi488++) {
- printf("INFO: running scale iteration...\n");
- int jxi = jxi488 + 1;
- fprintf(output, "========== JXI =%3d ====================\n", jxi);
- double xi = sconf->scale_vec[jxi488];
- 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 = 0; i132 < nsph; i132++) {
- int i = i132 + 1;
- int iogi = sconf->iog_vec[i132];
- if (iogi != i) {
- for (int l123 = 0; l123 < gconf->l_max; l123++) {
- c1->rmi[l123][i132] = c1->rmi[l123][iogi - 1];
- c1->rei[l123][i132] = c1->rei[l123][iogi - 1];
- }
- continue; // i132
- }
- int nsh = sconf->nshl_vec[i132];
- int ici = (nsh + 1) / 2;
- if (sconf->idfc == 0) {
- for (int ic = 0; ic < ici; ic++)
- c2->dc0[ic] = sconf->dc0_matrix[ic][i132][0]; // WARNING: IDFC=0 is not tested!
- } else { // IDFC != 0
- if (jxi == 1) {
- for (int ic = 0; ic < ici; ic++) {
- c2->dc0[ic] = sconf->dc0_matrix[ic][i132][jxi488];
- }
- }
- }
- if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
- int jer = 0;
- int lcalc = 0;
- dme(
- gconf->l_max, i, 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 jxi == gconf->jwtm) {
- // This is the condition that writes the transition matrix to output.
- int is = 1111;
- fstream ttms;
- string ttms_name = output_path + "/c_TTMS_sph";
- ttms.open(ttms_name.c_str(), 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();
- }
- }
- 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);
- rabas(gconf->in_pol, gconf->l_max, nsph, c1, tqse, tqspe, tqss, tqsps);
- for (int i170 = 0; i170 < nsph; i170++) {
- int i = i170 + 1;
- if (c1->iog[i170] >= i) {
- double albeds = c1->sscs[i170] / c1->sexs[i170];
- c1->sqscs[i170] *= sqsfi;
- c1->sqabs[i170] *= sqsfi;
- c1->sqexs[i170] *= sqsfi;
- fprintf(output, " SPHERE %2d\n", i);
- if (c1->nshl[i170] != 1) {
- fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i170]);
- } else {
- fprintf(
- output,
- " SIZE=%15.7lE, REFRACTIVE INDEX=(%15.7lE,%15.7lE)\n",
- c2->vsz[i170],
- c2->vkt[i170].real(),
- c2->vkt[i170].imag()
- );
- }
- fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
- fprintf(
- output,
- " %14.7lE%15.7lE%15.7lE%15.7lE\n",
- c1->sscs[i170], c1->sabs[i170],
- c1->sexs[i170], albeds
- );
- fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
- fprintf(
- output,
- " %14.7lE%15.7lE%15.7lE\n",
- c1->sqscs[i170], c1->sqabs[i170],
- c1->sqexs[i170]
- );
- fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i170].real(), c1->fsas[i170].imag());
- double csch = 2.0 * vk * sqsfi / c1->gcsv[i170];
- s0 = c1->fsas[i170] * exri;
- double qschu = csch * s0.imag();
- double pschu = csch * s0.real();
- double s0mag = cs0 * abs(s0);
- fprintf(
- output,
- " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
- qschu, pschu, s0mag
- );
- double rapr = c1->sexs[i170] - gaps[i170];
- double cosav = gaps[i170] / c1->sscs[i170];
- 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], tqss[0][i170]);
- fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170], tqss[1][i170]);
- tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170])), sizeof(double));
- tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170])), sizeof(double));
- double val = tqspe[0][i170].real();
- tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqspe[0][i170].imag();
- tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[0][i170].real();
- tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[0][i170].imag();
- tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- tppoan.write(reinterpret_cast<char *>(&(tqse[1][i170])), sizeof(double));
- tppoan.write(reinterpret_cast<char *>(&(tqss[1][i170])), sizeof(double));
- val = tqspe[1][i170].real();
- tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqspe[1][i170].imag();
- tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[1][i170].real();
- tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[1][i170].imag();
- tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- } // End if iog[i170] >= i
- } // 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 * abs(s0);
- fprintf(
- output,
- " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
- qschu, pschu, s0mag
- );
- }
- th = th1;
- for (int jth486 = 0; jth486 < nth; jth486++) { // OpenMP parallelizable section
- int jth = jth486 + 1;
- ph = ph1;
- for (int jph484 = 0; jph484 < nph; jph484++) {
- int jph = jph484 + 1;
- bool goto182 = (nk == 1) && (jxi > 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 = 0; jths482 < nths; jths482++) {
- int jths = jths482 + 1;
- 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 = 0; jphs480 < nphs; jphs480++) {
- int jphs = jphs480 + 1;
- if (gconf->meridional_type >= 1) {
- phs = ph + phsph;
- if (phs >= 360.0) phs -= 360.0;
- }
- bool goto190 = (nks == 1) && ((jxi > 1) || (jth > 1) || (jph > 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 || jxi == 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 = 0; i193 < 3; i193++) {
- un[i193] = unmp[i193];
- uns[i193] = unsmp[i193];
- }
- }
- 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",
- jth, jph, jths, jphs
- );
- 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 = 0; ns226 < nsph; ns226++) {
- int ns = ns226 + 1;
- fprintf(output, " SPHERE %2d\n", ns);
- fprintf(
- output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
- c1->sas[ns226][0][0].real(), c1->sas[ns226][0][0].imag(),
- c1->sas[ns226][1][0].real(), c1->sas[ns226][1][0].imag()
- );
- fprintf(
- output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
- c1->sas[ns226][0][1].real(), c1->sas[ns226][0][1].imag(),
- c1->sas[ns226][1][1].real(), c1->sas[ns226][1][1].imag()
- );
- if (jths == 1 && jphs == 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][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
- printf("INFO: done scale.\n");
- } //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;
- for (int zi = gconf->l_max - 1; zi > -1; zi--) {
- for (int zj = 0; zj < 3; zj++) {
- for (int zk = 0; zk < 2; zk++) {
- delete[] zpv[zi][zj][zk];
- }
- delete[] zpv[zi][zj];
- }
- delete[] zpv[zi];
- }
- delete[] zpv;
- 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 = 3; i > -1; i--) {
- delete[] cmul[i];
- delete[] cmullr[i];
- }
- delete[] cmul;
- delete[] cmullr;
- for (int ti = 1; ti > -1; ti--) {
- delete[] tqse[ti];
- delete[] tqss[ti];
- delete[] tqspe[ti];
- delete[] tqsps[ti];
+ 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]);
}
- delete[] tqse;
- delete[] tqss;
- delete[] tqspe;
- delete[] tqsps;
- printf("Finished: output written to %s.\n", (output_path + "/c_OSPH").c_str());
- } else { // NSPH mismatch between geometry and scatterer configurations.
- throw UnrecognizedConfigurationException(
- "Inconsistent geometry and scatterer configurations."
- );
+ sscr2(nsph, gconf->l_max, vk, exri, c1);
+ for (int ns226 = 0; ns226 < nsph; ns226++) {
+ int ns = ns226 + 1;
+ fprintf(output, " SPHERE %2d\n", ns);
+ fprintf(
+ output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
+ c1->sas[ns226][0][0].real(), c1->sas[ns226][0][0].imag(),
+ c1->sas[ns226][1][0].real(), c1->sas[ns226][1][0].imag()
+ );
+ fprintf(
+ output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
+ c1->sas[ns226][0][1].real(), c1->sas[ns226][0][1].imag(),
+ c1->sas[ns226][1][1].real(), c1->sas[ns226][1][1].imag()
+ );
+ if (jths == 1 && jphs == 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][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
+ printf("INFO: done scale.\n");
+ } //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;
+ for (int zi = gconf->l_max - 1; zi > -1; zi--) {
+ for (int zj = 0; zj < 3; zj++) {
+ for (int zk = 0; zk < 2; zk++) {
+ delete[] zpv[zi][zj][zk];
}
- delete sconf;
- delete gconf;
+ delete[] zpv[zi][zj];
+ }
+ delete[] zpv[zi];
+ }
+ delete[] zpv;
+ 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 = 3; i > -1; i--) {
+ delete[] cmul[i];
+ delete[] cmullr[i];
+ }
+ delete[] cmul;
+ delete[] cmullr;
+ for (int ti = 1; ti > -1; ti--) {
+ delete[] tqse[ti];
+ delete[] tqss[ti];
+ delete[] tqspe[ti];
+ delete[] tqsps[ti];
+ }
+ delete[] tqse;
+ delete[] tqss;
+ delete[] tqspe;
+ delete[] tqsps;
+ printf("Finished: output written to %s.\n", (output_path + "/c_OSPH").c_str());
+ } else { // NSPH mismatch between geometry and scatterer configurations.
+ throw UnrecognizedConfigurationException(
+ "Inconsistent geometry and scatterer configurations."
+ );
+ }
+ delete sconf;
+ delete gconf;
}