diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 0517eb3406ac1118e56dcac50ce6c044d1930b6b..2e8242e27d8d8d75520fead5f939d5c103031c1e 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -174,7 +174,7 @@ sanity_stage:
- building_stage
artifacts:
paths:
- - build/testing/valgrind.log
+ - build/testing/valgrind_*.log
exclude:
- ".git*"
- ".git/**/*"
@@ -190,14 +190,16 @@ sanity_stage:
- cd ../testing
- echo "Running memory sanity check for ParticleDescriptor"
- chmod +x test_ParticleDescriptor
- - valgrind --leak-check=full --log-file=valgrind.log ./test_ParticleDescriptor
- - grep "0 errors from 0 contexts" valgrind.log
+ - valgrind --leak-check=full --log-file=valgrind_Particle.log ./test_ParticleDescriptor
+ - grep "0 errors from 0 contexts" valgrind_Particle.log
- echo "Running memory sanity check for output classes"
- - chmod +x test_outputs
- - rm valgrind.log
- - valgrind --leak-check=full --log-file=valgrind.log ./test_outputs
- - grep "0 errors from 0 contexts" valgrind.log
+ - chmod +x test_cluster_outputs
+ - valgrind --leak-check=full --log-file=valgrind_cluster.log ./test_cluster_outputs
+ - grep "0 errors from 0 contexts" valgrind_cluster.log
- rm -rf c_OCLU_24
+ - chmod +x test_inclusion_outputs
+ - valgrind --leak-check=full --log-file=valgrind_inclusion.log ./test_inclusion_outputs
+ - grep "0 errors from 0 contexts" valgrind_inclusion.log
running_stage:
stage: run
@@ -294,10 +296,16 @@ testing_stage:
- export FFILE=../../test_data/cluster/OCLU_48
- python3 ../../src/scripts/pycompare.py --no-progress --ffile $FFILE --cfile c_OCLU --html=pycompare_48.html
- cd ../testing
- - echo "Checking consistency of HDF5 output"
- - chmod u+x test_outputs
- - ./test_outputs
+ - echo "Checking consistency of HDF5 cluster output"
+ - chmod u+x test_cluster_outputs
+ - ./test_cluster_outputs
- export FFILE=../../test_data/cluster/OCLU_24
- python3 ../../src/scripts/pycompare.py --no-progress --ffile $FFILE --cfile c_OCLU_24
- rm -rf c_OCLU_24
+ - echo "Checking consistency of HDF5 incluson output"
+ - chmod u+x test_inclusion_outputs
+ - ./test_inclusion_outputs
+ - export FFILE=../../test_data/inclusion/OINCLU
+ - python3 ../../src/scripts/pycompare.py --no-progress --ffile $FFILE --cfile c_OINCLU
+ - rm -rf c_OINCLU
\ No newline at end of file
diff --git a/build/Makefile b/build/Makefile
index 10f9d355d7f4072decd2e46c4e9ef78865cd4827..cba30075e68395cbbb87ea589cb98b5d691a00da 100644
--- a/build/Makefile
+++ b/build/Makefile
@@ -29,8 +29,8 @@ NP_SPHERE_BINS=sphere/np_sphere
NP_TRAPPING_SRCS=../src/trapping/np_trapping.cpp ../src/trapping/cfrfme.cpp ../src/trapping/clffft.cpp
NP_TRAPPING_OBJS=../src/trapping/np_trapping.o ../src/trapping/cfrfme.o ../src/trapping/clffft.o
NP_TRAPPING_BINS=trapping/np_trapping
-NP_TESTING_OBJS=../src/testing/test_outputs.o ../src/testing/test_ParticleDescriptor.o ../src/testing/test_TEDF.o ../src/testing/test_TTMS.o
-NP_TESTING_BINS=testing/test_outputs testing/test_ParticleDescriptor testing/test_TEDF testing/test_TTMS
+NP_TESTING_OBJS=../src/testing/test_cluster_outputs.o ../src/testing/test_inclusion_outputs.o ../src/testing/test_ParticleDescriptor.o ../src/testing/test_TEDF.o ../src/testing/test_TTMS.o
+NP_TESTING_BINS=testing/test_cluster_outputs testing/test_inclusion_outputs testing/test_ParticleDescriptor testing/test_TEDF testing/test_TTMS
all: $(NPTM_LIB) $(FORTRAN_BINS) $(NP_CLUSTER_BINS) $(NP_INCLUSION_BINS) $(NP_SPHERE_BINS) $(NP_TRAPPING_BINS) $(NP_TESTING_BINS)
@@ -83,8 +83,11 @@ trapping/np_trapping: $(NPTM_LIB) $(NP_TRAPPING_OBJS)
testing/test_ParticleDescriptor: $(NPTM_LIB) ../src/testing/test_ParticleDescriptor.o
$(CXX) $(CXXFLAGS) ../src/testing/test_ParticleDescriptor.o -o $@ $(CXXLDFLAGS)
-testing/test_outputs: $(NPTM_LIB) ../src/testing/test_outputs.o
- $(CXX) $(CXXFLAGS) ../src/testing/test_outputs.o -o $@ $(CXXLDFLAGS)
+testing/test_cluster_outputs: $(NPTM_LIB) ../src/testing/test_cluster_outputs.o
+ $(CXX) $(CXXFLAGS) ../src/testing/test_cluster_outputs.o -o $@ $(CXXLDFLAGS)
+
+testing/test_inclusion_outputs: $(NPTM_LIB) ../src/testing/test_inclusion_outputs.o
+ $(CXX) $(CXXFLAGS) ../src/testing/test_inclusion_outputs.o -o $@ $(CXXLDFLAGS)
testing/test_TEDF: $(NPTM_LIB) ../src/testing/test_TEDF.o
$(CXX) $(CXXFLAGS) ../src/testing/test_TEDF.o -o $@ $(CXXLDFLAGS)
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index ae4fae9fed32536d117bd990fb9c80ca30ac5c83..7d2533eac1e1f313fcd14b29e8805c0aeb229534 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -112,7 +112,7 @@ using namespace std;
* \param gconf: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` object.
* \param sa: `ScatteringAngles *` Pointer to a `ScatteringAngles` object.
* \param cid: `ClusterIterationData *` Pointer to a `ClusterIterationData` object.
- * \param output: `ClusterOutputInfo *` Pointer to a `ClusterOutputInfo` object.
+ * \param oi: `ClusterOutputInfo *` Pointer to a `ClusterOutputInfo` object.
* \param output_path: `const string &` Path to the output directory.
* \param vtppoanp: `VirtualBinaryFile *` Pointer to a `VirtualBinaryFile` object.
*/
@@ -326,6 +326,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
return;
}
+ //==================================================
+ // do the first outputs here, so that I open here the new files, afterwards I only append
+ //==================================================
vtppoanp->write_to_disk(output_path + "/c_TPPOAN");
delete vtppoanp;
@@ -420,8 +423,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
// the parallel loop over MPI processes covers a different set of indices for each thread
#pragma omp barrier
int myjxi488 = ixi488+myompthread;
- vtppoanp_2 = new VirtualBinaryFile();
// each thread opens new virtual files and stores their pointers in the shared array
+ vtppoanp_2 = new VirtualBinaryFile();
// each thread puts a copy of the pointers to its virtual files in the shared arrays
vtppoanarray[myompthread] = vtppoanp_2;
#pragma omp barrier
@@ -521,8 +524,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException(
- "Inconsistent geometry and scatterer configurations."
- );
+ "Inconsistent geometry and scatterer configurations."
+ );
}
delete sconf;
@@ -1150,7 +1153,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
output->vec_qschuc2[jindex - 1] = qschum;
output->vec_pschuc2[jindex - 1] = pschum;
output->vec_s0magc2[jindex - 1] = s0magm;
- }
+ } // ipol polarization switch
if (ilr210 == 1 && jlr == 2) {
output->vec_fsac11[jindex - 1] = cid->c1->fsacm[0][0];
output->vec_fsac21[jindex - 1] = cid->c1->fsacm[1][0];
@@ -1489,14 +1492,14 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
}
// label 318
for (int i = 0; i < 4; i++) {
- oindex = 16 * (jindex - 1) + 4 * i;
+ oindex = 16 * (jindex - 1) + 4 * i; // if IAVM fails, try adding directions
output->vec_dir_mulc[oindex] = cid->cmul[i][0];
output->vec_dir_mulc[oindex + 1] = cid->cmul[i][1];
output->vec_dir_mulc[oindex + 2] = cid->cmul[i][2];
output->vec_dir_mulc[oindex + 3] = cid->cmul[i][3];
}
for (int i = 0; i < 4; i++) {
- oindex = 16 * (jindex - 1) + 4 * i;
+ oindex = 16 * (jindex - 1) + 4 * i; // if IAVM fails, try adding directions
output->vec_dir_mulclr[oindex] = cid->cmullr[i][0];
output->vec_dir_mulclr[oindex + 1] = cid->cmullr[i][1];
output->vec_dir_mulclr[oindex + 2] = cid->cmullr[i][2];
diff --git a/src/include/algebraic.h b/src/include/algebraic.h
index d19db95674e0225e5fb0c4dc259384b3f1db2005..c59dbae2c43a452313970d1e3cca215b407cc55a 100644
--- a/src/include/algebraic.h
+++ b/src/include/algebraic.h
@@ -39,8 +39,9 @@ using namespace std;
* \param maxrefiters: `int &` Reference to the maximum number of refinement iterations.
* \param accuracygoal: `double &` Reference to the requested accuracy level.
* \param refinemode: `int` Flag for refinement mode selection.
- * \param max_size: `np_int` The maximum expected size (required by some call-backs,
- * optional, defaults to 0).
+ * \param output_path: `const string &` Path where the output needs to be placed.
+ * \param jxi488: `int` Index of the current wavelength calculation.
+ * \param max_size: `np_int` The maximum expected size (required by some call-backs, optional, defaults to 0).
* \param target_device: `int` ID of target GPU, if available (defaults to 0).
*/
void invert_matrix(dcomplex **mat, np_int size, int &ier, int &maxrefiters, double &accuracygoal, int refinemode, const string& output_path, int jxi488, np_int max_size=0, int target_device=0);
diff --git a/src/include/magma_calls.h b/src/include/magma_calls.h
index 32052810b5ea125ad40648f6e5747fa19643fa23..679b07f37c70b0cc88036f0af4e4cca9163ada78 100644
--- a/src/include/magma_calls.h
+++ b/src/include/magma_calls.h
@@ -59,20 +59,27 @@ void magma_zinvert1(dcomplex * &inva, np_int n, int &jer, int device_id);
* \param accuracygoal: `double &` Accuracy to achieve in iterative refinement, defined as the module of the maximum difference between the identity matrix and the matrix product of the (approximate) inverse times the original matrix. On return, it contains the actually achieved accuracy.
* \param refinemode: `int` Flag to control the refinement mode.
* \param device_id: `int` ID of the device for matrix inversion offloading.
+ * \param output_path: `const string &` Path where the output needs to be placed.
+ * \param jxi488: `int` Index of the current wavelength calculation.
*/
void magma_zinvert_and_refine(dcomplex **mat, np_int n, int &jer, int &maxrefiters, double &accuracygoal, int refinemode, int device_id, const string& output_path, int jxi488);
-/*! \brief apply iterative refinement of the solution of a matrix inversion
+/*! \brief Apply iterative refinement of the solution of a matrix inversion.
*
- * iteratively compute and apply a correction to the inverse inva of the complex matrix aorig, for a maximum number of maxiters times, or until achieving a maximum residual better than accuracygoal
+ * Iteratively compute and apply a correction to the inverse `inva` of the complex
+ * matrix `aorig`, for a maximum number of `maxiters` times, or until achieving a
+ * maximum residual better than `accuracygoal`.
*
* \param aorig: pointer to the first element of the matrix of complex to be inverted.
* \param inva: pointer to the first element of inverse.
* \param n: `np_int` The number of rows and columns of the [n x n] matrices.
* \param jer: `int &` Reference to an integer return flag.
* \param maxrefiters: `int` Maximum number of refinement iterations to apply.
- * \param accuracygoal: `double` Accuracy to achieve in iterative refinement, defined as the module of the maximum difference between the identity matrix and the matrix product of the (approximate) inverse times the original matrix. On return, it contains the actually achieved accuracy
+ * \param accuracygoal: `double` Accuracy to achieve in iterative refinement, defined as the module of the maximum difference between the identity matrix and the matrix product of the (approximate) inverse times the original matrix. On return, it contains the actually achieved accuracy.
+ * \param refinemode: `int` Flag for refinement mode selection.
* \param device_id: `int` ID of the device for matrix inversion offloading.
+ * \param output_path: `const string &` Path where the output needs to be placed.
+ * \param jxi488: `int` Index of the current wavelength calculation.
*/
void magma_refine(dcomplex *aorig, dcomplex *inva, np_int n, int &jer, int &maxrefiters, double &accuracygoal, int refinemode, int device_id, const string& output_path, int jxi488);
diff --git a/src/include/outputs.h b/src/include/outputs.h
index 1370745b2074769193ab76d39bda0d30ba40758b..b508e79d02e6510fcbfac6b9fc3d96e8b82a5a6e 100644
--- a/src/include/outputs.h
+++ b/src/include/outputs.h
@@ -21,6 +21,7 @@
#ifndef INCLUDE_OUTPUTS_H_
#define INCLUDE_OUTPUTS_H_
+// >>> OUTPUT FOR CLUSTER <<<
/*! \brief Class to collect output information for scattering from clusters.
*
* The results of the calculation can be saved in different formats.
@@ -54,12 +55,7 @@ protected:
*
* This function takes care of writing the output using the legacy
* formatted ASCII structure. If the output file does not exist, it
- * is created. If it exists, the new content is appended to the old
- * one, without overwriting. This is useful for a simpler handling
- * of multi-process output combination, provided that only one process
- * writes to the output file at any given time. Simultaneous write
- * operations by different processes is not supported by the legacy
- * ASCII format.
+ * is created. If it exists, the new content overwritten.
*
* \param output: `const string &` Path to the output to be written.
* \return result: `int` Exit code (0 if successful).
@@ -500,7 +496,7 @@ public:
/*! \brief Insert in the current output data the data of another block.
*
- * \param rhs: `const ClusterOutputData &` Reference to the source data block.
+ * \param rhs: `const ClusterOutputInfo &` Reference to the source data block.
* \return result: `int` Exit code (0 if successful).
*/
int insert(const ClusterOutputInfo &rhs);
@@ -541,4 +537,436 @@ public:
int mpisend(const mixMPI *mpidata);
#endif // MPI_VERSION
};
+// >>> END OF OUTPUT FOR CLUSTER <<<
+
+// >>> OUTPUT FOR INCLUSION <<<
+/*! \brief Class to collect output information for scattering from particle with inclusions.
+ *
+ * The results of the calculation can be saved in different formats.
+ * It is therefore convenient to have a proper memory structure that
+ * allows for storing the results and flushing them in any of the
+ * permitted formats with just one operation. The purpose of the
+ * `InclusionOutputInfo` class is to provide a wrapper for the output
+ * of the particle with inclusions scattering solver.
+ */
+class InclusionOutputInfo {
+protected:
+ //! \brief Number of incident azimuth calculations.
+ int _num_theta;
+ //! \brief Number of scattered azimuth calculations.
+ int _num_thetas;
+ //! \brief Number of incident elevation calculations.
+ int _num_phi;
+ //! \brief Number of scattered elevation calculations.
+ int _num_phis;
+ //! \brief ID of the first computed wavelength
+ int _first_xi;
+
+ /*! \brief Write the output to a HDF5 file.
+ *
+ * \param file_name: `const string &` Path to the output to be written.
+ * \return result: `int` Exit code (0 if successful).
+ */
+ int write_hdf5(const std::string &file_name);
+
+ /*! \brief Write the output to a legacy text file.
+ *
+ * This function takes care of writing the output using the legacy
+ * formatted ASCII structure. If the output file does not exist, it
+ * is created. If it exists, the new content is overwritten.
+ *
+ * \param output: `const string &` Path to the output to be written.
+ * \return result: `int` Exit code (0 if successful).
+ */
+ int write_legacy(const std::string &output);
+
+public:
+ //! \brief Read-only view on the ID of the first scale
+ const int &first_xi = _first_xi;
+ //! \brief Number of spheres in the aggregate.
+ int nsph;
+ //! \brief Maximum internal field expansion order.
+ int li;
+ //! \brief Maximum external field expansion order.
+ int le;
+ //! \brief Maximum field expansion order.
+ int lm;
+ //! \brief Maximum coefficient matrix dimension.
+ np_int mxndm;
+ //! \brief Incident polarization flag.
+ int inpol;
+ //! \brief Number of points for transition layer integration.
+ int npnt;
+ //! \brief Number of points for non-transition layer integration.
+ int npntts;
+ //! \brief Flag for intensity.
+ int iavm;
+ //! \brief Flag for reference to meridional plane.
+ int isam;
+ //! \brief Flag for dielectric function definition.
+ int idfc;
+ //! \brief Vector of spherical components X Cartesian coordinates.
+ double *vec_x_coords;
+ //! \brief Vector of spherical components Y Cartesian coordinates.
+ double *vec_y_coords;
+ //! \brief Vector of spherical components Z Cartesian coordinates.
+ double *vec_z_coords;
+ //! \brief First incident radiation azimuth angle.
+ double th;
+ //! \brief Incident radiation azimuth angle step.
+ double thstp;
+ //! \brief Last incident radiation azimuth angle.
+ double thlst;
+ //! \brief First scattered radiation azimuth angle.
+ double ths;
+ //! \brief Scattered radiation azimuth angle step.
+ double thsstp;
+ //! \brief Last scattered radiation azimuth angle.
+ double thslst;
+ //! \brief First incident radiation elevation angle.
+ double ph;
+ //! \brief Incident radiation elevation angle step.
+ double phstp;
+ //! \brief Last incident radiation elevation angle.
+ double phlst;
+ //! \brief First scattered radiation elevation angle.
+ double phs;
+ //! \brief Scattered radiation elevation angle step.
+ double phsstp;
+ //! \brief Last scattered radiation elevation angle.
+ double phslst;
+ //! \brief Number of directions to be explicitly solved.
+ int ndirs;
+ //! \brief Refractive index of external medium.
+ double exri;
+ //! \brief Number of scales (wavelengths)
+ int nxi;
+ //! \brief Number of scales handled by the current process.
+ int xi_block_size;
+ //! \brief Index of the wavelength for T-matrix output.
+ int jwtm;
+ //! \brief Vector of scale (wavelength) indices.
+ int *vec_jxi;
+ //! \brief Vector of error severities (0 - success, 1 - HJV, 2 - DME).
+ short *vec_ier;
+ //! \brief Vector of vacuum wave numbers.
+ double *vec_vk;
+ //! \brief Vector of computed scales.
+ double *vec_xi;
+ //! \brief Number of sphere configurations.
+ int configurations;
+ //! \brief Vector of sphere sizes (all configurations for every scale).
+ double *vec_sphere_sizes;
+ //! \brief Vector of sphere refractive indices (all configurations for every scale).
+ dcomplex *vec_sphere_ref_indices;
+ //! \brief Vector of particle scattering cross-sections (parallel polarization).
+ double *vec_scs1;
+ //! \brief Vector of particle scattering cross-sections (perpendicular polarization).
+ double *vec_scs2;
+ //! \brief Vector of particle absorption cross-sections (parallel polarization).
+ double *vec_abs1;
+ //! \brief Vector of particle absorption cross-sections (perpendicular polarization).
+ double *vec_abs2;
+ //! \brief Vector of particle extinction cross-sections (parallel polarization).
+ double *vec_exs1;
+ //! \brief Vector of particle extinction cross-sections (perpendicular polarization).
+ double *vec_exs2;
+ //! \brief Vector of particle albedos (parallel polarization).
+ double *vec_albeds1;
+ //! \brief Vector of particle albedos (perpendicular polarization).
+ double *vec_albeds2;
+ //! \brief Vector of particle scattering-to-geometric cross-sections (parallel polarization).
+ double *vec_scsrt1;
+ //! \brief Vector of particle scattering-to-geometric cross-sections (perpendicular polarization).
+ double *vec_scsrt2;
+ //! \brief Vector of particle absorption-to-geometric cross-sections (parallel polarization).
+ double *vec_absrt1;
+ //! \brief Vector of particle absorption-to-geometric cross-sections (perpendicular polarization).
+ double *vec_absrt2;
+ //! \brief Vector of particle extinction-to-geometric cross-sections (parallel polarization).
+ double *vec_exsrt1;
+ //! \brief Vector of particle extinction-to-geometric cross-sections (perpendicular polarization).
+ double *vec_exsrt2;
+ //! \brief Vector of particle QSCHU (parallel polarization).
+ double *vec_qschu1;
+ //! \brief Vector of particle QSCHU (perpendicular polarization).
+ double *vec_qschu2;
+ //! \brief Vector of particle PSCHU (parallel polarization).
+ double *vec_pschu1;
+ //! \brief Vector of particle PSCHU (perpendicular polarization).
+ double *vec_pschu2;
+ //! \brief Vector of particle S0MAG (parallel polarization).
+ double *vec_s0mag1;
+ //! \brief Vector of particle S0MAG (perpendicular polarization).
+ double *vec_s0mag2;
+ //! \brief Vector of particle average asymmetry parameter (parallel polarization).
+ double *vec_cosav1;
+ //! \brief Vector of particle average asymmetry parameter (perpendicular polarization).
+ double *vec_cosav2;
+ //! \brief Vector of particle average radiation pressure force (N - parallel polarization).
+ double *vec_raprs1;
+ //! \brief Vector of particle average radiation pressure force (N - perpendicular polarization).
+ double *vec_raprs2;
+ //! \brief Vector of particle average radiation force along incidence direction (N - parallel polarization).
+ double *vec_fk1;
+ //! \brief Vector of particle average radiation force along incidence direction (N - perpendicular polarization).
+ double *vec_fk2;
+ //! \brief Vector of forward scattering amplitudes for polarization parallel to incidence (one per scale).
+ dcomplex *vec_fsas11;
+ //! \brief Vector of forward scattering amplitudes for polarization perpendicular to incidence (one per scale).
+ dcomplex *vec_fsas21;
+ //! \brief Vector of forward scattering amplitudes for polarization parallel to incidence (one per scale).
+ dcomplex *vec_fsas22;
+ //! \brief Vector of forward scattering amplitudes for polarization perpendicular to incidence (one per scale).
+ dcomplex *vec_fsas12;
+ //! \brief Vector of incidence azimuth directions (one per incidence azimuth).
+ double *vec_dir_tidg;
+ //! \brief Vector of incidence elevation directions (one per incidence elevation).
+ double *vec_dir_pidg;
+ //! \brief Vector of scattering azimuth directions (one per scattering azimuth).
+ double *vec_dir_tsdg;
+ //! \brief Vector of scattering elevation directions (one per scattering elevation).
+ double *vec_dir_psdg;
+ //! \brief Vector of scattering angles (one per direction).
+ double *vec_dir_scand;
+ //! \brief Control parameter for incidence plane referred to meridional plane (one per direction).
+ double *vec_dir_cfmp;
+ //! \brief Control parameter for scattering plane referred to meridional plane (one per direction).
+ double *vec_dir_sfmp;
+ //! \brief Control parameter for incidence plane referred to scattering plane (one per direction).
+ double *vec_dir_cfsp;
+ //! \brief Control parameter for scattering plane referred to scattering plane (one per direction).
+ double *vec_dir_sfsp;
+ //! \brief Components of the unitary vector perpendicular to incidence plane (three per direction).
+ double *vec_dir_un;
+ //! \brief Components of the unitary vector perpendicular to scattering plane (three per direction).
+ double *vec_dir_uns;
+ //! \brief Vector of particle differential scattering cross-sections (parallel polarization).
+ double *vec_dir_scs1;
+ //! \brief Vector of particle differential scattering cross-sections (perpendicular polarization).
+ double *vec_dir_scs2;
+ //! \brief Vector of particle differential absorption cross-sections (parallel polarization).
+ double *vec_dir_abs1;
+ //! \brief Vector of particle differential absorption cross-sections (perpendicular polarization).
+ double *vec_dir_abs2;
+ //! \brief Vector of particle differential extinction cross-sections (parallel polarization).
+ double *vec_dir_exs1;
+ //! \brief Vector of particle differential extinction cross-sections (perpendicular polarization).
+ double *vec_dir_exs2;
+ //! \brief Vector of particle differential albedos (parallel polarization).
+ double *vec_dir_albeds1;
+ //! \brief Vector of particle differential albedos (perpendicular polarization).
+ double *vec_dir_albeds2;
+ //! \brief Vector of particle differential scattering-to-geometric cross-sections (parallel polarization).
+ double *vec_dir_scsrt1;
+ //! \brief Vector of particle differential scattering-to-geometric cross-sections (perpendicular polarization).
+ double *vec_dir_scsrt2;
+ //! \brief Vector of particle differential absorption-to-geometric cross-sections (parallel polarization).
+ double *vec_dir_absrt1;
+ //! \brief Vector of particle differential absorption-to-geometric cross-sections (perpendicular polarization).
+ double *vec_dir_absrt2;
+ //! \brief Vector of particle differential extinction-to-geometric cross-sections (parallel polarization).
+ double *vec_dir_exsrt1;
+ //! \brief Vector of particle differential extinction-to-geometric cross-sections (perpendicular polarization).
+ double *vec_dir_exsrt2;
+ //! \brief Vector of particle differential forward scattering amplitude with polarization parallel to parallel incidence field.
+ dcomplex *vec_dir_fsas11;
+ //! \brief Vector of particle differential forward scattering amplitude with polarization perpendicular to the parallel incidence field.
+ dcomplex *vec_dir_fsas21;
+ //! \brief Vector of particle differential forward scattering amplitude with polarization perpendicular to perpendicular incidence field.
+ dcomplex *vec_dir_fsas12;
+ //! \brief Vector of particle differential forward scattering amplitude with polarization parallel the perpendicular incidence field.
+ dcomplex *vec_dir_fsas22;
+ //! \brief Vector of particle differential scattering amplitude with polarization parallel to parallel incidence field.
+ dcomplex *vec_dir_sas11;
+ //! \brief Vector of particle differential scattering amplitude with polarization perpendicular to the parallel incidence field.
+ dcomplex *vec_dir_sas21;
+ //! \brief Vector of particle differential scattering amplitude with polarization perpendicular to perpendicular incidence field.
+ dcomplex *vec_dir_sas12;
+ //! \brief Vector of particle differential scattering amplitude with polarization parallel the perpendicular incidence field.
+ dcomplex *vec_dir_sas22;
+ //! \brief Vector of differential particle QSCHU (parallel polarization).
+ double *vec_dir_qschu1;
+ //! \brief Vector of differential particle QSCHU (perpendicular polarization).
+ double *vec_dir_qschu2;
+ //! \brief Vector of differential particle PSCHU (parallel polarization).
+ double *vec_dir_pschu1;
+ //! \brief Vector of differential particle PSCHU (perpendicular polarization).
+ double *vec_dir_pschu2;
+ //! \brief Vector of particle differential S0MAG (parallel polarization).
+ double *vec_dir_s0mag1;
+ //! \brief Vector of particle differential S0MAG (perpendicular polarization).
+ double *vec_dir_s0mag2;
+ //! \brief Vector of differential particle asymmetry parameters (parallel polarization).
+ double *vec_dir_cosav1;
+ //! \brief Vector of differential particle asymmetry parameters (perpendicular polarization).
+ double *vec_dir_cosav2;
+ //! \brief Vector of differential particle radiation pressure forces (1).
+ double *vec_dir_rapr1;
+ //! \brief Vector of differential particle radiation pressure forces (1).
+ double *vec_dir_rapr2;
+ //! \brief Vector of differential radiation pressure force components along the polarization direction (parallel polarization).
+ double *vec_dir_fl1;
+ //! \brief Vector of differential radiation pressure force components along the polarization direction (perpendicular polarization).
+ double *vec_dir_fl2;
+ //! \brief Vector of differential radiation pressure force components perpendicular to the polarization direction (parallel polarization).
+ double *vec_dir_fr1;
+ //! \brief Vector of differential radiation pressure force components perpendicular to the polarization direction (perpendicular polarization).
+ double *vec_dir_fr2;
+ //! \brief Vector of differential radiation pressure force components along the incidence direction (parallel polarization).
+ double *vec_dir_fk1;
+ //! \brief Vector of differential radiation pressure force components along the incidence direction (perpendicular polarization).
+ double *vec_dir_fk2;
+ //! \brief Vector of differential radiation pressure force components along the X axis (parallel polarization).
+ double *vec_dir_fx1;
+ //! \brief Vector of differential radiation pressure force components along the X axis (perpendicular polarization).
+ double *vec_dir_fx2;
+ //! \brief Vector of differential radiation pressure force components along the Y axis (parallel polarization).
+ double *vec_dir_fy1;
+ //! \brief Vector of differential radiation pressure force components along the Y axis (perpendicular polarization).
+ double *vec_dir_fy2;
+ //! \brief Vector of differential radiation pressure force components along the Z axis (parallel polarization).
+ double *vec_dir_fz1;
+ //! \brief Vector of differential radiation pressure force components along the Z axis (perpendicular polarization).
+ double *vec_dir_fz2;
+ //! \brief Vector of differential extinction contribution to radiation torque components along the polarization direction (parallel polarization).
+ double *vec_dir_tqel1;
+ //! \brief Vector of differential extinction contribution to radiation torque components along the polarization direction (perpendicular polarization).
+ double *vec_dir_tqel2;
+ //! \brief Vector of differential extinction contribution to radiation torque components perpendicular to the polarization direction (parallel polarization).
+ double *vec_dir_tqer1;
+ //! \brief Vector of differential extinction contribution to radiation torque components perpendicular to the polarization direction (perpendicular polarization).
+ double *vec_dir_tqer2;
+ //! \brief Vector of differential extinction contribution to radiation torque components along the incidence direction (parallel polarization).
+ double *vec_dir_tqek1;
+ //! \brief Vector of differential extinction contribution to radiation torque components along the incidence direction (perpendicular polarization).
+ double *vec_dir_tqek2;
+ //! \brief Vector of differential extinction contribution to radiation torque components along the X axis (parallel polarization).
+ double *vec_dir_tqex1;
+ //! \brief Vector of differential extinction contribution to radiation torque components along the X axis (perpendicular polarization).
+ double *vec_dir_tqex2;
+ //! \brief Vector of differential extinction contribution to radiation torque components along the Y axis (parallel polarization).
+ double *vec_dir_tqey1;
+ //! \brief Vector of differential extinction contribution to radiation torque components along the Y axis (perpendicular polarization).
+ double *vec_dir_tqey2;
+ //! \brief Vector of differential extinction contribution to radiation torque components along the Z axis (parallel polarization).
+ double *vec_dir_tqez1;
+ //! \brief Vector of differential extinction contribution to radiation torque components along the Z axis (perpendicular polarization).
+ double *vec_dir_tqez2;
+ //! \brief Vector of differential scattering contribution to radiation torque components along the polarization direction (parallel polarization).
+ double *vec_dir_tqsl1;
+ //! \brief Vector of differential scattering contribution to radiation torque components along the polarization direction (perpendicular polarization).
+ double *vec_dir_tqsl2;
+ //! \brief Vector of differential scattering contribution to radiation torque components perpendicular to the polarization direction (parallel polarization).
+ double *vec_dir_tqsr1;
+ //! \brief Vector of differential scattering contribution to radiation torque components perpendicular to the polarization direction (perpendicular polarization).
+ double *vec_dir_tqsr2;
+ //! \brief Vector of differential scattering contribution to radiation torque components along the incidence direction (parallel polarization).
+ double *vec_dir_tqsk1;
+ //! \brief Vector of differential scattering contribution to radiation torque components along the incidence direction (perpendicular polarization).
+ double *vec_dir_tqsk2;
+ //! \brief Vector of differential scattering contribution to radiation torque components along X axis (parallel polarization).
+ double *vec_dir_tqsx1;
+ //! \brief Vector of differential scattering contribution to radiation torque components along X axis (perpendicular polarization).
+ double *vec_dir_tqsx2;
+ //! \brief Vector of differential scattering contribution to radiation torque components along Y axis (parallel polarization).
+ double *vec_dir_tqsy1;
+ //! \brief Vector of differential scattering contribution to radiation torque components along Y axis (perpendicular polarization).
+ double *vec_dir_tqsy2;
+ //! \brief Vector of differential scattering contribution to radiation torque components along Z axis (parallel polarization).
+ double *vec_dir_tqsz1;
+ //! \brief Vector of differential scattering contribution to radiation torque components along Z axis (perpendicular polarization).
+ double *vec_dir_tqsz2;
+ //! \brief Vector of cluster Mueller transormation matrices referred to meridional plane (16 per direction per scale).
+ double *vec_dir_mull;
+ //! \brief Vector of cluster Mueller transormation matrices referred to scattering plane (16 per direction per scale).
+ double *vec_dir_mulllr;
+
+ /*! \brief `InclusionOutputInfo` default instance constructor.
+ *
+ * \param sc: `ScattererConfiguration *` Pointer to a `ScattererConfiguration` instance.
+ * \param gc: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` instance.
+ * \param mpidata: `const mixMPI*` Pointer to a mixMPI instance.
+ * \param first_xi: `int` Index of the first scale in output (optional, default is 1).
+ * \param xi_length: `int` Number of scales tobe included in output (optional, default is 0, meaning all).
+ */
+ InclusionOutputInfo(
+ ScattererConfiguration *sc, GeometryConfiguration *gc,
+ const mixMPI *mpidata, int first_xi = 1, int xi_length = 0
+ );
+
+ /*! \brief `InclusionOutputInfo` constructor from HDF5 input.
+ *
+ * \param hdf5_file_name: `const string &` Path to the HDF5 file to be read.
+ */
+ InclusionOutputInfo(const std::string &hdf5_file_name);
+
+ /*! \brief `InclusionOutputInfo` instance destroyer.
+ */
+ ~InclusionOutputInfo();
+
+ /*! \brief Estimate the size of the structure that would be built for given input.
+ *
+ * \param sc: `ScattererConfiguration *` Pointer to a `ScattererConfiguration` instance.
+ * \param gc: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` instance.
+ * \param first_xi: `int` Index of the first scale in output (optional, default is 0).
+ * \param xi_length: `int` Number of scales tobe included in output (optional, default is all).
+ * \return size: `long` Estimated instance size in bytes.
+ */
+ static long compute_size(
+ ScattererConfiguration *sc, GeometryConfiguration *gc,
+ int first_xi = 0, int xi_length = 0
+ );
+
+ /*! \brief Get the size of a `ClusterOutputInfo` instance in bytes.
+ *
+ * \return size: `long` Estimated instance size in bytes.
+ */
+ long compute_size();
+
+ /*! \brief Insert in the current output data the data of another block.
+ *
+ * \param rhs: `const InclusionOutputInfo &` Reference to the source data block.
+ * \return result: `int` Exit code (0 if successful).
+ */
+ int insert(const InclusionOutputInfo &rhs);
+
+ /*! \brief Write the output to a file.
+ *
+ * \param output: `const string &` Path to the output to be written.
+ * \param format: `const string &` Output format (one of LEGACY or HDF5).
+ * \return result: `int` Exit code (0 if successful).
+ */
+ int write(const std::string &output, const std::string &format);
+
+#ifdef MPI_VERSION
+ /*! \brief Receive output data from worker processes.
+ *
+ * This function is invoked by the MPI rank-0 process to fetch the
+ * output data produced by higher rank processes. When calling this
+ * function, process 0 halts until a valid data chunk is transmitted
+ * by the queried process.
+ *
+ * \param mpidata: `const mixMPI*` Pointer to a `mixMPI` instance.
+ * \param pid: `int` Rank of the process that is transmitting data.
+ * \return result: `int` An exit code (0 for success).
+ */
+ int mpireceive(const mixMPI* mpidata, int pid);
+
+ /*! \brief Send output data to process 0.
+ *
+ * This function is invoked by non-zero ranked MPI processes when
+ * they are ready to send back the output data. When a process meets
+ * this function call, it halts until MPI process 0 asks for the
+ * data transmission.
+ *
+ * \param mpidata: `const mixMPI*` Pointer to a `mixMPI` instance.
+ * \param pid: `int` Rank of the process that is transmitting data.
+ * \return result: `int` An exit code (0 for success).
+ */
+ int mpisend(const mixMPI *mpidata);
+#endif // MPI_VERSION
+};
+// >>> END OF OUTPUT FOR INCLUSION <<<
+
#endif // INCLUDE_OUTPUTS_H_
diff --git a/src/inclusion/inclusion.cpp b/src/inclusion/inclusion.cpp
index 2d052e38bf87d1148446ee3673967bf90935b92b..25b6880f7e69e83fb1101863d91bd95f85df8c58 100644
--- a/src/inclusion/inclusion.cpp
+++ b/src/inclusion/inclusion.cpp
@@ -91,6 +91,10 @@
#include "../include/utils.h"
#endif
+#ifndef INCLUDE_OUTPUTS_H_
+#include "../include/outputs.h"
+#endif
+
using namespace std;
// >>> InclusionIterationData header <<< //
@@ -816,11 +820,11 @@ InclusionIterationData::~InclusionIterationData() {
* \param gconf: `GeometryConfiguration *` Pointer to a `GeometryConfiguration` object.
* \param sa: `ScatteringAngles *` Pointer to a `ScatteringAngles` object.
* \param cid: `InclusionIterationData *` Pointer to an `InclusionIterationData` object.
- * \param output: `VirtualAsciiFile *` Pointer to a `VirtualAsciiFile` object.
+ * \param oi: `InclusionOutputInfo *` Pointer to an `InclusionOutputInfo` object.
* \param output_path: `const string &` Path to the output directory.
* \param vtppoanp: `VirtualBinaryFile *` Pointer to a `VirtualBinaryFile` object.
*/
-int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, VirtualAsciiFile *output, const string& output_path, VirtualBinaryFile *vtppoanp);
+int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, InclusionOutputInfo *output, const string& output_path, VirtualBinaryFile *vtppoanp);
/*! \brief C++ implementation of INCLU
*
@@ -929,81 +933,17 @@ void inclusion(const string& config_file, const string& data_file, const string&
ScatteringAngles *p_scattering_angles = new ScatteringAngles(gconf);
double wp = sconf->wp;
// Open empty virtual ascii file for output
- VirtualAsciiFile *p_output = new VirtualAsciiFile();
- char virtual_line[256];
+ InclusionOutputInfo *p_output = new InclusionOutputInfo(sconf, gconf, mpidata);
InclusionIterationData *cid = new InclusionIterationData(gconf, sconf, mpidata, device_count);
const np_int ndi = cid->c1->nsph * cid->c1->nlim;
const np_int ndit = 2 * ndi;
logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)cid->c1->ndm) + " x " + to_string((int64_t)cid->c1->ndm) +".\n");
time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)cid->c1->ndm) + " x " + to_string((int64_t)cid->c1->ndm) +".\n");
- //==========================
- // Write a block of info to the ascii output file
- //==========================
- sprintf(virtual_line, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
- p_output->append_line(virtual_line);
-#ifdef USE_ILP64
- sprintf(virtual_line, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
- nsph, cid->c1->li, cid->c1->le, gconf->mxndm, gconf->in_pol, gconf->npnt,
- gconf->npntts, gconf->iavm, gconf->iavm
- );
-#else
- sprintf(virtual_line, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
- nsph, cid->c1->li, cid->c1->le, gconf->mxndm, gconf->in_pol, gconf->npnt,
- gconf->npntts, gconf->iavm, gconf->iavm
- );
-#endif // USE_ILP64
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
- p_output->append_line(virtual_line);
- for (int ri = 0; ri < nsph; ri++) {
- sprintf(virtual_line, "%17.8lE%17.8lE%17.8lE\n",
- gconf->get_sph_x(ri), gconf->get_sph_y(ri), gconf->get_sph_z(ri)
- );
- p_output->append_line(virtual_line);
- }
- sprintf(virtual_line, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
- p_output->append_line(virtual_line);
- sprintf(
- virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
- p_scattering_angles->th, p_scattering_angles->thstp,
- p_scattering_angles->thlst, p_scattering_angles->ths,
- p_scattering_angles->thsstp, p_scattering_angles->thslst
- );
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
- p_output->append_line(virtual_line);
- sprintf(
- virtual_line, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
- p_scattering_angles->ph, p_scattering_angles->phstp,
- p_scattering_angles->phlst, p_scattering_angles->phs,
- p_scattering_angles->phsstp, p_scattering_angles->phslst
- );
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " READ(IR,*)JWTM\n");
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " %5d\n", gconf->jwtm);
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " READ(ITIN)NSPHT\n");
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " READ(ITIN)(IOG(I),I=1,NSPH)\n");
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " READ(ITIN)(XIV(I),I=1,NXI)\n");
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " READ(ITIN)NSHL(I),ROS(I)\n");
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " \n");
- p_output->append_line(virtual_line);
instr(sconf, cid->c1);
thdps(cid->c1->lm, cid->zpv);
double exdc = sconf->exdc;
double exri = sqrt(exdc);
- sprintf(virtual_line, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
- p_output->append_line(virtual_line);
// Create an empty bynary file
VirtualBinaryFile *vtppoanp = new VirtualBinaryFile();
@@ -1037,10 +977,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
vtppoanp->append_line(VirtualBinaryLine(nphs));
if (sconf->idfc < 0) {
cid->vk = cid->xip * cid->wn;
- sprintf(virtual_line, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", cid->vk);
- p_output->append_line(virtual_line);
- sprintf(virtual_line, " \n");
- p_output->append_line(virtual_line);
+ p_output->vec_vk[0] = cid->vk;
}
// do the first iteration on jxi488 separately, since it seems to be different from the others
@@ -1092,8 +1029,6 @@ void inclusion(const string& config_file, const string& data_file, const string&
//==================================================
// do the first outputs here, so that I open here the new files, afterwards I only append
//==================================================
- p_output->write_to_disk(output_path + "/c_OINCLU");
- delete p_output;
vtppoanp->write_to_disk(output_path + "/c_TPPOAN");
delete vtppoanp;
@@ -1113,7 +1048,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
int myMPIstride = ompnumthreads;
int myMPIblock = ompnumthreads;
// Define here shared arrays of virtual ascii and binary files, so that thread 0 will be able to access them all later
- VirtualAsciiFile **p_outarray = NULL;
+ InclusionOutputInfo **p_outarray = NULL;
VirtualBinaryFile **vtppoanarray = NULL;
#ifdef USE_NVTX
@@ -1136,7 +1071,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
if (myompthread == 0) {
// Initialise some shared variables only on thread 0
- p_outarray = new VirtualAsciiFile*[ompnumthreads];
+ p_outarray = new InclusionOutputInfo*[ompnumthreads];
vtppoanarray = new VirtualBinaryFile*[ompnumthreads];
myMPIblock = ompnumthreads;
myMPIstride = myMPIblock;
@@ -1165,11 +1100,14 @@ void inclusion(const string& config_file, const string& data_file, const string&
// To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
InclusionIterationData *cid_2 = NULL;
- VirtualAsciiFile *p_output_2 = NULL;
+ InclusionOutputInfo *p_output_2 = NULL;
VirtualBinaryFile *vtppoanp_2 = NULL;
// for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
if (myompthread == 0) {
cid_2 = cid;
+ // OMP thread 0 of MPI process 0 holds the pointer to the full output structure
+ p_output_2 = p_output;
+ p_outarray[0] = p_output_2;
} else {
// this is not thread 0, so do create fresh copies of all local variables
cid_2 = new InclusionIterationData(*cid);
@@ -1185,16 +1123,24 @@ void inclusion(const string& config_file, const string& data_file, const string&
#pragma omp barrier
int myjxi488 = ixi488+myompthread;
// each thread opens new virtual files and stores their pointers in the shared array
- p_output_2 = new VirtualAsciiFile();
vtppoanp_2 = new VirtualBinaryFile();
// each thread puts a copy of the pointers to its virtual files in the shared arrays
- p_outarray[myompthread] = p_output_2;
vtppoanarray[myompthread] = vtppoanp_2;
#pragma omp barrier
// each MPI process handles a number of contiguous scales corresponding to its number of OMP threads at this omp level of parallelism
if (myjxi488 <= cid_2->number_of_scales) {
+ if (myompthread > 0) {
+ // UPDATE: non-0 threads need to allocate memory for one scale at a time.
+ p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, 1);
+ p_outarray[myompthread] = p_output_2;
+ }
int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
+ } else {
+ if (myompthread > 0) {
+ // If there is no input for this thread, set output pointer to NULL.
+ p_outarray[myompthread] = NULL;
+ }
}
#pragma omp barrier
@@ -1205,8 +1151,11 @@ void inclusion(const string& config_file, const string& data_file, const string&
// threads different from 0 append their virtual files to the one of thread 0, and delete them
if (myompthread == 0) {
for (int ti=1; ti<ompnumthreads; ti++) {
- p_outarray[0]->append(*(p_outarray[ti]));
- delete p_outarray[ti];
+ if (p_outarray[ti] != NULL) {
+ p_outarray[0]->insert(*(p_outarray[ti]));
+ delete p_outarray[ti];
+ p_outarray[ti] = NULL;
+ }
vtppoanarray[0]->append(*(vtppoanarray[ti]));
delete vtppoanarray[ti];
}
@@ -1217,8 +1166,8 @@ void inclusion(const string& config_file, const string& data_file, const string&
//==============================================
if (myompthread == 0) {
// thread 0 writes its virtual files, now including contributions from all threads, to disk, and deletes them
- p_outarray[0]->append_to_disk(output_path + "/c_OINCLU");
- delete p_outarray[0];
+ // p_outarray[0]->append_to_disk(output_path + "/c_OINCLU");
+ // delete p_outarray[0];
vtppoanarray[0]->append_to_disk(output_path + "/c_TPPOAN");
delete vtppoanarray[0];
@@ -1226,11 +1175,14 @@ void inclusion(const string& config_file, const string& data_file, const string&
if (mpidata->mpirunning) {
// only go through this if MPI has been actually used
for (int rr=1; rr<mpidata->nprocs; rr++) {
+ // get the data from process rr by receiving it in total memory structure
+ p_outarray[0]->mpireceive(mpidata, rr);
// get the data from process rr, creating a new virtual ascii file
- VirtualAsciiFile *p_output = new VirtualAsciiFile(mpidata, rr);
+ // VirtualAsciiFile *p_output = new VirtualAsciiFile(mpidata, rr);
// append to disk and delete virtual ascii file
- p_output->append_to_disk(output_path + "/c_OINCLU");
- delete p_output;
+ // p_output->append_to_disk(output_path + "/c_OINCLU");
+ // delete p_output;
+
// get the data from process rr, creating a new virtual binary file
VirtualBinaryFile *vtppoanp = new VirtualBinaryFile(mpidata, rr);
// append to disk and delete virtual binary file
@@ -1264,12 +1216,15 @@ void inclusion(const string& config_file, const string& data_file, const string&
delete cid_2;
}
delete p_scattering_angles;
- }
+ p_output->write(output_path + "/c_OINCLU.hd5", "HDF5");
+ p_output->write(output_path + "/c_OINCLU", "LEGACY");
+ delete p_output;
+ } // closes s_nsph == nsph check
else { // Sphere number inconsistency error.
throw UnrecognizedConfigurationException(
- "Inconsistent geometry and scatterer configurations."
- );
+ "Inconsistent geometry and scatterer configurations."
+ );
}
delete sconf;
@@ -1302,7 +1257,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
// Create this variable and initialise it with a default here, so that it is defined anyway, with or without OpenMP support enabled
int ompnumthreads = 1;
- VirtualAsciiFile **p_outarray = NULL;
+ InclusionOutputInfo **p_outarray = NULL;
VirtualBinaryFile **vtppoanarray = NULL;
int myjxi488startoffset;
int myMPIstride = ompnumthreads;
@@ -1325,13 +1280,13 @@ void inclusion(const string& config_file, const string& data_file, const string&
// receive myMPIstride sent by MPI process 0 to all processes
MPI_Bcast(&myMPIstride, 1, MPI_INT, 0, MPI_COMM_WORLD);
// allocate virtual files for each thread
- p_outarray = new VirtualAsciiFile*[ompnumthreads];
+ p_outarray = new InclusionOutputInfo*[ompnumthreads];
vtppoanarray = new VirtualBinaryFile*[ompnumthreads];
}
#pragma omp barrier
// To test parallelism, I will now start feeding this function with "clean" copies of the parameters, so that they will not be changed by previous iterations, and each one will behave as the first one. Define all (empty) variables here, so they have the correct scope, then they get different definitions depending on thread number
InclusionIterationData *cid_2 = NULL;
- VirtualAsciiFile *p_output_2 = NULL;
+ InclusionOutputInfo *p_output_2 = NULL;
VirtualBinaryFile *vtppoanp_2 = NULL;
// PLACEHOLDER
// for threads other than the 0, create distinct copies of all relevant data, while for thread 0 just define new references / pointers to the original ones
@@ -1349,25 +1304,41 @@ void inclusion(const string& config_file, const string& data_file, const string&
#pragma omp barrier
int myjxi488 = ixi488 + myjxi488startoffset + myompthread;
// each thread opens new virtual files and stores their pointers in the shared array
- p_output_2 = new VirtualAsciiFile();
vtppoanp_2 = new VirtualBinaryFile();
// each thread puts a copy of the pointers to its virtual files in the shared arrays
- p_outarray[myompthread] = p_output_2;
vtppoanarray[myompthread] = vtppoanp_2;
#pragma omp barrier
if (myompthread==0) logger->log("Syncing OpenMP threads and starting the loop on wavelengths\n");
// ok, now I can actually start the parallel calculations
// each MPI process handles a number of contiguous scales corresponding to its number of OMP threads at this omp level of parallelism
if (myjxi488 <= cid_2->number_of_scales) {
+ if (myompthread > 0) {
+ // UPDATE: non-0 threads need to allocate memory for one scale at a time.
+ p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, 1);
+ p_outarray[myompthread] = p_output_2;
+ } else {
+ // Thread 0 of non-zero MPI processes needs to allocate memory for the
+ // output of all threads.
+ p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, ompnumthreads);
+ p_outarray[0] = p_output_2;
+ }
int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
- } // close the OMP parallel for loop
+ } else {
+ if (myompthread > 0) {
+ // If there is no input for this thread, set the output pointer to NULL.
+ p_outarray[myompthread] = NULL;
+ }
+ }
#pragma omp barrier
// threads different from 0 append their virtual files to the one of thread 0, and delete them
if (myompthread == 0) {
for (int ti=1; ti<ompnumthreads; ti++) {
- p_outarray[0]->append(*(p_outarray[ti]));
- delete p_outarray[ti];
+ if (p_outarray[ti] != NULL) {
+ p_outarray[0]->insert(*(p_outarray[ti]));
+ delete p_outarray[ti];
+ p_outarray[ti] = NULL;
+ }
vtppoanarray[0]->append(*(vtppoanarray[ti]));
delete vtppoanarray[ti];
}
@@ -1407,7 +1378,8 @@ void inclusion(const string& config_file, const string& data_file, const string&
#endif
}
-int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, VirtualAsciiFile *output, const string& output_path, VirtualBinaryFile *vtppoanp) {
+int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, InclusionIterationData *cid, InclusionOutputInfo *output, const string& output_path, VirtualBinaryFile *vtppoanp) {
+ const dcomplex cc0 = 0.0 + I * 0.0;
int nxi = sconf->number_of_scales;
char virtual_line[256];
string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n";
@@ -1434,12 +1406,14 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
int last_configuration;
dcomplex ent, entn;
double enti;
+ int num_configs = sconf->configurations;
+ int ndirs = sa->nkks;
+ int oindex = -1;
+ int jindex = jxi488 - output->first_xi + 1;
#ifdef USE_NVTX
nvtxRangePush("Prepare matrix calculation");
#endif
- sprintf(virtual_line, "========== JXI =%3d ====================\n", jxi488);
- output->append_line(virtual_line);
double xi = sconf->get_scale(jxi488 - 1);
double exdc = sconf->exdc;
double exri = sqrt(exdc);
@@ -1448,14 +1422,14 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
if (idfc >= 0) {
cid->vk = xi * cid->wn;
vkarg = cid->vk;
- sprintf(virtual_line, " VK=%15.7lE, XI=%15.7lE\n", cid->vk, xi);
- output->append_line(virtual_line);
+ output->vec_vk[jindex - 1] = cid->vk;
+ output->vec_xi[jindex - 1] = xi;
// goes to 120
} else { // label 119
vkarg = xi * cid->vk;
cid->sqsfi = 1.0 / (xi * xi);
- sprintf(virtual_line, " XI=%15.7lE\n", xi);
- output->append_line(virtual_line);
+ output->vec_vk[jindex - 1] = cid->vk;
+ output->vec_xi[jindex - 1] = xi;
}
// label 120
double sze = vkarg * cid->extr;
@@ -1494,8 +1468,8 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
}
indme(i133, npnt, npntts, vkarg, ent, enti, entn, jer, lcalc, cid->arg, cid->c1);
if (jer != 0) {
- sprintf(virtual_line, " STOP IN INDME\n");
- output->append_line(virtual_line);
+ output->vec_ier[jindex - 1] = 1;
+ output->vec_jxi[jindex - 1] = -jxi488;
message = "ERROR: indme failed with error code " + to_string(jer) + ".\n";
logger->log(message, LOG_ERRO);
delete logger;
@@ -1506,8 +1480,8 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
} // i133 loop
ospv(cid->c1, vkarg, sze, exri, entn, enti, jer, lcalc, cid->arg);
if (jer != 0) {
- sprintf(virtual_line, " STOP IN OSPV\n");
- output->append_line(virtual_line);
+ output->vec_ier[jindex - 1] = 2;
+ output->vec_jxi[jindex - 1] = -jxi488;
message = "ERROR: ospv failed with error code " + to_string(jer) + ".\n";
logger->log(message, LOG_ERRO);
delete logger;
@@ -1610,19 +1584,24 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
}
}
// label 156: continue from here
+ last_configuration = 0;
for (int i168 = 1; i168 <= nsph; i168++) {
if (cid->c1->iog[i168 - 1] >= i168) {
+ int i = i168 - 1;
+ last_configuration++;
+ oindex = (jindex - 1) * (num_configs + 1) + last_configuration - 1;
if (cid->c1->nshl[i168 - 1] != 1) {
- sprintf(virtual_line, " SPHERE N.%2d: SIZE=%15.7lE\n", i168, cid->c1->vsz[i168 - 1]);
- output->append_line(virtual_line);
+ output->vec_sphere_ref_indices[oindex] = cc0;
+ output->vec_sphere_sizes[oindex] = cid->c1->vsz[i];
} else {
- sprintf(virtual_line, " SPHERE N.%2d: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", i168, cid->c1->vsz[i168 - 1], real(cid->c1->vkt[i168 - 1]), imag(cid->c1->vkt[i168 - 1]));
- output->append_line(virtual_line);
+ output->vec_sphere_ref_indices[oindex] = cid->c1->vkt[i];
+ output->vec_sphere_sizes[oindex] = cid->c1->vsz[i];
}
}
} // i168 loop
- sprintf(virtual_line, " EXT. SPHERE: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", sze, real(entn), imag(entn));
- output->append_line(virtual_line);
+ oindex = (jindex - 1) * (num_configs + 1) + num_configs;
+ output->vec_sphere_sizes[oindex] = sze;
+ output->vec_sphere_ref_indices[oindex] = entn;
// label 160
double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0);
double csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcs;
@@ -1642,6 +1621,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
nvtxRangePush("Angle loop");
#endif
double th = sa->th;
+ int done_dirs = 0;
for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable?
double ph = sa->ph;
double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
@@ -1767,8 +1747,6 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
vtppoanp->append_line(VirtualBinaryLine(value));
}
}
- sprintf(virtual_line, " ENSEMBLE AVERAGE, MODE%2d\n", iavm);
- output->append_line(virtual_line);
int jlr = 2;
for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
int ipol = (ilr210 % 2 == 0) ? 1 : -1;
@@ -1784,68 +1762,72 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
double qschum = imag(s0m) * csch;
double pschum = real(s0m) * csch;
double s0magm = cabs(s0m) * cs0;
- if (inpol == 0) {
- sprintf(virtual_line, " LIN %2d\n", ipol);
- output->append_line(virtual_line);
- } else { // label 206
- sprintf(virtual_line, " CIRC %2d\n", ipol);
- output->append_line(virtual_line);
- }
+ // if (inpol == 0) {
+ // sprintf(virtual_line, " LIN %2d\n", ipol);
+ // output->append_line(virtual_line);
+ // } else { // label 206
+ // sprintf(virtual_line, " CIRC %2d\n", ipol);
+ // output->append_line(virtual_line);
+ // }
// label 208
- sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
- output->append_line(virtual_line);
- sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n", scasm, abssm, extsm, albdm);
- output->append_line(virtual_line);
- double alamb = 2.0 * 3.141592653589793238 / cid->vk;
- sprintf(virtual_line, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n", ipol, alamb, scasm, abssm, extsm);
- output->append_line(virtual_line);
- sprintf(virtual_line, " ---- SCS/GS -- ABC/GS -- EXS/GS ---\n");
- output->append_line(virtual_line);
- sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n", qscam, qabsm, qextm);
- output->append_line(virtual_line);
- sprintf(
- virtual_line, " FSAS(%1d,%1d)=%15.7lE%15.7lE FSAS(%1d,%1d)=%15.7lE%15.7lE\n",
- ilr210, ilr210, real(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]),
- imag(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,
- real(cid->c1->fsacm[jlr - 1][ilr210 - 1]), imag(cid->c1->fsacm[jlr - 1][ilr210 - 1])
- );
- output->append_line(virtual_line);
- sprintf(virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschum, pschum, s0magm);
- output->append_line(virtual_line);
double rapr = cid->c1->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1];
double cosav = cid->gapm[2][ilr210 - 1] / cid->c1->scscm[ilr210 - 1];
double fz = rapr;
- sprintf(virtual_line, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
- output->append_line(virtual_line);
- sprintf(virtual_line, " Fk=%15.7lE\n", fz);
- output->append_line(virtual_line);
+ if (ipol == -1) {
+ output->vec_scs1[jindex - 1] = scasm;
+ output->vec_abs1[jindex - 1] = abssm;
+ output->vec_exs1[jindex - 1] = extsm;
+ output->vec_albeds1[jindex - 1] = albdm;
+ output->vec_scsrt1[jindex - 1] = qscam;
+ output->vec_absrt1[jindex - 1] = qabsm;
+ output->vec_exsrt1[jindex - 1] = qextm;
+ output->vec_fsas11[jindex - 1] = cid->c1->fsacm[0][0];
+ output->vec_fsas21[jindex - 1] = cid->c1->fsacm[1][0];
+ output->vec_qschu1[jindex -1] = qschum;
+ output->vec_pschu1[jindex -1] = pschum;
+ output->vec_s0mag1[jindex -1] = s0magm;
+ output->vec_cosav1[jindex -1] = cosav;
+ output->vec_raprs1[jindex -1] = rapr;
+ output->vec_fk1[jindex -1] = fz;
+ } else if (ipol == 1) {
+ output->vec_scs2[jindex - 1] = scasm;
+ output->vec_abs2[jindex - 1] = abssm;
+ output->vec_exs2[jindex - 1] = extsm;
+ output->vec_albeds2[jindex - 1] = albdm;
+ output->vec_scsrt2[jindex - 1] = qscam;
+ output->vec_absrt2[jindex - 1] = qabsm;
+ output->vec_exsrt2[jindex - 1] = qextm;
+ output->vec_fsas22[jindex - 1] = cid->c1->fsacm[1][1];
+ output->vec_fsas12[jindex - 1] = cid->c1->fsacm[0][1];
+ output->vec_qschu2[jindex -1] = qschum;
+ output->vec_pschu2[jindex -1] = pschum;
+ output->vec_s0mag2[jindex -1] = s0magm;
+ output->vec_cosav2[jindex -1] = cosav;
+ output->vec_raprs2[jindex -1] = rapr;
+ output->vec_fk2[jindex -1] = fz;
+ }
} // ilr210 loop
- double rmbrif = (real(cid->c1->fsacm[0][0]) - real(cid->c1->fsacm[1][1])) / real(cid->c1->fsacm[0][0]);
- double rmdchr = (imag(cid->c1->fsacm[0][0]) - imag(cid->c1->fsacm[1][1])) / imag(cid->c1->fsacm[0][0]);
- sprintf(virtual_line, " (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n", rmbrif);
- output->append_line(virtual_line);
- sprintf(virtual_line, " (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n", rmdchr);
- output->append_line(virtual_line);
}
// label 212
- sprintf(virtual_line, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n", jth486, jph484, jths, jphs);
- output->append_line(virtual_line);
- sprintf(virtual_line, " TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n", th, ph, ths, phs);
- output->append_line(virtual_line);
- sprintf(virtual_line, " SCAND=%10.3lE\n", cid->scan);
- output->append_line(virtual_line);
- sprintf(virtual_line, " CFMP=%15.7lE, SFMP=%15.7lE\n", cid->cfmp, cid->sfmp);
- output->append_line(virtual_line);
- sprintf(virtual_line, " CFSP=%15.7lE, SFSP=%15.7lE\n", cid->cfsp, cid->sfsp);
- output->append_line(virtual_line);
+ output->vec_dir_scand[done_dirs] = cid->scan;
+ output->vec_dir_cfmp[done_dirs] = cid->cfmp;
+ output->vec_dir_sfmp[done_dirs] = cid->sfmp;
+ output->vec_dir_cfsp[done_dirs] = cid->cfsp;
+ output->vec_dir_sfsp[done_dirs] = cid->sfsp;
if (isam >= 0) {
- sprintf(virtual_line, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n", cid->un[0], cid->un[1], cid->un[2]);
- output->append_line(virtual_line);
- sprintf(virtual_line, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n", cid->uns[0], cid->uns[1], cid->uns[2]);
- output->append_line(virtual_line);
+ output->vec_dir_un[3 * done_dirs] = cid->un[0];
+ output->vec_dir_un[3 * done_dirs + 1] = cid->un[1];
+ output->vec_dir_un[3 * done_dirs + 2] = cid->un[2];
+ output->vec_dir_uns[3 * done_dirs] = cid->uns[0];
+ output->vec_dir_uns[3 * done_dirs + 1] = cid->uns[1];
+ output->vec_dir_uns[3 * done_dirs + 2] = cid->uns[2];
} else { // label 214
- sprintf(virtual_line, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", cid->un[0], cid->un[1], cid->un[2]);
- output->append_line(virtual_line);
+ output->vec_dir_un[3 * done_dirs] = cid->un[0];
+ output->vec_dir_un[3 * done_dirs + 1] = cid->un[1];
+ output->vec_dir_un[3 * done_dirs + 2] = cid->un[2];
+ output->vec_dir_uns[3 * done_dirs] = 0.0;
+ output->vec_dir_uns[3 * done_dirs + 1] = 0.0;
+ output->vec_dir_uns[3 * done_dirs + 2] = 0.0;
}
// label 220
pcros(cid->vk, exri, cid->c1);
@@ -1926,8 +1908,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
vtppoanp->append_line(VirtualBinaryLine(value));
}
}
- sprintf(virtual_line, " SINGLE SCATTERER\n");
- output->append_line(virtual_line);
+ oindex = (jindex - 1) * ndirs + done_dirs;
int jlr = 2;
for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
int ipol = (ilr290 % 2 == 0) ? 1 : -1;
@@ -1943,53 +1924,38 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
double qschu = imag(s0) * csch;
double pschu = real(s0) * csch;
double s0mag = cabs(s0) * cs0;
- if (inpol == 0) {
- sprintf(virtual_line, " LIN %2d\n", ipol);
- output->append_line(virtual_line);
- } else { // label 273
- sprintf(virtual_line, " CIRC %2d\n", ipol);
- output->append_line(virtual_line);
- }
// label 275
- sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
- output->append_line(virtual_line);
- sprintf(
- virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
- scasec, abssec, extsec, albedc
- );
- output->append_line(virtual_line);
- double alam = 2.0 * 3.141592653589793238 / cid->vk;
- sprintf(virtual_line, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n", ipol, alam, scasec, abssec, extsec);
- sprintf(virtual_line, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
- output->append_line(virtual_line);
- sprintf(
- virtual_line, " %14.7lE%15.7lE%15.7lE\n",
- qsca, qabs, qext
- );
- output->append_line(virtual_line);
- sprintf(
- virtual_line,
- " FSAS(%1d,%1d)=%15.7lE%15.7lE FSAS(%1d,%1d)=%15.7lE%15.7lE\n",
- ilr290, ilr290, real(cid->c1->fsac[ilr290 - 1][ilr290 - 1]),
- imag(cid->c1->fsac[ilr290 - 1][ilr290 - 1]), jlr, ilr290,
- real(cid->c1->fsac[jlr - 1][ilr290 - 1]),
- imag(cid->c1->fsac[jlr - 1][ilr290 - 1])
- );
- output->append_line(virtual_line);
- sprintf(
- virtual_line,
- " SAS(%1d,%1d)=%15.7lE%15.7lE SAS(%1d,%1d)=%15.7lE%15.7lE\n",
- ilr290, ilr290, real(cid->c1->sac[ilr290 - 1][ilr290 - 1]),
- imag(cid->c1->sac[ilr290 - 1][ilr290 - 1]), jlr, ilr290,
- real(cid->c1->sac[jlr - 1][ilr290 - 1]),
- imag(cid->c1->sac[jlr - 1][ilr290 - 1])
- );
- output->append_line(virtual_line);
- sprintf(
- virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
- qschu, pschu, s0mag
- );
- output->append_line(virtual_line);
+ if (ipol == -1) {
+ output->vec_dir_scs1[oindex] = scasec;
+ output->vec_dir_abs1[oindex] = abssec;
+ output->vec_dir_exs1[oindex] = extsec;
+ output->vec_dir_albeds1[oindex] = albedc;
+ output->vec_dir_scsrt1[oindex] = qsca;
+ output->vec_dir_absrt1[oindex] = qabs;
+ output->vec_dir_exsrt1[oindex] = qext;
+ output->vec_dir_fsas11[oindex] = cid->c1->fsac[0][0];
+ output->vec_dir_fsas21[oindex] = cid->c1->fsac[1][0];
+ output->vec_dir_sas11[oindex] = cid->c1->sac[0][0];
+ output->vec_dir_sas21[oindex] = cid->c1->sac[1][0];
+ output->vec_dir_qschu1[oindex] = qschu;
+ output->vec_dir_pschu1[oindex] = pschu;
+ output->vec_dir_s0mag1[oindex] = s0mag;
+ } else if (ipol == 1) {
+ output->vec_dir_scs2[oindex] = scasec;
+ output->vec_dir_abs2[oindex] = abssec;
+ output->vec_dir_exs2[oindex] = extsec;
+ output->vec_dir_albeds2[oindex] = albedc;
+ output->vec_dir_scsrt2[oindex] = qsca;
+ output->vec_dir_absrt2[oindex] = qabs;
+ output->vec_dir_exsrt2[oindex] = qext;
+ output->vec_dir_fsas22[oindex] = cid->c1->fsac[1][1];
+ output->vec_dir_fsas12[oindex] = cid->c1->fsac[0][1];
+ output->vec_dir_sas22[oindex] = cid->c1->sac[1][1];
+ output->vec_dir_sas12[oindex] = cid->c1->sac[0][1];
+ output->vec_dir_qschu2[oindex] = qschu;
+ output->vec_dir_pschu2[oindex] = pschu;
+ output->vec_dir_s0mag2[oindex] = s0mag;
+ }
bool goto290 = isam >= 0 && (jths > 1 || jphs > 1);
if (!goto290) {
cid->gapv[0] = cid->gap[0][ilr290 - 1];
@@ -1999,12 +1965,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
double scatts = cid->c1->scsc[ilr290 - 1];
double rapr, cosav, fp, fn, fk, fx, fy, fz;
rftr(cid->u, cid->up, cid->un, cid->gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz);
- sprintf(virtual_line, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
- output->append_line(virtual_line);
- sprintf(virtual_line, " Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n", fp, fn, fk);
- output->append_line(virtual_line);
- sprintf(virtual_line, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n", fx, fy, fz);
- output->append_line(virtual_line);
+ if (ipol == -1) {
+ output->vec_dir_cosav1[oindex] = cosav;
+ output->vec_dir_rapr1[oindex] = rapr;
+ output->vec_dir_fl1[oindex] = fp;
+ output->vec_dir_fr1[oindex] = fn;
+ output->vec_dir_fk1[oindex] = fk;
+ output->vec_dir_fx1[oindex] = fx;
+ output->vec_dir_fy1[oindex] = fy;
+ output->vec_dir_fz1[oindex] = fz;
+ } else if (ipol == 1) {
+ output->vec_dir_cosav2[oindex] = cosav;
+ output->vec_dir_rapr2[oindex] = rapr;
+ output->vec_dir_fl2[oindex] = fp;
+ output->vec_dir_fr2[oindex] = fn;
+ output->vec_dir_fk2[oindex] = fk;
+ output->vec_dir_fx2[oindex] = fx;
+ output->vec_dir_fy2[oindex] = fy;
+ output->vec_dir_fz2[oindex] = fz;
+ }
cid->tqev[0] = cid->tqce[ilr290 - 1][0];
cid->tqev[1] = cid->tqce[ilr290 - 1][1];
cid->tqev[2] = cid->tqce[ilr290 - 1][2];
@@ -2013,45 +1992,48 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
cid->tqsv[2] = cid->tqcs[ilr290 - 1][2];
double tep, ten, tek, tsp, tsn, tsk;
tqr(cid->u, cid->up, cid->un, cid->tqev, cid->tqsv, tep, ten, tek, tsp, tsn, tsk);
- sprintf(virtual_line, " TQEl=%15.7lE, TQEr=%15.7lE, TQEk=%15.7lE\n", tep, ten, tek);
- output->append_line(virtual_line);
- sprintf(virtual_line, " TQSl=%15.7lE, TQSr=%15.7lE, TQSk=%15.7lE\n", tsp, tsn, tsk);
- output->append_line(virtual_line);
- sprintf(
- virtual_line, " TQEx=%15.7lE, TQEy=%15.7lE, TQEz=%15.7lE\n",
- cid->tqce[ilr290 - 1][0], cid->tqce[ilr290 - 1][1], cid->tqce[ilr290 - 1][2]
- );
- output->append_line(virtual_line);
- sprintf(
- virtual_line, " TQSx=%15.7lE, TQSy=%15.7lE, TQSz=%15.7lE\n",
- cid->tqcs[ilr290 - 1][0], cid->tqcs[ilr290 - 1][1], cid->tqcs[ilr290 - 1][2]
- );
- output->append_line(virtual_line);
+ if (ipol == -1) {
+ output->vec_dir_tqel1[oindex] = tep;
+ output->vec_dir_tqer1[oindex] = ten;
+ output->vec_dir_tqek1[oindex] = tek;
+ output->vec_dir_tqsl1[oindex] = tsp;
+ output->vec_dir_tqsr1[oindex] = tsn;
+ output->vec_dir_tqsk1[oindex] = tsk;
+ output->vec_dir_tqex1[oindex] = cid->tqce[0][0];
+ output->vec_dir_tqey1[oindex] = cid->tqce[0][1];
+ output->vec_dir_tqez1[oindex] = cid->tqce[0][2];
+ output->vec_dir_tqsx1[oindex] = cid->tqcs[0][0];
+ output->vec_dir_tqsy1[oindex] = cid->tqcs[0][1];
+ output->vec_dir_tqsz1[oindex] = cid->tqcs[0][2];
+ } else if (ipol == 1) {
+ output->vec_dir_tqel2[oindex] = tep;
+ output->vec_dir_tqer2[oindex] = ten;
+ output->vec_dir_tqek2[oindex] = tek;
+ output->vec_dir_tqsl2[oindex] = tsp;
+ output->vec_dir_tqsr2[oindex] = tsn;
+ output->vec_dir_tqsk2[oindex] = tsk;
+ output->vec_dir_tqex2[oindex] = cid->tqce[1][0];
+ output->vec_dir_tqey2[oindex] = cid->tqce[1][1];
+ output->vec_dir_tqez2[oindex] = cid->tqce[1][2];
+ output->vec_dir_tqsx2[oindex] = cid->tqcs[1][0];
+ output->vec_dir_tqsy2[oindex] = cid->tqcs[1][1];
+ output->vec_dir_tqsz2[oindex] = cid->tqcs[1][2];
+ }
}
} //ilr290 loop
- double rbirif = (real(cid->c1->fsac[0][0]) - real(cid->c1->fsac[1][1])) / real(cid->c1->fsac[0][0]);
- double rdichr = (imag(cid->c1->fsac[0][0]) - imag(cid->c1->fsac[1][1])) / imag(cid->c1->fsac[0][0]);
- sprintf(virtual_line, " (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n", rbirif);
- output->append_line(virtual_line);
- sprintf(virtual_line, " (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n", rdichr);
- output->append_line(virtual_line);
- sprintf(virtual_line, " MULL\n");
- output->append_line(virtual_line);
for (int i = 0; i < 4; i++) {
- sprintf(
- virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
- cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->cmul[i][3]
- );
- output->append_line(virtual_line);
+ oindex = 16 * (jindex - 1) * ndirs + 16 * done_dirs + 4 * i;
+ output->vec_dir_mull[oindex] = cid->cmul[i][0];
+ output->vec_dir_mull[oindex + 1] = cid->cmul[i][1];
+ output->vec_dir_mull[oindex + 2] = cid->cmul[i][2];
+ output->vec_dir_mull[oindex + 3] = cid->cmul[i][3];
}
- sprintf(virtual_line, " MULLLR\n");
- output->append_line(virtual_line);
for (int i = 0; i < 4; i++) {
- sprintf(
- virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
- cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3]
- );
- output->append_line(virtual_line);
+ oindex = 16 * (jindex - 1) * ndirs + 16 * done_dirs + 4 * i;
+ output->vec_dir_mulllr[oindex] = cid->cmullr[i][0];
+ output->vec_dir_mulllr[oindex + 1] = cid->cmullr[i][1];
+ output->vec_dir_mulllr[oindex + 2] = cid->cmullr[i][2];
+ output->vec_dir_mulllr[oindex + 3] = cid->cmullr[i][3];
}
if (iavm != 0) {
mmulc(cid->c1->vintm, cid->cmullr, cid->cmul);
@@ -2069,37 +2051,25 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
vtppoanp->append_line(VirtualBinaryLine(value));
}
}
- sprintf(virtual_line, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
- output->append_line(virtual_line);
- if (inpol == 0) {
- sprintf(virtual_line, " LIN\n");
- output->append_line(virtual_line);
- } else { // label 316
- sprintf(virtual_line, " CIRC\n");
- output->append_line(virtual_line);
- }
// label 318
- sprintf(virtual_line, " MULC\n");
- output->append_line(virtual_line);
for (int i = 0; i < 4; i++) {
- sprintf(
- virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
- cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->cmul[i][3]
- );
- output->append_line(virtual_line);
+ oindex = 16 * (jindex - 1) + 4 * i; // if IAVM fails, try adding directions
+ output->vec_dir_mull[oindex] = cid->cmul[i][0];
+ output->vec_dir_mull[oindex + 1] = cid->cmul[i][1];
+ output->vec_dir_mull[oindex + 2] = cid->cmul[i][2];
+ output->vec_dir_mull[oindex + 3] = cid->cmul[i][3];
}
- sprintf(virtual_line, " MULCLR\n");
- output->append_line(virtual_line);
for (int i = 0; i < 4; i++) {
- sprintf(
- virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
- cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3]
- );
- output->append_line(virtual_line);
+ oindex = 16 * (jindex - 1) + 4 * i; // if IAVM fails, try adding directions
+ output->vec_dir_mulllr[oindex] = cid->cmul[i][0];
+ output->vec_dir_mulllr[oindex + 1] = cid->cmullr[i][1];
+ output->vec_dir_mulllr[oindex + 2] = cid->cmullr[i][2];
+ output->vec_dir_mulllr[oindex + 3] = cid->cmullr[i][3];
}
}
// label 420, continues jphs loop
if (isam < 1) phs += sa->phsstp;
+ done_dirs++;
} // jphs loop, labeled 480
if (isam <= 1) thsl += sa->thsstp;
} // jths loop, labeled 482
@@ -2110,6 +2080,7 @@ int inclusion_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryCo
#ifdef USE_NVTX
nvtxRangePop();
#endif
+ if (jer == 0) output->vec_jxi[jindex - 1] = jxi488;
interval_end = chrono::high_resolution_clock::now();
elapsed = interval_end - interval_start;
message = "INFO: angle loop for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
diff --git a/src/libnptm/outputs.cpp b/src/libnptm/outputs.cpp
index 1f50f1a06954763593390f29cb65b96e34c64caf..4d6e61a2fb6273f9e665e038098d6f76ad0a868b 100644
--- a/src/libnptm/outputs.cpp
+++ b/src/libnptm/outputs.cpp
@@ -56,7 +56,7 @@
using namespace std;
-// >> ClusterOutputInfo CLASS IMPLEMENTATION <<
+// >>> ClusterOutputInfo CLASS IMPLEMENTATION <<<
ClusterOutputInfo::ClusterOutputInfo(
ScattererConfiguration *sc, GeometryConfiguration *gc,
const mixMPI *mpidata, int first_xi, int xi_length
@@ -873,7 +873,7 @@ long ClusterOutputInfo::compute_size(
) {
long result = sizeof(np_int);
result += 21 * sizeof(int);
- result += 56 * sizeof(double);
+ result += 14 * sizeof(double);
result += 121 * sizeof(long);
int _nsph = gc->number_of_spheres;
double _th = gc->in_theta_start;
@@ -928,7 +928,7 @@ long ClusterOutputInfo::compute_size(
long ClusterOutputInfo::compute_size() {
long result = sizeof(np_int);
result += 21 * sizeof(int);
- result += 56 * sizeof(double);
+ result += 14 * sizeof(double);
result += 121 * sizeof(long);
result += 3 * nsph * sizeof(double); // sphere coordinate vectors
result += xi_block_size * sizeof(int); // scale index vector
@@ -2739,8 +2739,2039 @@ int ClusterOutputInfo::mpisend(const mixMPI *mpidata) {
MPI_Send(vec_dir_mulc, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
MPI_Send(vec_dir_mulclr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
}
- result = 0;
return result;
}
#endif //MPI_VERSION
-// >> END OF ClusterOutputInfo CLASS IMPLEMENTATION <<
+// >>> END OF ClusterOutputInfo CLASS IMPLEMENTATION <<<
+
+// >>> InclusionOutputInfo CLASS IMPLEMENTATION <<<
+InclusionOutputInfo::InclusionOutputInfo(
+ ScattererConfiguration *sc, GeometryConfiguration *gc,
+ const mixMPI *mpidata, int first_xi, int xi_length
+) {
+ nsph = gc->number_of_spheres;
+ li = gc->li;
+ le = gc->le;
+ lm = gc->l_max;
+ mxndm = gc->mxndm;
+ inpol = gc->in_pol;
+ npnt = gc->npnt;
+ npntts = gc->npntts;
+ iavm = gc->iavm;
+ isam = gc->isam;
+ jwtm = gc->jwtm;
+ // Get spherical constituent coordinates
+ vec_x_coords = new double[nsph];
+ vec_y_coords = new double[nsph];
+ vec_z_coords = new double[nsph];
+ for (int nsi = 0; nsi < nsph; nsi++) {
+ vec_x_coords[nsi] = gc->get_sph_x(nsi);
+ vec_y_coords[nsi] = gc->get_sph_y(nsi);
+ vec_z_coords[nsi] = gc->get_sph_z(nsi);
+ }
+ // Get directional information
+ th = gc->in_theta_start;
+ thstp = gc->in_theta_step;
+ thlst = gc->in_theta_end;
+ ths = gc->sc_theta_start;
+ thsstp = gc->sc_theta_step;
+ thslst = gc->sc_theta_end;
+ _num_theta = (thstp == 0.0) ? 1 : 1 + int((thlst - th) / thstp);
+ _num_thetas = (thsstp == 0.0) ? 1 : 1 + int((thslst - ths) / thsstp);
+ ph = gc->in_phi_start;
+ phstp = gc->in_phi_step;
+ phlst = gc->in_phi_end;
+ phs = gc->sc_phi_start;
+ phsstp = gc->sc_phi_step;
+ phslst = gc->sc_phi_end;
+ _num_phi = (phstp == 0.0) ? 1 : 1 + int((phlst - ph) / phstp);
+ _num_phis = (phsstp == 0.0) ? 1 : 1 + int((phslst - phs) / phsstp);
+ ndirs = _num_theta * _num_thetas * _num_phi * _num_phis;
+ // Get scattering problem configuration
+ double exdc = sc->exdc;
+ exri = sqrt(exdc);
+ idfc = sc->idfc;
+ nxi = sc->number_of_scales;
+ _first_xi = first_xi;
+ xi_block_size = (xi_length == 0) ? nxi - (_first_xi - 1) : xi_length;
+ vec_jxi = new int[xi_block_size]();
+ vec_ier = new short[xi_block_size]();
+ vec_vk = new double[xi_block_size]();
+ vec_xi = new double[xi_block_size]();
+ configurations = sc->configurations;
+ vec_sphere_sizes = new double[xi_block_size * (configurations + 1)]();
+ vec_sphere_ref_indices = new dcomplex[xi_block_size * (configurations + 1)]();
+ vec_scs1 = new double[xi_block_size]();
+ vec_scs2 = new double[xi_block_size]();
+ vec_abs1 = new double[xi_block_size]();
+ vec_abs2 = new double[xi_block_size]();
+ vec_exs1 = new double[xi_block_size]();
+ vec_exs2 = new double[xi_block_size]();
+ vec_albeds1 = new double[xi_block_size]();
+ vec_albeds2 = new double[xi_block_size]();
+ vec_scsrt1 = new double[xi_block_size]();
+ vec_scsrt2 = new double[xi_block_size]();
+ vec_absrt1 = new double[xi_block_size]();
+ vec_absrt2 = new double[xi_block_size]();
+ vec_exsrt1 = new double[xi_block_size]();
+ vec_exsrt2 = new double[xi_block_size]();
+ vec_qschu1 = new double[xi_block_size]();
+ vec_qschu2 = new double[xi_block_size]();
+ vec_pschu1 = new double[xi_block_size]();
+ vec_pschu2 = new double[xi_block_size]();
+ vec_s0mag1 = new double[xi_block_size]();
+ vec_s0mag2 = new double[xi_block_size]();
+ vec_cosav1 = new double[xi_block_size]();
+ vec_cosav2 = new double[xi_block_size]();
+ vec_raprs1 = new double[xi_block_size]();
+ vec_raprs2 = new double[xi_block_size]();
+ vec_fk1 = new double[xi_block_size]();
+ vec_fk2 = new double[xi_block_size]();
+ vec_fsas11 = new dcomplex[xi_block_size]();
+ vec_fsas21 = new dcomplex[xi_block_size]();
+ vec_fsas22 = new dcomplex[xi_block_size]();
+ vec_fsas12 = new dcomplex[xi_block_size]();
+ vec_dir_tidg = new double[_num_theta];
+ vec_dir_pidg = new double[_num_phi];
+ vec_dir_tsdg = new double[_num_thetas];
+ vec_dir_psdg = new double[_num_phis];
+ // Initialize directions (they are scale-independent)
+ double cti = th, cpi = ph, cts = ths, cps = phs;
+ for (int di = 0; di < _num_theta; di++) {
+ vec_dir_tidg[di] = cti;
+ cti += thstp;
+ }
+ for (int di = 0; di < _num_thetas; di++) {
+ vec_dir_tsdg[di] = cts;
+ cts += thsstp;
+ }
+ for (int di = 0; di < _num_phi; di++) {
+ vec_dir_pidg[di] = cpi;
+ cpi += phstp;
+ }
+ for (int di = 0; di < _num_phis; di++) {
+ vec_dir_psdg[di] = cps;
+ cps += phsstp;
+ }
+ vec_dir_scand = new double[ndirs]();
+ vec_dir_cfmp = new double[ndirs]();
+ vec_dir_sfmp = new double[ndirs]();
+ vec_dir_cfsp = new double[ndirs]();
+ vec_dir_sfsp = new double[ndirs]();
+ vec_dir_un = new double[3 * ndirs]();
+ vec_dir_uns = new double[3 * ndirs]();
+ vec_dir_scs1 = new double[ndirs * xi_block_size]();
+ vec_dir_scs2 = new double[ndirs * xi_block_size]();
+ vec_dir_abs1 = new double[ndirs * xi_block_size]();
+ vec_dir_abs2 = new double[ndirs * xi_block_size]();
+ vec_dir_exs1 = new double[ndirs * xi_block_size]();
+ vec_dir_exs2 = new double[ndirs * xi_block_size]();
+ vec_dir_albeds1 = new double[ndirs * xi_block_size]();
+ vec_dir_albeds2 = new double[ndirs * xi_block_size]();
+ vec_dir_scsrt1 = new double[ndirs * xi_block_size]();
+ vec_dir_scsrt2 = new double[ndirs * xi_block_size]();
+ vec_dir_absrt1 = new double[ndirs * xi_block_size]();
+ vec_dir_absrt2 = new double[ndirs * xi_block_size]();
+ vec_dir_exsrt1 = new double[ndirs * xi_block_size]();
+ vec_dir_exsrt2 = new double[ndirs * xi_block_size]();
+ vec_dir_fsas11 = new dcomplex[ndirs * xi_block_size]();
+ vec_dir_fsas21 = new dcomplex[ndirs * xi_block_size]();
+ vec_dir_fsas12 = new dcomplex[ndirs * xi_block_size]();
+ vec_dir_fsas22 = new dcomplex[ndirs * xi_block_size]();
+ vec_dir_sas11 = new dcomplex[ndirs * xi_block_size]();
+ vec_dir_sas21 = new dcomplex[ndirs * xi_block_size]();
+ vec_dir_sas12 = new dcomplex[ndirs * xi_block_size]();
+ vec_dir_sas22 = new dcomplex[ndirs * xi_block_size]();
+ vec_dir_qschu1 = new double[ndirs * xi_block_size]();
+ vec_dir_qschu2 = new double[ndirs * xi_block_size]();
+ vec_dir_pschu1 = new double[ndirs * xi_block_size]();
+ vec_dir_pschu2 = new double[ndirs * xi_block_size]();
+ vec_dir_s0mag1 = new double[ndirs * xi_block_size]();
+ vec_dir_s0mag2 = new double[ndirs * xi_block_size]();
+ vec_dir_cosav1 = new double[ndirs * xi_block_size]();
+ vec_dir_cosav2 = new double[ndirs * xi_block_size]();
+ vec_dir_rapr1 = new double[ndirs * xi_block_size]();
+ vec_dir_rapr2 = new double[ndirs * xi_block_size]();
+ vec_dir_fl1 = new double[ndirs * xi_block_size]();
+ vec_dir_fl2 = new double[ndirs * xi_block_size]();
+ vec_dir_fr1 = new double[ndirs * xi_block_size]();
+ vec_dir_fr2 = new double[ndirs * xi_block_size]();
+ vec_dir_fk1 = new double[ndirs * xi_block_size]();
+ vec_dir_fk2 = new double[ndirs * xi_block_size]();
+ vec_dir_fx1 = new double[ndirs * xi_block_size]();
+ vec_dir_fx2 = new double[ndirs * xi_block_size]();
+ vec_dir_fy1 = new double[ndirs * xi_block_size]();
+ vec_dir_fy2 = new double[ndirs * xi_block_size]();
+ vec_dir_fz1 = new double[ndirs * xi_block_size]();
+ vec_dir_fz2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqel1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqel2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqer1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqer2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqek1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqek2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqex1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqex2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqey1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqey2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqez1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqez2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsl1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsl2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsr1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsr2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsk1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsk2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsx1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsx2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsy1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsy2 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsz1 = new double[ndirs * xi_block_size]();
+ vec_dir_tqsz2 = new double[ndirs * xi_block_size]();
+ vec_dir_mull = new double[16 * ndirs * xi_block_size]();
+ vec_dir_mulllr = new double[16 * ndirs * xi_block_size]();
+}
+
+InclusionOutputInfo::InclusionOutputInfo(const std::string &hdf5_file_name) {
+ unsigned int flags = H5F_ACC_RDONLY;
+ HDFFile *hdf_file = new HDFFile(hdf5_file_name, flags);
+ herr_t status = hdf_file->get_status();
+ string str_name, str_type;
+ if (status == 0) {
+ status = hdf_file->read("NSPH", "INT32_(1)", &nsph);
+ status = hdf_file->read("LI", "INT32_(1)", &li);
+ status = hdf_file->read("LE", "INT32_(1)", &le);
+ status = hdf_file->read("LM", "INT32_(1)", &lm);
+ long tmp;
+ status = hdf_file->read("MXNDM", "INT64_(1)", &tmp);
+ mxndm = (np_int)tmp;
+ status = hdf_file->read("INPOL", "INT32_(1)", &inpol);
+ status = hdf_file->read("NPNT", "INT32_(1)", &npnt);
+ status = hdf_file->read("NPNTTS", "INT32_(1)", &npntts);
+ status = hdf_file->read("IAVM", "INT32_(1)", &iavm);
+ status = hdf_file->read("ISAM", "INT32_(1)", &isam);
+ status = hdf_file->read("JWTM", "INT32_(1)", &jwtm);
+ str_type = "FLOAT64_(" + to_string(nsph) + ")";
+ vec_x_coords = new double[nsph];
+ vec_y_coords = new double[nsph];
+ vec_z_coords = new double[nsph];
+ status = hdf_file->read("VEC_SPH_X", str_type, vec_x_coords);
+ status = hdf_file->read("VEC_SPH_Y", str_type, vec_y_coords);
+ status = hdf_file->read("VEC_SPH_Z", str_type, vec_z_coords);
+ status = hdf_file->read("TH_START", "FLOAT64_(1)", &th);
+ status = hdf_file->read("TH_STEP", "FLOAT64_(1)", &thstp);
+ status = hdf_file->read("TH_END", "FLOAT64_(1)", &thlst);
+ _num_theta = (thstp == 0.0) ? 1 : 1 + int((thlst - th) / thstp);
+ status = hdf_file->read("THS_START", "FLOAT64_(1)", &ths);
+ status = hdf_file->read("THS_STEP", "FLOAT64_(1)", &thsstp);
+ status = hdf_file->read("THS_END", "FLOAT64_(1)", &thslst);
+ _num_thetas = (thsstp == 0.0) ? 1 : 1 + int((thslst - ths) / thsstp);
+ status = hdf_file->read("PH_START", "FLOAT64_(1)", &ph);
+ status = hdf_file->read("PH_STEP", "FLOAT64_(1)", &phstp);
+ status = hdf_file->read("PH_END", "FLOAT64_(1)", &phlst);
+ _num_phi = (phstp == 0.0) ? 1 : 1 + int((phlst - ph) / phstp);
+ status = hdf_file->read("PHS_START", "FLOAT64_(1)", &phs);
+ status = hdf_file->read("PHS_STEP", "FLOAT64_(1)", &phsstp);
+ status = hdf_file->read("PHS_END", "FLOAT64_(1)", &phslst);
+ _num_phis = (phsstp == 0.0) ? 1 : 1 + int((phslst - phs) / phsstp);
+ ndirs = _num_theta * _num_thetas * _num_phi * _num_phis;
+ status = hdf_file->read("EXRI", "FLOAT64_(1)", &exri);
+ status = hdf_file->read("IDFC", "INT32_(1)", &idfc);
+ status = hdf_file->read("XI1", "INT32_(1)", &_first_xi);
+ status = hdf_file->read("NXI", "INT32_(1)", &xi_block_size);
+ nxi = (_first_xi == 1) ? xi_block_size : xi_block_size + _first_xi;
+ str_type = "INT32_(" + to_string(xi_block_size) + ")";
+ vec_jxi = new int[xi_block_size];
+ status = hdf_file->read("VEC_JXI", str_type, vec_jxi);
+ str_type = "INT16_(" + to_string(xi_block_size) + ")";
+ vec_ier = new short[xi_block_size];
+ status = hdf_file->read("VEC_IER", str_type, vec_ier);
+ str_type = "FLOAT64_(" + to_string(xi_block_size) + ")";
+ vec_vk = new double[xi_block_size];
+ status = hdf_file->read("VEC_VK", str_type, vec_vk);
+ vec_xi = new double[xi_block_size];
+ status = hdf_file->read("VEC_VK", str_type, vec_xi);
+ status = hdf_file->read("NCONF", "INT32_(1)", &configurations);
+ str_type = "FLOAT64_(" + to_string(xi_block_size * (configurations + 1)) + ")";
+ vec_sphere_sizes = new double[xi_block_size * (configurations + 1)];
+ status = hdf_file->read("VEC_SPH_SIZES", str_type, vec_sphere_sizes);
+ str_type = "FLOAT64_(" + to_string(2 * xi_block_size * (configurations + 1)) + ")";
+ vec_sphere_ref_indices = new dcomplex[xi_block_size * (configurations + 1)];
+ status = hdf_file->read("VEC_SPH_REFRI", str_type, vec_sphere_ref_indices);
+ str_type = "FLOAT64_(" + to_string(xi_block_size) + ")";
+ vec_scs1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_SCS1", str_type, vec_scs1);
+ vec_scs2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_SCS2", str_type, vec_scs2);
+ vec_abs1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_ABS1", str_type, vec_abs1);
+ vec_abs2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_ABS2", str_type, vec_abs2);
+ vec_exs1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_EXS1", str_type, vec_exs1);
+ vec_exs2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_EXS2", str_type, vec_exs2);
+ vec_albeds1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_ALBEDS1", str_type, vec_albeds1);
+ vec_albeds2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_ALBEDS2", str_type, vec_albeds2);
+ vec_scsrt1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_SCSRT1", str_type, vec_scsrt1);
+ vec_scsrt2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_SCSRT2", str_type, vec_scsrt2);
+ vec_absrt1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_ABSRT1", str_type, vec_absrt1);
+ vec_absrt2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_ABSRT2", str_type, vec_absrt2);
+ vec_exsrt1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_EXSRT1", str_type, vec_exsrt1);
+ vec_exsrt2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_EXSRT2", str_type, vec_exsrt2);
+ vec_qschu1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_QSCHU1", str_type, vec_qschu1);
+ vec_qschu2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_QSCHU2", str_type, vec_qschu2);
+ vec_pschu1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_PSCHU1", str_type, vec_pschu1);
+ vec_pschu2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_PSCHU2", str_type, vec_pschu2);
+ vec_s0mag1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_S0MAG1", str_type, vec_s0mag1);
+ vec_s0mag2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_S0MAG2", str_type, vec_s0mag2);
+ vec_cosav1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_COSAV1", str_type, vec_cosav1);
+ vec_cosav2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_COSAV2", str_type, vec_cosav2);
+ vec_raprs1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_RAPRS1", str_type, vec_raprs1);
+ vec_raprs2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_RAPRS2", str_type, vec_raprs2);
+ vec_fk1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_FK1", str_type, vec_fk1);
+ vec_fk2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_FK2", str_type, vec_fk2);
+ str_type = "FLOAT64_(" + to_string(2 * xi_block_size) + ")";
+ vec_fsas11 = new dcomplex[xi_block_size];
+ status = hdf_file->read("VEC_FSAS11", str_type, vec_fsas11);
+ vec_fsas21 = new dcomplex[xi_block_size];
+ status = hdf_file->read("VEC_FSAS21", str_type, vec_fsas21);
+ vec_fsas22 = new dcomplex[xi_block_size];
+ status = hdf_file->read("VEC_FSAS22", str_type, vec_fsas22);
+ vec_fsas12 = new dcomplex[xi_block_size];
+ status = hdf_file->read("VEC_FSAS12", str_type, vec_fsas12);
+ str_type = "FLOAT64_(" + to_string(_num_theta) + ")";
+ vec_dir_tidg = new double[_num_theta];
+ status = hdf_file->read("VEC_DIR_TIDG", str_type, vec_dir_tidg);
+ str_type = "FLOAT64_(" + to_string(_num_phi) + ")";
+ vec_dir_pidg = new double[_num_phi];
+ status = hdf_file->read("VEC_DIR_PIDG", str_type, vec_dir_pidg);
+ str_type = "FLOAT64_(" + to_string(_num_thetas) + ")";
+ vec_dir_tsdg = new double[_num_thetas];
+ status = hdf_file->read("VEC_DIR_TSDG", str_type, vec_dir_tsdg);
+ str_type = "FLOAT64_(" + to_string(_num_phis) + ")";
+ vec_dir_psdg = new double[_num_phis];
+ status = hdf_file->read("VEC_DIR_PSDG", str_type, vec_dir_psdg);
+ str_type = "FLOAT64_(" + to_string(ndirs) + ")";
+ vec_dir_scand = new double[ndirs];
+ status = hdf_file->read("VEC_DIR_SCAND", str_type, vec_dir_scand);
+ vec_dir_cfmp = new double[ndirs];
+ status = hdf_file->read("VEC_DIR_CFMP", str_type, vec_dir_cfmp);
+ vec_dir_sfmp = new double[ndirs];
+ status = hdf_file->read("VEC_DIR_SFMP", str_type, vec_dir_sfmp);
+ vec_dir_cfsp = new double[ndirs];
+ status = hdf_file->read("VEC_DIR_CFSP", str_type, vec_dir_cfsp);
+ vec_dir_sfsp = new double[ndirs];
+ status = hdf_file->read("VEC_DIR_SFSP", str_type, vec_dir_sfsp);
+ str_type = "FLOAT64_(" + to_string(3 * ndirs) + ")";
+ vec_dir_un = new double[3 * ndirs];
+ status = hdf_file->read("VEC_DIR_UN", str_type, vec_dir_un);
+ vec_dir_uns = new double[3 * ndirs];
+ status = hdf_file->read("VEC_DIR_UNS", str_type, vec_dir_uns);
+ str_type = "FLOAT64_(" + to_string(ndirs * xi_block_size) + ")";
+ vec_dir_scs1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SCS1", str_type, vec_dir_scs1);
+ vec_dir_scs2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SCS2", str_type, vec_dir_scs2);
+ vec_dir_abs1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_ABS1", str_type, vec_dir_abs1);
+ vec_dir_abs2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_ABS2", str_type, vec_dir_abs2);
+ vec_dir_exs1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_EXS1", str_type, vec_dir_exs1);
+ vec_dir_exs2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_EXS2", str_type, vec_dir_exs2);
+ vec_dir_albeds1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_ALBEDS1", str_type, vec_dir_albeds1);
+ vec_dir_albeds2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_ALBEDS2", str_type, vec_dir_albeds2);
+ vec_dir_scsrt1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SCSRT1", str_type, vec_dir_scsrt1);
+ vec_dir_scsrt2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SCSRT2", str_type, vec_dir_scsrt2);
+ vec_dir_absrt1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_ABSRT1", str_type, vec_dir_absrt1);
+ vec_dir_absrt2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_ABSRT2", str_type, vec_dir_absrt2);
+ vec_dir_exsrt1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_EXSRT1", str_type, vec_dir_exsrt1);
+ vec_dir_exsrt2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_EXSRT2", str_type, vec_dir_exsrt2);
+ str_type = "FLOAT64_(" + to_string(2 * ndirs * xi_block_size) + ")";
+ vec_dir_fsas11 = new dcomplex[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FSAS11", str_type, vec_dir_fsas11);
+ vec_dir_fsas21 = new dcomplex[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FSAS21", str_type, vec_dir_fsas21);
+ vec_dir_fsas12 = new dcomplex[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FSAS12", str_type, vec_dir_fsas12);
+ vec_dir_fsas22 = new dcomplex[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FSAS22", str_type, vec_dir_fsas22);
+ vec_dir_sas11 = new dcomplex[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SAS11", str_type, vec_dir_sas11);
+ vec_dir_sas21 = new dcomplex[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SAS21", str_type, vec_dir_sas21);
+ vec_dir_sas12 = new dcomplex[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SAS12", str_type, vec_dir_sas12);
+ vec_dir_sas22 = new dcomplex[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SAS22", str_type, vec_dir_sas22);
+ str_type = "FLOAT64_(" + to_string(ndirs * xi_block_size) + ")";
+ vec_dir_qschu1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_QSCHU1", str_type, vec_dir_qschu1);
+ vec_dir_qschu2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_QSCHU2", str_type, vec_dir_qschu2);
+ vec_dir_pschu1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_PSCHU1", str_type, vec_dir_pschu1);
+ vec_dir_pschu2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_PSCHU2", str_type, vec_dir_pschu2);
+ vec_dir_s0mag1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_S0MAG1", str_type, vec_dir_s0mag1);
+ vec_dir_s0mag2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_S0MAG2", str_type, vec_dir_s0mag2);
+ vec_dir_cosav1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_COSAV1", str_type, vec_dir_cosav1);
+ vec_dir_cosav2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_COSAV2", str_type, vec_dir_cosav2);
+ vec_dir_rapr1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_RAPR1", str_type, vec_dir_rapr1);
+ vec_dir_rapr2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_RAPR2", str_type, vec_dir_rapr2);
+ vec_dir_fl1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FL1", str_type, vec_dir_fl1);
+ vec_dir_fl2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FL2", str_type, vec_dir_fl2);
+ vec_dir_fr1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FR1", str_type, vec_dir_fr1);
+ vec_dir_fr2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FR2", str_type, vec_dir_fr2);
+ vec_dir_fk1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FK1", str_type, vec_dir_fk1);
+ vec_dir_fk2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FK2", str_type, vec_dir_fk2);
+ vec_dir_fx1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FX1", str_type, vec_dir_fx1);
+ vec_dir_fx2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FX2", str_type, vec_dir_fx2);
+ vec_dir_fy1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FY1", str_type, vec_dir_fy1);
+ vec_dir_fy2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FY2", str_type, vec_dir_fy2);
+ vec_dir_fz1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FZ1", str_type, vec_dir_fz1);
+ vec_dir_fz2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FZ2", str_type, vec_dir_fz2);
+ vec_dir_tqel1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQEL1", str_type, vec_dir_tqel1);
+ vec_dir_tqel2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQEL2", str_type, vec_dir_tqel2);
+ vec_dir_tqer1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQER1", str_type, vec_dir_tqer1);
+ vec_dir_tqer2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQER2", str_type, vec_dir_tqer2);
+ vec_dir_tqek1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQEK1", str_type, vec_dir_tqek1);
+ vec_dir_tqek2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQEK2", str_type, vec_dir_tqek2);
+ vec_dir_tqex1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQEX1", str_type, vec_dir_tqex1);
+ vec_dir_tqex2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQEX2", str_type, vec_dir_tqex2);
+ vec_dir_tqey1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQEY1", str_type, vec_dir_tqey1);
+ vec_dir_tqey2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQEY2", str_type, vec_dir_tqey2);
+ vec_dir_tqez1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQEZ1", str_type, vec_dir_tqez1);
+ vec_dir_tqez2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQEZ2", str_type, vec_dir_tqez2);
+ vec_dir_tqsl1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSL1", str_type, vec_dir_tqsl1);
+ vec_dir_tqsl2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSL2", str_type, vec_dir_tqsl2);
+ vec_dir_tqsr1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSR1", str_type, vec_dir_tqsr1);
+ vec_dir_tqsr2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSR2", str_type, vec_dir_tqsr2);
+ vec_dir_tqsk1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSK1", str_type, vec_dir_tqsk1);
+ vec_dir_tqsk2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSK2", str_type, vec_dir_tqsk2);
+ vec_dir_tqsx1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSX1", str_type, vec_dir_tqsx1);
+ vec_dir_tqsx2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSX2", str_type, vec_dir_tqsx2);
+ vec_dir_tqsy1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSY1", str_type, vec_dir_tqsy1);
+ vec_dir_tqsy2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSY2", str_type, vec_dir_tqsy2);
+ vec_dir_tqsz1 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSZ1", str_type, vec_dir_tqsz1);
+ vec_dir_tqsz2 = new double[ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_TQSZ2", str_type, vec_dir_tqsz2);
+ str_type = "FLOAT64_(" + to_string(16 * ndirs * xi_block_size) + ")";
+ vec_dir_mull = new double[16 * ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_MULL", str_type, vec_dir_mull);
+ vec_dir_mulllr = new double[16 * ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_MULLLR", str_type, vec_dir_mulllr);
+ status = hdf_file->close();
+ delete hdf_file;
+ } else {
+ if (hdf_file != NULL) delete hdf_file;
+ UnrecognizedFormatException ex("Error: " + hdf5_file_name + " not recognized as a valid HDF5 file!");
+ throw ex;
+ }
+}
+
+InclusionOutputInfo::~InclusionOutputInfo() {
+ delete[] vec_x_coords;
+ delete[] vec_y_coords;
+ delete[] vec_z_coords;
+ delete[] vec_jxi;
+ delete[] vec_ier;
+ delete[] vec_vk;
+ delete[] vec_xi;
+ delete[] vec_sphere_sizes;
+ delete[] vec_sphere_ref_indices;
+ delete[] vec_scs1;
+ delete[] vec_scs2;
+ delete[] vec_abs1;
+ delete[] vec_abs2;
+ delete[] vec_exs1;
+ delete[] vec_exs2;
+ delete[] vec_albeds1;
+ delete[] vec_albeds2;
+ delete[] vec_scsrt1;
+ delete[] vec_scsrt2;
+ delete[] vec_absrt1;
+ delete[] vec_absrt2;
+ delete[] vec_exsrt1;
+ delete[] vec_exsrt2;
+ delete[] vec_qschu1;
+ delete[] vec_qschu2;
+ delete[] vec_pschu1;
+ delete[] vec_pschu2;
+ delete[] vec_s0mag1;
+ delete[] vec_s0mag2;
+ delete[] vec_cosav1;
+ delete[] vec_cosav2;
+ delete[] vec_raprs1;
+ delete[] vec_raprs2;
+ delete[] vec_fk1;
+ delete[] vec_fk2;
+ delete[] vec_fsas11;
+ delete[] vec_fsas21;
+ delete[] vec_fsas22;
+ delete[] vec_fsas12;
+ delete[] vec_dir_tidg;
+ delete[] vec_dir_pidg;
+ delete[] vec_dir_tsdg;
+ delete[] vec_dir_psdg;
+ delete[] vec_dir_scand;
+ delete[] vec_dir_cfmp;
+ delete[] vec_dir_sfmp;
+ delete[] vec_dir_cfsp;
+ delete[] vec_dir_sfsp;
+ delete[] vec_dir_un;
+ delete[] vec_dir_uns;
+ delete[] vec_dir_scs1;
+ delete[] vec_dir_scs2;
+ delete[] vec_dir_abs1;
+ delete[] vec_dir_abs2;
+ delete[] vec_dir_exs1;
+ delete[] vec_dir_exs2;
+ delete[] vec_dir_albeds1;
+ delete[] vec_dir_albeds2;
+ delete[] vec_dir_scsrt1;
+ delete[] vec_dir_scsrt2;
+ delete[] vec_dir_absrt1;
+ delete[] vec_dir_absrt2;
+ delete[] vec_dir_exsrt1;
+ delete[] vec_dir_exsrt2;
+ delete[] vec_dir_fsas11;
+ delete[] vec_dir_fsas21;
+ delete[] vec_dir_fsas12;
+ delete[] vec_dir_fsas22;
+ delete[] vec_dir_sas11;
+ delete[] vec_dir_sas21;
+ delete[] vec_dir_sas12;
+ delete[] vec_dir_sas22;
+ delete[] vec_dir_qschu1;
+ delete[] vec_dir_qschu2;
+ delete[] vec_dir_pschu1;
+ delete[] vec_dir_pschu2;
+ delete[] vec_dir_s0mag1;
+ delete[] vec_dir_s0mag2;
+ delete[] vec_dir_cosav1;
+ delete[] vec_dir_cosav2;
+ delete[] vec_dir_rapr1;
+ delete[] vec_dir_rapr2;
+ delete[] vec_dir_fl1;
+ delete[] vec_dir_fl2;
+ delete[] vec_dir_fr1;
+ delete[] vec_dir_fr2;
+ delete[] vec_dir_fk1;
+ delete[] vec_dir_fk2;
+ delete[] vec_dir_fx1;
+ delete[] vec_dir_fx2;
+ delete[] vec_dir_fy1;
+ delete[] vec_dir_fy2;
+ delete[] vec_dir_fz1;
+ delete[] vec_dir_fz2;
+ delete[] vec_dir_tqel1;
+ delete[] vec_dir_tqel2;
+ delete[] vec_dir_tqer1;
+ delete[] vec_dir_tqer2;
+ delete[] vec_dir_tqek1;
+ delete[] vec_dir_tqek2;
+ delete[] vec_dir_tqex1;
+ delete[] vec_dir_tqex2;
+ delete[] vec_dir_tqey1;
+ delete[] vec_dir_tqey2;
+ delete[] vec_dir_tqez1;
+ delete[] vec_dir_tqez2;
+ delete[] vec_dir_tqsl1;
+ delete[] vec_dir_tqsl2;
+ delete[] vec_dir_tqsr1;
+ delete[] vec_dir_tqsr2;
+ delete[] vec_dir_tqsk1;
+ delete[] vec_dir_tqsk2;
+ delete[] vec_dir_tqsx1;
+ delete[] vec_dir_tqsx2;
+ delete[] vec_dir_tqsy1;
+ delete[] vec_dir_tqsy2;
+ delete[] vec_dir_tqsz1;
+ delete[] vec_dir_tqsz2;
+ delete[] vec_dir_mull;
+ delete[] vec_dir_mulllr;
+}
+
+long InclusionOutputInfo::compute_size() {
+ long result = sizeof(np_int);
+ result += 21 * sizeof(int);
+ result += 13 * sizeof(double);
+ result += 120 * sizeof(long); // vector pointers
+ result += xi_block_size * (sizeof(int) + sizeof(short)); // vec_jxi, vec_ier
+ result += 3 * nsph * sizeof(double); // coordinate vectors
+ result += 29 * xi_block_size * sizeof(double); // wavelength dependent real values
+ result += 5 * xi_block_size * sizeof(dcomplex); // wavelength dependent complex values
+ result += 15 * ndirs * sizeof(double); // values depending only on directions
+ result += 92 * ndirs * xi_block_size * sizeof(double); // values depending on directions and wavelengths
+ result += 8 * ndirs * xi_block_size * sizeof(dcomplex); // values depending on directions and wavelengths
+ return result;
+}
+
+long InclusionOutputInfo::compute_size(
+ ScattererConfiguration *sc, GeometryConfiguration *gc,
+ int first_xi, int xi_length
+) {
+ long result = sizeof(np_int);
+ result += 21 * sizeof(int);
+ result += 13 * sizeof(double);
+ result += 120 * sizeof(long); // vector pointers
+ int _nsph = gc->number_of_spheres;
+ double _th = gc->in_theta_start;
+ double _thstp = gc->in_theta_step;
+ double _thlst = gc->in_theta_end;
+ double _ths = gc->sc_theta_start;
+ double _thsstp = gc->sc_theta_step;
+ double _thslst = gc->sc_theta_end;
+ int num_theta = 1 + int((_thlst - _th) / _thstp);
+ int num_thetas = 1 + int((_thslst - _ths) / _thsstp);
+ double _ph = gc->in_phi_start;
+ double _phstp = gc->in_phi_step;
+ double _phlst = gc->in_phi_end;
+ double _phs = gc->sc_phi_start;
+ double _phsstp = gc->sc_phi_step;
+ double _phslst = gc->sc_phi_end;
+ int num_phi = 1 + int((_phlst - _ph) / _phstp);
+ int num_phis = 1 + int((_phslst - _phs) / _phsstp);
+ int _ndirs = num_theta * num_thetas * num_phi * num_phis;
+ int _nxi = sc->number_of_scales;
+ int _xi_block_size = (xi_length == 0) ? _nxi : xi_length;
+ int _configurations = sc->configurations;
+ result += _xi_block_size * (sizeof(int) + sizeof(short)); // vec_jxi, vec_ier
+ result += 3 * _nsph * sizeof(double); // coordinate vectors
+ result += 29 * _xi_block_size * sizeof(double); // wavelength dependent real values
+ result += 5 * _xi_block_size * sizeof(dcomplex); // wavelength dependent complex values
+ result += 15 * _ndirs * sizeof(double); // values depending only on directions
+ result += 92 * _ndirs * _xi_block_size * sizeof(double); // values depending on directions and wavelengths
+ result += 8 * _ndirs * _xi_block_size * sizeof(dcomplex); // values depending on directions and wavelengths
+ return result;
+}
+
+int InclusionOutputInfo::insert(const InclusionOutputInfo &rhs) {
+ int result = 0;
+ result += (rhs.inpol == inpol) ? 0 : 1;
+ result += (rhs.iavm == iavm) ? 0 : 1;
+ result += (rhs.isam == isam) ? 0 : 1;
+ result += (rhs._num_theta == _num_theta) ? 0 : 1;
+ result += (rhs._num_thetas == _num_thetas) ? 0 : 1;
+ result += (rhs._num_phi == _num_phi) ? 0 : 1;
+ result += (rhs._num_phis == _num_phis) ? 0 : 1;
+ result += (rhs.ndirs == ndirs) ? 0 : 1;
+ result += (rhs.exri == exri) ? 0 : 1;
+ result += (rhs.idfc == idfc) ? 0 : 1;
+ result += (rhs.configurations == configurations) ? 0 : 1;
+ if (result == 0) {
+ int offset, chunk_size, xi1;
+ xi1 = rhs._first_xi;
+ // Insert vectors whose size depends on wavelengths
+ offset = xi1 - _first_xi;
+ chunk_size = rhs.xi_block_size;
+ memcpy(vec_jxi + offset, rhs.vec_jxi, chunk_size * sizeof(int));
+ memcpy(vec_ier + offset, rhs.vec_ier, chunk_size * sizeof(short));
+ memcpy(vec_vk + offset, rhs.vec_vk, chunk_size * sizeof(double));
+ memcpy(vec_xi + offset, rhs.vec_xi, chunk_size * sizeof(double));
+ memcpy(vec_scs1 + offset, rhs.vec_scs1, chunk_size * sizeof(double));
+ memcpy(vec_scs2 + offset, rhs.vec_scs2, chunk_size * sizeof(double));
+ memcpy(vec_abs1 + offset, rhs.vec_abs1, chunk_size * sizeof(double));
+ memcpy(vec_abs2 + offset, rhs.vec_abs2, chunk_size * sizeof(double));
+ memcpy(vec_exs1 + offset, rhs.vec_exs1, chunk_size * sizeof(double));
+ memcpy(vec_exs2 + offset, rhs.vec_exs2, chunk_size * sizeof(double));
+ memcpy(vec_albeds1 + offset, rhs.vec_albeds1, chunk_size * sizeof(double));
+ memcpy(vec_albeds2 + offset, rhs.vec_albeds2, chunk_size * sizeof(double));
+ memcpy(vec_scsrt1 + offset, rhs.vec_scsrt1, chunk_size * sizeof(double));
+ memcpy(vec_scsrt2 + offset, rhs.vec_scsrt2, chunk_size * sizeof(double));
+ memcpy(vec_absrt1 + offset, rhs.vec_absrt1, chunk_size * sizeof(double));
+ memcpy(vec_absrt2 + offset, rhs.vec_absrt2, chunk_size * sizeof(double));
+ memcpy(vec_exsrt1 + offset, rhs.vec_exsrt1, chunk_size * sizeof(double));
+ memcpy(vec_exsrt2 + offset, rhs.vec_exsrt2, chunk_size * sizeof(double));
+ memcpy(vec_qschu1 + offset, rhs.vec_qschu1, chunk_size * sizeof(double));
+ memcpy(vec_qschu2 + offset, rhs.vec_qschu2, chunk_size * sizeof(double));
+ memcpy(vec_pschu1 + offset, rhs.vec_pschu1, chunk_size * sizeof(double));
+ memcpy(vec_pschu2 + offset, rhs.vec_pschu2, chunk_size * sizeof(double));
+ memcpy(vec_s0mag1 + offset, rhs.vec_s0mag1, chunk_size * sizeof(double));
+ memcpy(vec_s0mag2 + offset, rhs.vec_s0mag2, chunk_size * sizeof(double));
+ memcpy(vec_cosav1 + offset, rhs.vec_cosav1, chunk_size * sizeof(double));
+ memcpy(vec_cosav2 + offset, rhs.vec_cosav2, chunk_size * sizeof(double));
+ memcpy(vec_raprs1 + offset, rhs.vec_raprs1, chunk_size * sizeof(double));
+ memcpy(vec_raprs2 + offset, rhs.vec_raprs2, chunk_size * sizeof(double));
+ memcpy(vec_fk1 + offset, rhs.vec_fk1, chunk_size * sizeof(double));
+ memcpy(vec_fk2 + offset, rhs.vec_fk2, chunk_size * sizeof(double));
+ memcpy(vec_fsas11 + offset, rhs.vec_fsas11, chunk_size * sizeof(dcomplex));
+ memcpy(vec_fsas21 + offset, rhs.vec_fsas21, chunk_size * sizeof(dcomplex));
+ memcpy(vec_fsas22 + offset, rhs.vec_fsas22, chunk_size * sizeof(dcomplex));
+ memcpy(vec_fsas12 + offset, rhs.vec_fsas12, chunk_size * sizeof(dcomplex));
+
+ // Insert vectors whose sizes depend on configurations
+ offset = (xi1 - _first_xi) * (configurations + 1);
+ chunk_size = rhs.xi_block_size * (configurations + 1);
+ memcpy(vec_sphere_sizes + offset, rhs.vec_sphere_sizes, chunk_size * sizeof(double));
+ memcpy(vec_sphere_ref_indices + offset, rhs.vec_sphere_ref_indices, chunk_size * sizeof(dcomplex));
+
+ // Insert vectors whose sizes depend on directions and wavelengths
+ offset = (xi1 - _first_xi) * ndirs;
+ chunk_size = rhs.xi_block_size * ndirs;
+ memcpy(vec_dir_scs1 + offset, rhs.vec_dir_scs1, chunk_size * sizeof(double));
+ memcpy(vec_dir_scs2 + offset, rhs.vec_dir_scs2, chunk_size * sizeof(double));
+ memcpy(vec_dir_abs1 + offset, rhs.vec_dir_abs1, chunk_size * sizeof(double));
+ memcpy(vec_dir_abs2 + offset, rhs.vec_dir_abs2, chunk_size * sizeof(double));
+ memcpy(vec_dir_exs1 + offset, rhs.vec_dir_exs1, chunk_size * sizeof(double));
+ memcpy(vec_dir_exs2 + offset, rhs.vec_dir_exs2, chunk_size * sizeof(double));
+ memcpy(vec_dir_albeds1 + offset, rhs.vec_dir_albeds1, chunk_size * sizeof(double));
+ memcpy(vec_dir_albeds2 + offset, rhs.vec_dir_albeds2, chunk_size * sizeof(double));
+ memcpy(vec_dir_scsrt1 + offset, rhs.vec_dir_scsrt1, chunk_size * sizeof(double));
+ memcpy(vec_dir_scsrt2 + offset, rhs.vec_dir_scsrt2, chunk_size * sizeof(double));
+ memcpy(vec_dir_absrt1 + offset, rhs.vec_dir_absrt1, chunk_size * sizeof(double));
+ memcpy(vec_dir_absrt2 + offset, rhs.vec_dir_absrt2, chunk_size * sizeof(double));
+ memcpy(vec_dir_exsrt1 + offset, rhs.vec_dir_exsrt1, chunk_size * sizeof(double));
+ memcpy(vec_dir_exsrt2 + offset, rhs.vec_dir_exsrt2, chunk_size * sizeof(double));
+ memcpy(vec_dir_fsas11 + offset, rhs.vec_dir_fsas11, chunk_size * sizeof(dcomplex));
+ memcpy(vec_dir_fsas21 + offset, rhs.vec_dir_fsas21, chunk_size * sizeof(dcomplex));
+ memcpy(vec_dir_fsas12 + offset, rhs.vec_dir_fsas12, chunk_size * sizeof(dcomplex));
+ memcpy(vec_dir_fsas22 + offset, rhs.vec_dir_fsas22, chunk_size * sizeof(dcomplex));
+ memcpy(vec_dir_sas11 + offset, rhs.vec_dir_sas11, chunk_size * sizeof(dcomplex));
+ memcpy(vec_dir_sas21 + offset, rhs.vec_dir_sas21, chunk_size * sizeof(dcomplex));
+ memcpy(vec_dir_sas12 + offset, rhs.vec_dir_sas12, chunk_size * sizeof(dcomplex));
+ memcpy(vec_dir_sas22 + offset, rhs.vec_dir_sas22, chunk_size * sizeof(dcomplex));
+ memcpy(vec_dir_qschu1 + offset, rhs.vec_dir_qschu1, chunk_size * sizeof(double));
+ memcpy(vec_dir_qschu2 + offset, rhs.vec_dir_qschu2, chunk_size * sizeof(double));
+ memcpy(vec_dir_pschu1 + offset, rhs.vec_dir_pschu1, chunk_size * sizeof(double));
+ memcpy(vec_dir_pschu2 + offset, rhs.vec_dir_pschu2, chunk_size * sizeof(double));
+ memcpy(vec_dir_s0mag1 + offset, rhs.vec_dir_s0mag1, chunk_size * sizeof(double));
+ memcpy(vec_dir_s0mag2 + offset, rhs.vec_dir_s0mag2, chunk_size * sizeof(double));
+ memcpy(vec_dir_cosav1 + offset, rhs.vec_dir_cosav1, chunk_size * sizeof(double));
+ memcpy(vec_dir_cosav2 + offset, rhs.vec_dir_cosav2, chunk_size * sizeof(double));
+ memcpy(vec_dir_rapr1 + offset, rhs.vec_dir_rapr1, chunk_size * sizeof(double));
+ memcpy(vec_dir_rapr2 + offset, rhs.vec_dir_rapr2, chunk_size * sizeof(double));
+ memcpy(vec_dir_fl1 + offset, rhs.vec_dir_fl1, chunk_size * sizeof(double));
+ memcpy(vec_dir_fl2 + offset, rhs.vec_dir_fl2, chunk_size * sizeof(double));
+ memcpy(vec_dir_fr1 + offset, rhs.vec_dir_fr1, chunk_size * sizeof(double));
+ memcpy(vec_dir_fr2 + offset, rhs.vec_dir_fr2, chunk_size * sizeof(double));
+ memcpy(vec_dir_fk1 + offset, rhs.vec_dir_fk1, chunk_size * sizeof(double));
+ memcpy(vec_dir_fk2 + offset, rhs.vec_dir_fk2, chunk_size * sizeof(double));
+ memcpy(vec_dir_fx1 + offset, rhs.vec_dir_fx1, chunk_size * sizeof(double));
+ memcpy(vec_dir_fx2 + offset, rhs.vec_dir_fx2, chunk_size * sizeof(double));
+ memcpy(vec_dir_fy1 + offset, rhs.vec_dir_fy1, chunk_size * sizeof(double));
+ memcpy(vec_dir_fy2 + offset, rhs.vec_dir_fy2, chunk_size * sizeof(double));
+ memcpy(vec_dir_fz1 + offset, rhs.vec_dir_fz1, chunk_size * sizeof(double));
+ memcpy(vec_dir_fz2 + offset, rhs.vec_dir_fz2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqel1 + offset, rhs.vec_dir_tqel1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqel2 + offset, rhs.vec_dir_tqel2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqer1 + offset, rhs.vec_dir_tqer1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqer2 + offset, rhs.vec_dir_tqer2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqek1 + offset, rhs.vec_dir_tqek1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqek2 + offset, rhs.vec_dir_tqek2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqex1 + offset, rhs.vec_dir_tqex1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqex2 + offset, rhs.vec_dir_tqex2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqey1 + offset, rhs.vec_dir_tqey1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqey2 + offset, rhs.vec_dir_tqey2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqez1 + offset, rhs.vec_dir_tqez1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqez2 + offset, rhs.vec_dir_tqez2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsl1 + offset, rhs.vec_dir_tqsl1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsl2 + offset, rhs.vec_dir_tqsl2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsr1 + offset, rhs.vec_dir_tqsr1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsr2 + offset, rhs.vec_dir_tqsr2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsk1 + offset, rhs.vec_dir_tqsk1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsk2 + offset, rhs.vec_dir_tqsk2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsx1 + offset, rhs.vec_dir_tqsx1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsx2 + offset, rhs.vec_dir_tqsx2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsy1 + offset, rhs.vec_dir_tqsy1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsy2 + offset, rhs.vec_dir_tqsy2, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsz1 + offset, rhs.vec_dir_tqsz1, chunk_size * sizeof(double));
+ memcpy(vec_dir_tqsz2 + offset, rhs.vec_dir_tqsz2, chunk_size * sizeof(double));
+ memcpy(vec_dir_mull + 16 * offset, rhs.vec_dir_mull, 16 * chunk_size * sizeof(double));
+ memcpy(vec_dir_mulllr + 16 * offset, rhs.vec_dir_mulllr, 16 * chunk_size * sizeof(double));
+ }
+ return result;
+}
+
+int InclusionOutputInfo::write(const std::string &output, const std::string &format) {
+ int result = 0;
+ if (format.compare("LEGACY") == 0) {
+ result = write_legacy(output);
+ } else if (format.compare("HDF5") == 0) {
+ result = write_hdf5(output);
+ } else {
+ string message = "Unknown format mode: \"" + format + "\"";
+ throw UnrecognizedConfigurationException(message);
+ }
+ return result;
+}
+
+int InclusionOutputInfo::write_hdf5(const std::string &file_name) {
+ List<string> *rec_name_list = new List<string>(1);
+ List<string> *rec_type_list = new List<string>(1);
+ List<void *> *rec_ptr_list = new List<void *>(1);
+ string str_type, str_name;
+ rec_name_list->set(0, "NSPH");
+ rec_type_list->set(0, "INT32_(1)");
+ rec_ptr_list->set(0, &nsph);
+ rec_name_list->append("LI");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&li);
+ rec_name_list->append("LE");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&le);
+ rec_name_list->append("LM");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&lm);
+ rec_name_list->append("MXNDM");
+ rec_type_list->append("INT64_(1)");
+ rec_ptr_list->append(&mxndm);
+ rec_name_list->append("INPOL");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&inpol);
+ rec_name_list->append("NPNT");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&npnt);
+ rec_name_list->append("NPNTTS");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&npntts);
+ rec_name_list->append("IAVM");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&iavm);
+ rec_name_list->append("ISAM");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&isam);
+ rec_name_list->append("JWTM");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&jwtm);
+ rec_name_list->append("VEC_SPH_X");
+ rec_type_list->append("FLOAT64_(" + to_string(nsph) + ")");
+ rec_ptr_list->append(vec_x_coords);
+ rec_name_list->append("VEC_SPH_Y");
+ rec_type_list->append("FLOAT64_(" + to_string(nsph) + ")");
+ rec_ptr_list->append(vec_y_coords);
+ rec_name_list->append("VEC_SPH_Z");
+ rec_type_list->append("FLOAT64_(" + to_string(nsph) + ")");
+ rec_ptr_list->append(vec_z_coords);
+ rec_name_list->append("TH_START");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&th);
+ rec_name_list->append("TH_STEP");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&thstp);
+ rec_name_list->append("TH_END");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&thlst);
+ rec_name_list->append("THS_START");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&ths);
+ rec_name_list->append("THS_STEP");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&thsstp);
+ rec_name_list->append("THS_END");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&thslst);
+ rec_name_list->append("PH_START");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&ph);
+ rec_name_list->append("PH_STEP");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&phstp);
+ rec_name_list->append("PH_END");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&phlst);
+ rec_name_list->append("PHS_START");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&phs);
+ rec_name_list->append("PHS_STEP");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&phsstp);
+ rec_name_list->append("PHS_END");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&phslst);
+ rec_name_list->append("EXRI");
+ rec_type_list->append("FLOAT64_(1)");
+ rec_ptr_list->append(&exri);
+ rec_name_list->append("IDFC");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&idfc);
+ rec_name_list->append("XI1");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&_first_xi);
+ rec_name_list->append("NXI");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&xi_block_size);
+ rec_name_list->append("VEC_JXI");
+ rec_type_list->append("INT32_(" + to_string(xi_block_size) + ")");
+ rec_ptr_list->append(vec_jxi);
+ rec_name_list->append("VEC_IER");
+ rec_type_list->append("INT16_(" + to_string(xi_block_size) + ")");
+ rec_ptr_list->append(vec_ier);
+ rec_name_list->append("VEC_VK");
+ rec_type_list->append("FLOAT64_(" + to_string(xi_block_size) + ")");
+ rec_ptr_list->append(vec_vk);
+ rec_name_list->append("VEC_XI");
+ rec_type_list->append("FLOAT64_(" + to_string(xi_block_size) + ")");
+ rec_ptr_list->append(vec_xi);
+ rec_name_list->append("NCONF");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&configurations);
+ rec_name_list->append("VEC_SPH_SIZES");
+ rec_type_list->append("FLOAT64_(" + to_string(xi_block_size * (configurations + 1)) + ")");
+ rec_ptr_list->append(vec_sphere_sizes);
+ rec_name_list->append("VEC_SPH_REFRI");
+ rec_type_list->append("FLOAT64_(" + to_string(2 * xi_block_size * (configurations+ 1)) + ")");
+ rec_ptr_list->append(vec_sphere_ref_indices);
+ str_type = "FLOAT64_(" + to_string(xi_block_size) + ")";
+ rec_name_list->append("VEC_SCS1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_scs1);
+ rec_name_list->append("VEC_SCS2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_scs2);
+ rec_name_list->append("VEC_ABS1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_abs1);
+ rec_name_list->append("VEC_ABS2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_abs2);
+ rec_name_list->append("VEC_EXS1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_exs1);
+ rec_name_list->append("VEC_EXS2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_exs2);
+ rec_name_list->append("VEC_ALBEDS1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_albeds1);
+ rec_name_list->append("VEC_ALBEDS2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_albeds2);
+ rec_name_list->append("VEC_SCSRT1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_scsrt1);
+ rec_name_list->append("VEC_SCSRT2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_scsrt2);
+ rec_name_list->append("VEC_ABSRT1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_absrt1);
+ rec_name_list->append("VEC_ABSRT2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_absrt2);
+ rec_name_list->append("VEC_EXSRT1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_exsrt1);
+ rec_name_list->append("VEC_EXSRT2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_exsrt2);
+ rec_name_list->append("VEC_QSCHU1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_qschu1);
+ rec_name_list->append("VEC_QSCHU2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_qschu2);
+ rec_name_list->append("VEC_PSCHU1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_pschu1);
+ rec_name_list->append("VEC_PSCHU2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_pschu2);
+ rec_name_list->append("VEC_S0MAG1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_s0mag1);
+ rec_name_list->append("VEC_S0MAG2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_s0mag2);
+ rec_name_list->append("VEC_COSAV1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_cosav1);
+ rec_name_list->append("VEC_COSAV2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_cosav2);
+ rec_name_list->append("VEC_RAPRS1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_raprs1);
+ rec_name_list->append("VEC_RAPRS2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_raprs2);
+ rec_name_list->append("VEC_FK1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_fk1);
+ rec_name_list->append("VEC_FK2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_fk2);
+ str_type = "FLOAT64_(" + to_string(2 * xi_block_size) + ")";
+ rec_name_list->append("VEC_FSAS11");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_fsas11);
+ rec_name_list->append("VEC_FSAS21");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_fsas21);
+ rec_name_list->append("VEC_FSAS12");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_fsas12);
+ rec_name_list->append("VEC_FSAS22");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_fsas22);
+ str_type = "FLOAT64_(" + to_string(_num_theta) + ")";
+ rec_name_list->append("VEC_DIR_TIDG");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tidg);
+ str_type = "FLOAT64_(" + to_string(_num_phi) + ")";
+ rec_name_list->append("VEC_DIR_PIDG");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_pidg);
+ str_type = "FLOAT64_(" + to_string(_num_thetas) + ")";
+ rec_name_list->append("VEC_DIR_TSDG");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tsdg);
+ str_type = "FLOAT64_(" + to_string(_num_phis) + ")";
+ rec_name_list->append("VEC_DIR_PSDG");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_psdg);
+ str_type = "FLOAT64_(" + to_string(ndirs) + ")";
+ rec_name_list->append("VEC_DIR_SCAND");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_scand);
+ rec_name_list->append("VEC_DIR_CFMP");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_cfmp);
+ rec_name_list->append("VEC_DIR_CFSP");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_cfsp);
+ rec_name_list->append("VEC_DIR_SFMP");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_sfmp);
+ rec_name_list->append("VEC_DIR_SFSP");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_sfsp);
+ str_type = "FLOAT64_(" + to_string(3 * ndirs) + ")";
+ rec_name_list->append("VEC_DIR_UN");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_un);
+ rec_name_list->append("VEC_DIR_UNS");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_uns);
+ str_type = "FLOAT64_(" + to_string(ndirs * xi_block_size) + ")";
+ rec_name_list->append("VEC_DIR_SCS1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_scs1);
+ rec_name_list->append("VEC_DIR_SCS2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_scs2);
+ rec_name_list->append("VEC_DIR_ABS1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_abs1);
+ rec_name_list->append("VEC_DIR_ABS2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_abs2);
+ rec_name_list->append("VEC_DIR_EXS1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_exs1);
+ rec_name_list->append("VEC_DIR_EXS2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_exs2);
+ rec_name_list->append("VEC_DIR_ALBEDS1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_albeds1);
+ rec_name_list->append("VEC_DIR_ALBEDS2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_albeds2);
+ rec_name_list->append("VEC_DIR_SCSRT1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_scsrt1);
+ rec_name_list->append("VEC_DIR_SCSRT2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_scsrt2);
+ rec_name_list->append("VEC_DIR_ABSRT1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_absrt1);
+ rec_name_list->append("VEC_DIR_ABSRT2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_absrt2);
+ rec_name_list->append("VEC_DIR_EXSRT1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_exsrt1);
+ rec_name_list->append("VEC_DIR_EXSRT2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_exsrt2);
+ str_type = "FLOAT64_(" + to_string(2 * ndirs * xi_block_size) + ")";
+ rec_name_list->append("VEC_DIR_FSAS11");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fsas11);
+ rec_name_list->append("VEC_DIR_FSAS21");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fsas21);
+ rec_name_list->append("VEC_DIR_FSAS12");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fsas12);
+ rec_name_list->append("VEC_DIR_FSAS22");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fsas22);
+ rec_name_list->append("VEC_DIR_SAS11");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_sas11);
+ rec_name_list->append("VEC_DIR_SAS21");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_sas21);
+ rec_name_list->append("VEC_DIR_SAS12");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_sas12);
+ rec_name_list->append("VEC_DIR_SAS22");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_sas22);
+ str_type = "FLOAT64_(" + to_string(ndirs * xi_block_size) + ")";
+ rec_name_list->append("VEC_DIR_QSCHU1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_qschu1);
+ rec_name_list->append("VEC_DIR_QSCHU2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_qschu2);
+ rec_name_list->append("VEC_DIR_PSCHU1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_pschu1);
+ rec_name_list->append("VEC_DIR_PSCHU2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_pschu2);
+ rec_name_list->append("VEC_DIR_S0MAG1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_s0mag1);
+ rec_name_list->append("VEC_DIR_S0MAG2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_s0mag2);
+ rec_name_list->append("VEC_DIR_COSAV1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_cosav1);
+ rec_name_list->append("VEC_DIR_COSAV2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_cosav2);
+ rec_name_list->append("VEC_DIR_RAPR1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_rapr1);
+ rec_name_list->append("VEC_DIR_RAPR2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_rapr2);
+ rec_name_list->append("VEC_DIR_FL1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fl1);
+ rec_name_list->append("VEC_DIR_FL2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fl2);
+ rec_name_list->append("VEC_DIR_FR1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fr1);
+ rec_name_list->append("VEC_DIR_FR2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fr2);
+ rec_name_list->append("VEC_DIR_FK1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fk1);
+ rec_name_list->append("VEC_DIR_FK2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fk2);
+ rec_name_list->append("VEC_DIR_FX1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fx1);
+ rec_name_list->append("VEC_DIR_FX2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fx2);
+ rec_name_list->append("VEC_DIR_FY1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fy1);
+ rec_name_list->append("VEC_DIR_FY2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fy2);
+ rec_name_list->append("VEC_DIR_FZ1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fz1);
+ rec_name_list->append("VEC_DIR_FZ2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fz2);
+ rec_name_list->append("VEC_DIR_TQEL1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqel1);
+ rec_name_list->append("VEC_DIR_TQEL2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqel2);
+ rec_name_list->append("VEC_DIR_TQER1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqer1);
+ rec_name_list->append("VEC_DIR_TQER2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqer2);
+ rec_name_list->append("VEC_DIR_TQEK1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqek1);
+ rec_name_list->append("VEC_DIR_TQEK2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqek2);
+ rec_name_list->append("VEC_DIR_TQEX1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqex1);
+ rec_name_list->append("VEC_DIR_TQEX2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqex2);
+ rec_name_list->append("VEC_DIR_TQEY1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqey1);
+ rec_name_list->append("VEC_DIR_TQEY2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqey2);
+ rec_name_list->append("VEC_DIR_TQEZ1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqez1);
+ rec_name_list->append("VEC_DIR_TQEZ2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqez2);
+ rec_name_list->append("VEC_DIR_TQSL1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsl1);
+ rec_name_list->append("VEC_DIR_TQSL2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsl2);
+ rec_name_list->append("VEC_DIR_TQSR1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsr1);
+ rec_name_list->append("VEC_DIR_TQSR2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsr2);
+ rec_name_list->append("VEC_DIR_TQSK1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsk1);
+ rec_name_list->append("VEC_DIR_TQSK2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsk2);
+ rec_name_list->append("VEC_DIR_TQSX1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsx1);
+ rec_name_list->append("VEC_DIR_TQSX2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsx2);
+ rec_name_list->append("VEC_DIR_TQSY1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsy1);
+ rec_name_list->append("VEC_DIR_TQSY2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsy2);
+ rec_name_list->append("VEC_DIR_TQSZ1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsz1);
+ rec_name_list->append("VEC_DIR_TQSZ2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_tqsz2);
+ str_type = "FLOAT64_(" + to_string(16 * ndirs * xi_block_size) + ")";
+ rec_name_list->append("VEC_DIR_MULL");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_mull);
+ rec_name_list->append("VEC_DIR_MULLLR");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_mulllr);
+
+ // Convert the lists to arrays and write them to HDF5
+ string *rec_names = rec_name_list->to_array();
+ string *rec_types = rec_type_list->to_array();
+ void **rec_pointers = rec_ptr_list->to_array();
+ const int rec_num = rec_name_list->length();
+ FileSchema *schema = new FileSchema(rec_num, rec_types, rec_names);
+ HDFFile *hdf_file = HDFFile::from_schema(*schema, file_name, H5F_ACC_TRUNC);
+ for (int ri = 0; ri < rec_num; ri++)
+ hdf_file->write(rec_names[ri], rec_types[ri], rec_pointers[ri]);
+ hdf_file->close();
+
+ // Clean memory
+ delete rec_name_list;
+ delete rec_type_list;
+ delete rec_ptr_list;
+ delete[] rec_names;
+ delete[] rec_types;
+ delete[] rec_pointers;
+ delete schema;
+ delete hdf_file;
+ return 0;
+}
+
+int InclusionOutputInfo::write_legacy(const std::string &output) {
+ const dcomplex cc0 = 0.0 + I * 0.0;
+ int result = 0;
+ FILE *p_outfile = fopen(output.c_str(), "w");
+ if (p_outfile != NULL) {
+ if (vec_jxi[0] == 1) {
+ fprintf(p_outfile, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
+#ifdef USE_ILP64
+ fprintf(
+ p_outfile, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
+ nsph, li, le, mxndm, inpol, npnt, npntts,
+ iavm, iavm
+ );
+#else
+ fprintf(
+ p_outfile, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
+ nsph, li, le, mxndm, inpol, npnt, npntts,
+ iavm, iavm
+ );
+#endif // USE_ILP64
+ fprintf(p_outfile, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
+ for (int ri = 0; ri < nsph; ri++) {
+ fprintf(
+ p_outfile, "%17.8lE%17.8lE%17.8lE\n",
+ vec_x_coords[ri], vec_y_coords[ri], vec_z_coords[ri]
+ );
+ }
+ fprintf(p_outfile, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
+ fprintf(
+ p_outfile, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
+ th, thstp, thlst, ths, thsstp, thslst
+ );
+ fprintf(p_outfile, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
+ fprintf(
+ p_outfile, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
+ ph, phstp, phlst, phs, phsstp, phslst
+ );
+ fprintf(p_outfile, " READ(IR,*)JWTM\n");
+ fprintf(p_outfile, " %5d\n", jwtm);
+ fprintf(p_outfile, " READ(ITIN)NSPHT\n");
+ fprintf(p_outfile, " READ(ITIN)(IOG(I),I=1,NSPH)\n");
+ fprintf(p_outfile, " READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
+ fprintf(p_outfile, " READ(ITIN)(XIV(I),I=1,NXI)\n");
+ fprintf(p_outfile, " READ(ITIN)NSHL(I),ROS(I)\n");
+ fprintf(p_outfile, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n \n");
+ fprintf(p_outfile, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
+ if (idfc < 0) {
+ fprintf( p_outfile, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n \n", vec_vk[0]);
+ }
+ // End preamble writing
+ }
+ // Wavelength loop
+ for (int jxi = 0; jxi < xi_block_size; jxi++) {
+ int done_dirs = 0;
+ double alamb = 2.0 * 3.141592653589793238 / vec_vk[jxi];
+ fprintf(p_outfile, "========== JXI =%3d ====================\n", vec_jxi[jxi]);
+ if (idfc >= 0) {
+ fprintf(p_outfile, " VK=%15.7lE, XI=%15.7lE\n", vec_vk[jxi], vec_xi[jxi]);
+ } else {
+ fprintf(p_outfile, " XI=%15.7lE\n", vec_xi[jxi]);
+ }
+ if (vec_ier[jxi] == 1) {
+ fprintf(p_outfile, " STOP IN INDME\n");
+ break;
+ } else if (vec_ier[jxi] == 2) {
+ fprintf(p_outfile, " STOP IN OSPV\n");
+ break;
+ }
+ for (int i168 = 1; i168 <= configurations; i168++) {
+ int cindex = jxi * (configurations + 1) + i168 - 1;
+ if (vec_sphere_ref_indices[cindex] == cc0) {
+ fprintf(p_outfile, " SPHERE N.%2d: SIZE=%15.7lE\n", i168, vec_sphere_sizes[cindex]);
+ } else {
+ fprintf(
+ p_outfile, " SPHERE N.%2d: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
+ i168, vec_sphere_sizes[cindex], real(vec_sphere_ref_indices[cindex]),
+ imag(vec_sphere_ref_indices[cindex])
+ );
+ }
+ } // i168 configuration loop
+ int cindex = jxi * (configurations + 1) + configurations;
+ fprintf(
+ p_outfile, " EXT. SPHERE: SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
+ vec_sphere_sizes[cindex], real(vec_sphere_ref_indices[cindex]),
+ imag(vec_sphere_ref_indices[cindex])
+ );
+ fprintf(p_outfile, " ENSEMBLE AVERAGE, MODE%2d\n", iavm);
+ if (inpol == 0) fprintf(p_outfile, " LIN -1\n");
+ else fprintf(p_outfile, " CIRC -1\n");
+ fprintf(p_outfile, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+ fprintf(
+ p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_scs1[jxi], vec_abs1[jxi], vec_exs1[jxi], vec_albeds1[jxi]
+ );
+ fprintf(p_outfile, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n",
+ -1, alamb, vec_scs1[jxi], vec_abs1[jxi], vec_exs1[jxi]
+ );
+ fprintf(p_outfile, " ---- SCS/GS -- ABC/GS -- EXS/GS ---\n");
+ fprintf(
+ p_outfile, " %14.7lE%15.7lE%15.7lE\n",
+ vec_scsrt1[jxi], vec_absrt1[jxi], vec_exsrt1[jxi]
+ );
+ fprintf(
+ p_outfile, " FSAS(1,1)=%15.7lE%15.7lE FSAS(2,1)=%15.7lE%15.7lE\n",
+ real(vec_fsas11[jxi]), imag(vec_fsas11[jxi]),
+ real(vec_fsas21[jxi]), imag(vec_fsas21[jxi])
+ );
+ fprintf(
+ p_outfile, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ vec_qschu1[jxi], vec_pschu1[jxi], vec_s0mag1[jxi]
+ );
+ fprintf(
+ p_outfile, " COSAV=%15.7lE, RAPRS=%15.7lE\n",
+ vec_cosav1[jxi], vec_raprs1[jxi]
+ );
+ fprintf(p_outfile, " Fk=%15.7lE\n", vec_fk1[jxi]);
+ if (inpol == 0) fprintf(p_outfile, " LIN 1\n");
+ else fprintf(p_outfile, " CIRC 1\n");
+ fprintf(p_outfile, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+ fprintf(
+ p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_scs2[jxi], vec_abs2[jxi], vec_exs2[jxi], vec_albeds2[jxi]
+ );
+ fprintf(p_outfile, "INSERTION: SCASECM %5d%15.7E%15.7E%15.7E%15.7E\n",
+ 1, alamb, vec_scs2[jxi], vec_abs2[jxi], vec_exs2[jxi]
+ );
+ fprintf(p_outfile, " ---- SCS/GS -- ABC/GS -- EXS/GS ---\n");
+ fprintf(
+ p_outfile, " %14.7lE%15.7lE%15.7lE\n",
+ vec_scsrt2[jxi], vec_absrt2[jxi], vec_exsrt2[jxi]
+ );
+ fprintf(
+ p_outfile, " FSAS(2,2)=%15.7lE%15.7lE FSAS(1,2)=%15.7lE%15.7lE\n",
+ real(vec_fsas22[jxi]), imag(vec_fsas22[jxi]),
+ real(vec_fsas12[jxi]), imag(vec_fsas12[jxi])
+ );
+ fprintf(
+ p_outfile, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ vec_qschu2[jxi], vec_pschu2[jxi], vec_s0mag2[jxi]
+ );
+ fprintf(
+ p_outfile, " COSAV=%15.7lE, RAPRS=%15.7lE\n",
+ vec_cosav2[jxi], vec_raprs2[jxi]
+ );
+ fprintf(p_outfile, " Fk=%15.7lE\n", vec_fk2[jxi]);
+ fprintf(
+ p_outfile, " (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n",
+ (real(vec_fsas11[jxi]) - real(vec_fsas22[jxi])) / real(vec_fsas11[jxi])
+ );
+ fprintf(
+ p_outfile, " (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n",
+ (imag(vec_fsas11[jxi]) - imag(vec_fsas22[jxi])) / imag(vec_fsas11[jxi])
+ );
+ for (int jth = 0; jth < _num_theta; jth++) {
+ for (int jph = 0; jph < _num_phi; jph++) {
+ for (int jths = 0; jths < _num_thetas; jths++) {
+ for (int jphs = 0; jphs < _num_phis; jphs++) {
+ bool goto290 = isam >= 0 && (jths > 0 || jphs > 0);
+ int dir_index = ndirs * jxi + done_dirs;
+ fprintf(
+ p_outfile, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n",
+ jth + 1, jph + 1, jths + 1, jphs +1
+ );
+ fprintf(
+ p_outfile, " TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n",
+ th + jth * thstp, ph + jph * phstp, ths + jths * thsstp, phs + jphs * phsstp
+ );
+ fprintf(p_outfile, " SCAND=%10.3lE\n", vec_dir_scand[done_dirs]);
+ fprintf(
+ p_outfile, " CFMP=%15.7lE, SFMP=%15.7lE\n",
+ vec_dir_cfmp[done_dirs], vec_dir_sfmp[done_dirs]
+ );
+ fprintf(
+ p_outfile, " CFSP=%15.7lE, SFSP=%15.7lE\n",
+ vec_dir_cfsp[done_dirs], vec_dir_sfsp[done_dirs]
+ );
+ if (isam >= 0) {
+ fprintf(
+ p_outfile, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n",
+ vec_dir_un[done_dirs],
+ vec_dir_un[done_dirs + 1],
+ vec_dir_un[done_dirs + 2]
+ );
+ fprintf(
+ p_outfile, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n",
+ vec_dir_uns[done_dirs],
+ vec_dir_uns[done_dirs + 1],
+ vec_dir_uns[done_dirs + 2]
+ );
+ } else { // label 214
+ fprintf(
+ p_outfile, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n",
+ vec_dir_un[done_dirs],
+ vec_dir_un[done_dirs + 1],
+ vec_dir_un[done_dirs + 2]
+ );
+ }
+ fprintf(p_outfile, " SINGLE SCATTERER\n");
+ if (inpol == 0) {
+ fprintf(p_outfile, " LIN -1\n");
+ } else {
+ fprintf(p_outfile, " CIRC -1\n");
+ }
+ // label 275
+ fprintf(p_outfile, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+ fprintf(
+ p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_dir_scs1[dir_index],
+ vec_dir_abs1[dir_index],
+ vec_dir_exs1[dir_index],
+ vec_dir_albeds1[dir_index]
+ );
+ fprintf(
+ p_outfile, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n",
+ -1, alamb, vec_dir_scs1[dir_index], vec_dir_abs1[dir_index],
+ vec_dir_exs1[dir_index]
+ );
+ fprintf(p_outfile, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
+ fprintf(
+ p_outfile, " %14.7lE%15.7lE%15.7lE\n",
+ vec_dir_scsrt1[dir_index],
+ vec_dir_absrt1[dir_index],
+ vec_dir_exsrt1[dir_index]
+ );
+ fprintf(
+ p_outfile, " FSAS(1,1)=%15.7lE%15.7lE FSAS(2,1)=%15.7lE%15.7lE\n",
+ real(vec_dir_fsas11[dir_index]), imag(vec_dir_fsas11[dir_index]),
+ real(vec_dir_fsas21[dir_index]), imag(vec_dir_fsas21[dir_index])
+ );
+ fprintf(
+ p_outfile, " SAS(1,1)=%15.7lE%15.7lE SAS(2,1)=%15.7lE%15.7lE\n",
+ real(vec_dir_sas11[dir_index]), imag(vec_dir_sas11[dir_index]),
+ real(vec_dir_sas21[dir_index]), imag(vec_dir_sas21[dir_index])
+ );
+ fprintf(
+ p_outfile, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ vec_dir_qschu1[dir_index],
+ vec_dir_pschu1[dir_index],
+ vec_dir_s0mag1[dir_index]
+ );
+ if (!goto290) {
+ fprintf(
+ p_outfile, " COSAV=%15.7lE, RAPRS=%15.7lE\n",
+ vec_dir_cosav1[dir_index], vec_dir_rapr1[dir_index]
+ );
+ fprintf(
+ p_outfile, " Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n",
+ vec_dir_fl1[dir_index], vec_dir_fr1[dir_index],
+ vec_dir_fk1[dir_index]
+ );
+ fprintf(
+ p_outfile, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
+ vec_dir_fx1[dir_index], vec_dir_fy1[dir_index],
+ vec_dir_fz1[dir_index]
+ );
+ fprintf(
+ p_outfile, " TQEl=%15.7lE, TQEr=%15.7lE, TQEk=%15.7lE\n",
+ vec_dir_tqel1[dir_index], vec_dir_tqer1[dir_index],
+ vec_dir_tqek1[dir_index]
+ );
+ fprintf(
+ p_outfile, " TQSl=%15.7lE, TQSr=%15.7lE, TQSk=%15.7lE\n",
+ vec_dir_tqsl1[dir_index], vec_dir_tqsr1[dir_index],
+ vec_dir_tqsk1[dir_index]
+ );
+ fprintf(
+ p_outfile, " TQEx=%15.7lE, TQEy=%15.7lE, TQEz=%15.7lE\n",
+ vec_dir_tqex1[dir_index], vec_dir_tqey1[dir_index],
+ vec_dir_tqez1[dir_index]
+ );
+ fprintf(
+ p_outfile, " TQSx=%15.7lE, TQSy=%15.7lE, TQSz=%15.7lE\n",
+ vec_dir_tqsx1[dir_index], vec_dir_tqsy1[dir_index],
+ vec_dir_tqsz1[dir_index]
+ );
+ } // goto290 switch
+ if (inpol == 0) {
+ fprintf(p_outfile, " LIN 1\n");
+ } else {
+ fprintf(p_outfile, " CIRC 1\n");
+ }
+ // label 275
+ fprintf(p_outfile, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+ fprintf(
+ p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_dir_scs2[dir_index],
+ vec_dir_abs2[dir_index],
+ vec_dir_exs2[dir_index],
+ vec_dir_albeds2[dir_index]
+ );
+ fprintf(
+ p_outfile, "INSERTION: SCASEC %5d%14.7lE%14.7lE%14.7lE%14.7lE\n",
+ 1, alamb, vec_dir_scs2[dir_index], vec_dir_abs2[dir_index],
+ vec_dir_exs2[dir_index]
+ );
+ fprintf(p_outfile, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
+ fprintf(
+ p_outfile, " %14.7lE%15.7lE%15.7lE\n",
+ vec_dir_scsrt2[dir_index],
+ vec_dir_absrt2[dir_index],
+ vec_dir_exsrt2[dir_index]
+ );
+ fprintf(
+ p_outfile, " FSAS(2,2)=%15.7lE%15.7lE FSAS(1,2)=%15.7lE%15.7lE\n",
+ real(vec_dir_fsas22[dir_index]), imag(vec_dir_fsas22[dir_index]),
+ real(vec_dir_fsas12[dir_index]), imag(vec_dir_fsas12[dir_index])
+ );
+ fprintf(
+ p_outfile, " SAS(2,2)=%15.7lE%15.7lE SAS(1,2)=%15.7lE%15.7lE\n",
+ real(vec_dir_sas22[dir_index]), imag(vec_dir_sas22[dir_index]),
+ real(vec_dir_sas12[dir_index]), imag(vec_dir_sas12[dir_index])
+ );
+ fprintf(
+ p_outfile, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ vec_dir_qschu2[dir_index],
+ vec_dir_pschu2[dir_index],
+ vec_dir_s0mag2[dir_index]
+ );
+ if (!goto290) {
+ fprintf(
+ p_outfile, " COSAV=%15.7lE, RAPRS=%15.7lE\n",
+ vec_dir_cosav2[dir_index], vec_dir_rapr2[dir_index]
+ );
+ fprintf(
+ p_outfile, " Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n",
+ vec_dir_fl2[dir_index], vec_dir_fr2[dir_index],
+ vec_dir_fk2[dir_index]
+ );
+ fprintf(
+ p_outfile, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
+ vec_dir_fx2[dir_index], vec_dir_fy2[dir_index],
+ vec_dir_fz2[dir_index]
+ );
+ fprintf(
+ p_outfile, " TQEl=%15.7lE, TQEr=%15.7lE, TQEk=%15.7lE\n",
+ vec_dir_tqel2[dir_index], vec_dir_tqer2[dir_index],
+ vec_dir_tqek2[dir_index]
+ );
+ fprintf(
+ p_outfile, " TQSl=%15.7lE, TQSr=%15.7lE, TQSk=%15.7lE\n",
+ vec_dir_tqsl2[dir_index], vec_dir_tqsr2[dir_index],
+ vec_dir_tqsk2[dir_index]
+ );
+ fprintf(
+ p_outfile, " TQEx=%15.7lE, TQEy=%15.7lE, TQEz=%15.7lE\n",
+ vec_dir_tqex2[dir_index], vec_dir_tqey2[dir_index],
+ vec_dir_tqez2[dir_index]
+ );
+ fprintf(
+ p_outfile, " TQSx=%15.7lE, TQSy=%15.7lE, TQSz=%15.7lE\n",
+ vec_dir_tqsx2[dir_index], vec_dir_tqsy2[dir_index],
+ vec_dir_tqsz2[dir_index]
+ );
+ } // goto290 switch
+ fprintf(
+ p_outfile, " (RE(FSAS(1,1))-RE(FSAS(2,2)))/RE(FSAS(1,1))=%15.7lE\n",
+ (real(vec_dir_fsas11[dir_index]) - real(vec_dir_fsas22[dir_index])) / real(vec_dir_fsas11[dir_index])
+ );
+ fprintf(
+ p_outfile, " (IM(FSAS(1,1))-IM(FSAS(2,2)))/IM(FSAS(1,1))=%15.7lE\n",
+ (imag(vec_dir_fsas11[dir_index]) - imag(vec_dir_fsas22[dir_index])) / imag(vec_dir_fsas11[dir_index])
+ );
+ fprintf(p_outfile, " MULL\n");
+ for (int i = 0; i < 4; i++) {
+ int mull_index = 16 * dir_index + 4 * i;
+ fprintf(
+ p_outfile, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_dir_mull[mull_index],
+ vec_dir_mull[mull_index + 1],
+ vec_dir_mull[mull_index + 2],
+ vec_dir_mull[mull_index + 3]
+ );
+ } // i MULL loop
+ fprintf(p_outfile, " MULLLR\n");
+ for (int i = 0; i < 4; i++) {
+ int mull_index = 16 * dir_index + 4 * i;
+ fprintf(
+ p_outfile, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_dir_mulllr[mull_index],
+ vec_dir_mulllr[mull_index + 1],
+ vec_dir_mulllr[mull_index + 2],
+ vec_dir_mulllr[mull_index + 3]
+ );
+ } // i MULLR loop
+ if (iavm != 0) {
+ fprintf(p_outfile, " ENSEMBLE AVERAGE, MODE%2d\n", iavm);
+ if (inpol == 0) fprintf(p_outfile, " LIN\n");
+ else fprintf(p_outfile, " CIRC\n");
+ // label 318
+ fprintf(p_outfile, " MULL\n");
+ for (int i = 0; i < 4; i++) {
+ int mul_dir_index = 16 * dir_index + 4 * i;
+ fprintf(
+ p_outfile, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_dir_mull[mul_dir_index],
+ vec_dir_mull[mul_dir_index + 1],
+ vec_dir_mull[mul_dir_index + 2],
+ vec_dir_mull[mul_dir_index + 3]
+ );
+ } // i MULL loop
+ fprintf(p_outfile, " MULLLR\n");
+ for (int i = 0; i < 4; i++) {
+ int mul_dir_index = 16 * dir_index + 4 * i;
+ fprintf(
+ p_outfile, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_dir_mulllr[mul_dir_index],
+ vec_dir_mulllr[mul_dir_index + 1],
+ vec_dir_mulllr[mul_dir_index + 2],
+ vec_dir_mulllr[mul_dir_index + 3]
+ );
+ } // i MULLR loop
+ } // iavm != 0 block
+ done_dirs++;
+ } // jphs loop
+ } // jths loop
+ } // jph loop
+ } // jth loop
+ } // jxi wavelength loop
+ fclose(p_outfile);
+ } else {
+ result = -1;
+ }
+ return result;
+}
+
+#ifdef MPI_VERSION
+int InclusionOutputInfo::mpireceive(const mixMPI* mpidata, int pid) {
+ int result = 0;
+ int chk_inpol, chk_iavm, chk_isam, chk_num_theta, chk_num_thetas;
+ int chk_num_phi, chk_num_phis, chk_ndirs, chk_idfc, chk_configs;
+ double chk_exri;
+ MPI_Recv(&chk_inpol, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_iavm, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_isam, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_num_theta, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_num_thetas, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_num_phi, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_num_phis, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_ndirs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_exri, 1, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_idfc, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_configs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ result += (chk_inpol == inpol) ? 0 : 1;
+ result += (chk_iavm == iavm) ? 0 : 1;
+ result += (chk_isam == isam) ? 0 : 1;
+ result += (chk_num_theta == _num_theta) ? 0 : 1;
+ result += (chk_num_thetas == _num_thetas) ? 0 : 1;
+ result += (chk_num_phi == _num_phi) ? 0 : 1;
+ result += (chk_num_phis == _num_phis) ? 0 : 1;
+ result += (chk_ndirs == ndirs) ? 0 : 1;
+ result += (chk_exri == exri) ? 0 : 1;
+ result += (chk_idfc == idfc) ? 0 : 1;
+ result += (chk_configs == configurations) ? 0 : 1;
+ if (result == 0) {
+ int xi1, offset, chunk_size;
+ MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
+ MPI_Recv(&xi1, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ // Receive vectors of single values per scale
+ MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ offset = xi1 - _first_xi;
+ MPI_Recv(vec_jxi + offset, chunk_size, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_ier + offset, chunk_size, MPI_SHORT, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_vk + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_xi + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_scs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_scs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_abs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_abs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_exs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_exs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_albeds1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_albeds2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_scsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_scsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_absrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_absrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_exsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_exsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_qschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_qschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_pschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_pschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_s0mag1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_s0mag2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_cosav1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_cosav2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_raprs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_raprs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_fk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_fk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_fsas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_fsas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_fsas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_fsas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+ // Receive vectors whose sizes depend on configurations
+ MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ offset = (xi1 - _first_xi) * (configurations + 1);
+ MPI_Recv(vec_sphere_sizes + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_sphere_ref_indices + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+ // Receive vectors whose sizes depend on directions and wavelengths
+ MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ offset = (xi1 - _first_xi) * ndirs;
+ MPI_Recv(vec_dir_scs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_scs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_abs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_abs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_exs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_exs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_albeds1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_albeds2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_scsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_scsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_absrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_absrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_exsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_exsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fsas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fsas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fsas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fsas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_sas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_sas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_sas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_sas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_qschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_qschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_pschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_pschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_s0mag1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_s0mag2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_cosav1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_cosav2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_rapr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_rapr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+ MPI_Recv(vec_dir_fl1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fl2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fx1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fx2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fy1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fy2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fz1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fz2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqel1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqel2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqer1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqer2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqek1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqek2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqex1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqex2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqey1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqey2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqez1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqez2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsl1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsl2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsx1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsx2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsy1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsy2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsz1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_tqsz2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_mull + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_mulllr + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ } else {
+ MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
+ }
+ return result;
+}
+
+int InclusionOutputInfo::mpisend(const mixMPI *mpidata) {
+ int result = 0;
+ int chunk_size;
+ // Send output metadata for configuration cross-check
+ MPI_Send(&inpol, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&iavm, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&isam, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&_num_theta, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&_num_thetas, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&_num_phi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&_num_phis, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&ndirs, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&exri, 1, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&idfc, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&configurations, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ // Wait for process 0 to cross-check the configuration
+ MPI_Recv(&result, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ if (result == 0) {
+ // Process 0 confirmed the consistency of configuration. Send the data.
+ // Send vectors of single values per scale
+ MPI_Send(&_first_xi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&xi_block_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_jxi, xi_block_size, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_ier, xi_block_size, MPI_SHORT, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_vk, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_xi, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_scs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_scs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_abs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_abs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_exs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_exs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_albeds1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_albeds2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_scsrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_scsrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_absrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_absrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_exsrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_exsrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_qschu1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_qschu2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_pschu1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_pschu2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_s0mag1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_s0mag2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_cosav1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_cosav2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_raprs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_raprs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_fk1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_fk2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_fsas11, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_fsas21, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_fsas22, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_fsas12, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+
+ // Send vectors whose sizes depend on configurations
+ chunk_size = xi_block_size * (configurations + 1);
+ MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_sphere_sizes, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_sphere_ref_indices, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+
+ // Send vectors whose sizes depend on directions and wavelengths
+ chunk_size = xi_block_size * ndirs;
+ MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_scs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_scs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_abs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_abs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_exs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_exs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_albeds1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_albeds2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_scsrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_scsrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_absrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_absrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_exsrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_exsrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fsas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fsas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fsas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fsas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_sas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_sas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_sas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_sas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_qschu1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_qschu2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_pschu1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_pschu2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_s0mag1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_s0mag2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_cosav1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_cosav2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_rapr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_rapr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fl1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fl2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fx1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fx2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fy1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fy2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fz1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fz2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqel1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqel2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqer1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqer2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqek1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqek2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqex1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqex2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqey1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqey2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqez1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqez2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsl1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsl2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsx1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsx2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsy1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsy2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsz1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_tqsz2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_mull, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_mulllr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ }
+ return result;
+}
+#endif // MPI_VERSION
+// >>> END OF InclusionOutputInfo CLASS IMPLEMENTATION <<<
diff --git a/src/scripts/pycompare.py b/src/scripts/pycompare.py
index fe68f7338bd6721829dcbe4a1856b015fbef43f3..d8c7072f233e908517c2222c21a983ee8b460f63 100755
--- a/src/scripts/pycompare.py
+++ b/src/scripts/pycompare.py
@@ -141,7 +141,7 @@ def compare_files(config):
print("ERROR: C++ file is shorter than FORTRAN file.")
fortran_file.close()
c_file.close()
- mismatch_count['errors'] = line_count
+ mismatch_count['errors'] = line_count if line_count > 0 else 1
return mismatch_count
f_lines[0] = fortran_file.readline()
c_lines[0] = c_file.readline()
diff --git a/src/testing/test_outputs.cpp b/src/testing/test_cluster_outputs.cpp
similarity index 100%
rename from src/testing/test_outputs.cpp
rename to src/testing/test_cluster_outputs.cpp
diff --git a/src/testing/test_inclusion_outputs.cpp b/src/testing/test_inclusion_outputs.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ac93471acc09b4b03d09e294b08a69b6cab02884
--- /dev/null
+++ b/src/testing/test_inclusion_outputs.cpp
@@ -0,0 +1,65 @@
+#include <string>
+
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
+#ifndef INCLUDE_CONFIGURATION_H_
+#include "../include/Configuration.h"
+#endif
+
+#ifndef INCLUDE_COMMONS_H_
+#include "../include/Commons.h"
+#endif
+
+#ifndef INCLUDE_OUTPUTS_H_
+#include "../include/outputs.h"
+#endif
+
+using namespace std;
+
+int test_inclusion_hdf5_output();
+int test_inclusion_devel();
+
+int main() {
+ int result = 0;
+ result += test_inclusion_hdf5_output(); // 1 if failed
+ result += test_inclusion_devel(); // 10 if failed
+ return result;
+}
+
+int test_inclusion_hdf5_output() {
+ int result = 0;
+ try {
+ const string hdf5_file = "../../test_data/inclusion/c_OINCLU.hd5";
+ InclusionOutputInfo *oi = new InclusionOutputInfo(hdf5_file);
+ oi->write("c_OINCLU", "LEGACY");
+ delete oi;
+ } catch (const exception& ex) {
+ result = 1;
+ }
+ return result;
+}
+
+int test_inclusion_devel() {
+ int result = 0;
+ try {
+ const string geom_data_file = "../../test_data/inclusion/DINCLU";
+ const string scat_data_file = "../../test_data/inclusion/DEDFB";
+ mixMPI *mpidata = new mixMPI();
+ GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(geom_data_file);
+ ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(scat_data_file);
+ InclusionOutputInfo *oi = new InclusionOutputInfo(sconf, gconf, mpidata);
+ delete gconf;
+ delete sconf;
+ delete oi;
+ delete mpidata;
+ } catch (const exception& ex) {
+ result = 10;
+ }
+ return result;
+}
diff --git a/test_data/inclusion/c_OINCLU.hd5 b/test_data/inclusion/c_OINCLU.hd5
new file mode 100644
index 0000000000000000000000000000000000000000..c32b085708391b2755357c05f1fb32a63f5772a0
Binary files /dev/null and b/test_data/inclusion/c_OINCLU.hd5 differ