diff --git a/src/include/outputs.h b/src/include/outputs.h
index 39887cd096a9f1e79d9da2768228aaecb50dafb5..abea6a6f891c1e318c7c7c244f196cc469a38b48 100644
--- a/src/include/outputs.h
+++ b/src/include/outputs.h
@@ -467,9 +467,9 @@ public:
/*! \brief `ClusterOutputInfo` constructor from HDF5 input.
*
- * \param hdf5_file_name: `const string &` Path to the HDF5 file to be read.
+ * \param hdf5_name: `const string &` Path to the HDF5 file to be read.
*/
- ClusterOutputInfo(const std::string &hdf5_file_name);
+ ClusterOutputInfo(const std::string &hdf5_name);
/*! \brief `ClusterOutputInfo` instance destroyer.
*/
@@ -479,13 +479,13 @@ public:
*
* \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 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 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
+ int first_xi = 1, int xi_length = 0
);
/*! \brief Get the size of a `ClusterOutputInfo` instance in bytes.
@@ -647,7 +647,7 @@ public:
int jwtm;
//! \brief Vector of scale (wavelength) indices.
int *vec_jxi;
- //! \brief Vector of error severities (0 - success, 1 - HJV, 2 - DME).
+ //! \brief Vector of error severities (0 - success, 1 - INDME, 2 - OSPV).
short *vec_ier;
//! \brief Vector of vacuum wave numbers.
double *vec_vk;
@@ -897,9 +897,9 @@ public:
/*! \brief `InclusionOutputInfo` constructor from HDF5 input.
*
- * \param hdf5_file_name: `const string &` Path to the HDF5 file to be read.
+ * \param hdf5_name: `const string &` Path to the HDF5 file to be read.
*/
- InclusionOutputInfo(const std::string &hdf5_file_name);
+ InclusionOutputInfo(const std::string &hdf5_name);
/*! \brief `InclusionOutputInfo` instance destroyer.
*/
@@ -909,13 +909,13 @@ public:
*
* \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 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 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
+ int first_xi = 1, int xi_length = 0
);
/*! \brief Get the size of a `ClusterOutputInfo` instance in bytes.
@@ -1004,18 +1004,16 @@ protected:
public:
//! \brief Read-only view on the ID of the first scale
const int &first_xi = _first_xi;
+ //! \brief Number of spheres.
+ int nsph;
//! \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.
@@ -1054,18 +1052,106 @@ public:
int xi_block_size;
//! \brief Index of the wavelength for T-matrix output.
int jwtm;
+ //! \brief Number of sphere types.
+ int configurations;
+ //! \brief Highest expansion order achieved in calculations.
+ int lcalc;
+ //! \brief Harmonic functions argument.
+ dcomplex arg;
//! \brief Vector of scale (wavelength) indices.
int *vec_jxi;
- //! \brief Vector of error severities (0 - success, 1 - HJV, 2 - DME).
+ //! \brief Vector of error severities (0 - success, 1 - DME).
short *vec_ier;
//! \brief Vector of vacuum wave numbers.
double *vec_vk;
//! \brief Vector of computed scales.
double *vec_xi;
- //! \brief Vector of sphere sizes (one for every scale).
+ //! \brief Vector of sphere sizes (one for every configuration and scale).
double *vec_sphere_sizes;
- //! \brief Vector of sphere refractive indices (one for every scale).
+ //! \brief Vector of sphere refractive indices (one for every configuration and scale).
dcomplex *vec_sphere_ref_indices;
+ //! \brief Vector of sphere scattering cross-sections.
+ double *vec_scs;
+ //! \brief Vector of sphere absorption cross-sections.
+ double *vec_abs;
+ //! \brief Vector of sphere extinction cross-sections.
+ double *vec_exs;
+ //! \brief Vector of sphere albedos.
+ double *vec_albeds;
+ //! \brief Vector of sphere scattering-to-geometric cross-sections.
+ double *vec_scsrt;
+ //! \brief Vector of sphere absorption-to-geometric cross-sections.
+ double *vec_absrt;
+ //! \brief Vector of sphere extinction-to-geometric cross-sections.
+ double *vec_exsrt;
+ //! \brief Vector of sphere forward scattering amplitudes.
+ dcomplex *vec_fsas;
+ //! \brief Vector of sphere QSCHU.
+ double *vec_qschu;
+ //! \brief Vector of sphere PSCHU.
+ double *vec_pschu;
+ //! \brief Vector of sphere S0MAG.
+ double *vec_s0mag;
+ //! \brief Vector of sphere average asymmetry parameter.
+ double *vec_cosav;
+ //! \brief Vector of sphere average radiation pressure force (N).
+ double *vec_raprs;
+ //! \brief Vector of sphere average extinction torque along incidence direction (parallel polarization).
+ double *vec_tqek1;
+ //! \brief Vector of sphere average extinction torque along incidence direction (perpendicular polarization).
+ double *vec_tqek2;
+ //! \brief Vector of sphere average scattering torque along incidence direction (parallel polarization).
+ double *vec_tqsk1;
+ //! \brief Vector of sphere average scattering torque along incidence direction (perpendicular polarization).
+ double *vec_tqsk2;
+ //| \brief Vector of total forward scattering amplitudes.
+ dcomplex *vec_fsat;
+ //! \brief Vector of total QSCHU.
+ double *vec_qschut;
+ //! \brief Vector of total PSCHU.
+ double *vec_pschut;
+ //! \brief Vector of total S0MAG.
+ double *vec_s0magt;
+ //! \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 sphere differential scattering amplitude with polarization parallel to parallel incidence field.
+ dcomplex *vec_dir_sas11;
+ //! \brief Vector of sphere differential scattering amplitude with polarization perpendicular to the parallel incidence field.
+ dcomplex *vec_dir_sas21;
+ //! \brief Vector of sphere differential scattering amplitude with polarization perpendicular to perpendicular incidence field.
+ dcomplex *vec_dir_sas12;
+ //! \brief Vector of sphere differential scattering amplitude with polarization parallel the perpendicular incidence field.
+ dcomplex *vec_dir_sas22;
+ //! \brief Vector of differential radiation pressure force components along the X axis.
+ double *vec_dir_fx;
+ //! \brief Vector of differential radiation pressure force components along the Y axis.
+ double *vec_dir_fy;
+ //! \brief Vector of differential radiation pressure force components along the Z axis.
+ double *vec_dir_fz;
+ //! \brief Vector of sphere Mueller transormation matrices referred to meridional plane.
+ double *vec_dir_muls;
+ //! \brief Vector of sphere Mueller transormation matrices referred to scattering plane.
+ double *vec_dir_mulslr;
/*! \brief `SphereOutputInfo` default instance constructor.
*
@@ -1082,9 +1168,9 @@ public:
/*! \brief `SphereOutputInfo` constructor from HDF5 input.
*
- * \param hdf5_file_name: `const string &` Path to the HDF5 file to be read.
+ * \param hdf5_name: `const string &` Path to the HDF5 file to be read.
*/
- SphereOutputInfo(const std::string &hdf5_file_name);
+ SphereOutputInfo(const std::string &hdf5_name);
/*! \brief `InclusionOutputInfo` instance destroyer.
*/
@@ -1094,13 +1180,13 @@ public:
*
* \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 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 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
+ int first_xi = 1, int xi_length = 0
);
/*! \brief Get the size of a `ClusterOutputInfo` instance in bytes.
@@ -1136,7 +1222,7 @@ public:
* \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);
+ int mpireceive(const mixMPI *mpidata, int pid);
/*! \brief Send output data to process 0.
*
diff --git a/src/libnptm/outputs.cpp b/src/libnptm/outputs.cpp
index 4d6e61a2fb6273f9e665e038098d6f76ad0a868b..58f40cd5865b1455cd0a0f0b47c2ec099511c241 100644
--- a/src/libnptm/outputs.cpp
+++ b/src/libnptm/outputs.cpp
@@ -288,9 +288,9 @@ ClusterOutputInfo::ClusterOutputInfo(
vec_dir_mulclr = new double[16 * ndirs * xi_block_size]();
}
-ClusterOutputInfo::ClusterOutputInfo(const std::string &hdf5_file_name) {
+ClusterOutputInfo::ClusterOutputInfo(const std::string &hdf5_name) {
unsigned int flags = H5F_ACC_RDONLY;
- HDFFile *hdf_file = new HDFFile(hdf5_file_name, flags);
+ HDFFile *hdf_file = new HDFFile(hdf5_name, flags);
herr_t status = hdf_file->get_status();
string str_name, str_type;
if (status == 0) {
@@ -696,7 +696,7 @@ ClusterOutputInfo::ClusterOutputInfo(const std::string &hdf5_file_name) {
delete hdf_file;
} else {
if (hdf_file != NULL) delete hdf_file;
- UnrecognizedFormatException ex("Error: " + hdf5_file_name + " not recognized as a valid HDF5 file!");
+ UnrecognizedFormatException ex("Error: " + hdf5_name + " not recognized as a valid HDF5 file!");
throw ex;
}
}
@@ -961,6 +961,7 @@ long ClusterOutputInfo::compute_size() {
int ClusterOutputInfo::insert(const ClusterOutputInfo &rhs) {
int result = 0;
+ result += (rhs.nsph == nsph) ? 0 : 1;
result += (rhs.inpol == inpol) ? 0 : 1;
result += (rhs.iavm == iavm) ? 0 : 1;
result += (rhs.isam == isam) ? 0 : 1;
@@ -2355,9 +2356,10 @@ int ClusterOutputInfo::write_legacy(const std::string &output) {
#ifdef MPI_VERSION
int ClusterOutputInfo::mpireceive(const mixMPI *mpidata, int pid) {
int result = 0;
- int chk_inpol, chk_iavm, chk_isam, chk_num_theta, chk_num_thetas;
+ int chk_nsph, 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_nsph, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
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);
@@ -2369,6 +2371,7 @@ int ClusterOutputInfo::mpireceive(const mixMPI *mpidata, int pid) {
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_nsph == nsph) ? 0 : 1;
result += (chk_inpol == inpol) ? 0 : 1;
result += (chk_iavm == iavm) ? 0 : 1;
result += (chk_isam == isam) ? 0 : 1;
@@ -2559,6 +2562,7 @@ int ClusterOutputInfo::mpisend(const mixMPI *mpidata) {
int result = 0;
int chunk_size;
// Send output metadata for configuration cross-check
+ MPI_Send(&nsph, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
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);
@@ -2932,9 +2936,9 @@ InclusionOutputInfo::InclusionOutputInfo(
vec_dir_mulllr = new double[16 * ndirs * xi_block_size]();
}
-InclusionOutputInfo::InclusionOutputInfo(const std::string &hdf5_file_name) {
+InclusionOutputInfo::InclusionOutputInfo(const std::string &hdf5_name) {
unsigned int flags = H5F_ACC_RDONLY;
- HDFFile *hdf_file = new HDFFile(hdf5_file_name, flags);
+ HDFFile *hdf_file = new HDFFile(hdf5_name, flags);
herr_t status = hdf_file->get_status();
string str_name, str_type;
if (status == 0) {
@@ -3060,18 +3064,24 @@ InclusionOutputInfo::InclusionOutputInfo(const std::string &hdf5_file_name) {
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);
+ // 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;
+ }
str_type = "FLOAT64_(" + to_string(ndirs) + ")";
vec_dir_scand = new double[ndirs];
status = hdf_file->read("VEC_DIR_SCAND", str_type, vec_dir_scand);
@@ -3236,7 +3246,7 @@ InclusionOutputInfo::InclusionOutputInfo(const std::string &hdf5_file_name) {
delete hdf_file;
} else {
if (hdf_file != NULL) delete hdf_file;
- UnrecognizedFormatException ex("Error: " + hdf5_file_name + " not recognized as a valid HDF5 file!");
+ UnrecognizedFormatException ex("Error: " + hdf5_name + " not recognized as a valid HDF5 file!");
throw ex;
}
}
@@ -3420,6 +3430,7 @@ long InclusionOutputInfo::compute_size(
int InclusionOutputInfo::insert(const InclusionOutputInfo &rhs) {
int result = 0;
+ result += (rhs.nsph == nsph) ? 0 : 1;
result += (rhs.inpol == inpol) ? 0 : 1;
result += (rhs.iavm == iavm) ? 0 : 1;
result += (rhs.isam == isam) ? 0 : 1;
@@ -3776,22 +3787,6 @@ int InclusionOutputInfo::write_hdf5(const std::string &file_name) {
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);
@@ -4481,9 +4476,10 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
#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_nsph, 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_nsph, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
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);
@@ -4495,6 +4491,7 @@ int InclusionOutputInfo::mpireceive(const mixMPI* mpidata, int pid) {
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_nsph == nsph) ? 0 : 1;
result += (chk_inpol == inpol) ? 0 : 1;
result += (chk_iavm == iavm) ? 0 : 1;
result += (chk_isam == isam) ? 0 : 1;
@@ -4638,6 +4635,7 @@ int InclusionOutputInfo::mpisend(const mixMPI *mpidata) {
int result = 0;
int chunk_size;
// Send output metadata for configuration cross-check
+ MPI_Send(&nsph, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
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);
@@ -4775,3 +4773,1130 @@ int InclusionOutputInfo::mpisend(const mixMPI *mpidata) {
}
#endif // MPI_VERSION
// >>> END OF InclusionOutputInfo CLASS IMPLEMENTATION <<<
+
+// >>> SphereOutputInfo CLASS IMPLEMENTATION <<<
+SphereOutputInfo::SphereOutputInfo(
+ ScattererConfiguration *sc, GeometryConfiguration *gc,
+ const mixMPI *mpidata, int first_xi, int xi_length
+) {
+ _first_xi = first_xi;
+ nsph = gc->number_of_spheres;
+ lm = gc->l_max;
+ inpol = gc->in_pol;
+ npnt = gc->npnt;
+ npntts = gc->npntts;
+ isam = gc->isam;
+ idfc = sc->idfc;
+ th = gc->in_theta_start;
+ thstp = gc->in_theta_step;
+ thlst = gc->in_theta_end;
+ _num_theta = (thstp == 0.0) ? 1 : 1 + (int)((thlst - th) / thstp);
+ ths = gc->sc_theta_start;
+ thsstp = gc->sc_theta_step;
+ thslst = gc->sc_theta_end;
+ _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;
+ _num_phi = (phstp == 0.0) ? 1 : 1 + (int)((phlst - ph) / phstp);
+ phs = gc->sc_phi_start;
+ phsstp = gc->sc_phi_step;
+ phslst = gc->sc_phi_end;
+ _num_phis = (phsstp == 0.0) ? 1 : 1 + (int)((phslst - phs) / phsstp);
+ ndirs = _num_theta * _num_thetas * _num_phi * _num_phis;
+ configurations = sc->configurations;
+ double exdc = sc->exdc;
+ exri = sqrt(exdc);
+ nxi = sc->number_of_scales;
+ xi_block_size = (xi_length == 0) ? nxi : xi_length;
+ jwtm = gc->jwtm;
+ lcalc = 0;
+ arg = 0.0 + I * 0.0;
+ 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]();
+ vec_sphere_sizes = new double[configurations * xi_block_size]();
+ vec_sphere_ref_indices = new dcomplex[configurations * xi_block_size]();
+ vec_scs = new double[configurations * xi_block_size]();
+ vec_abs = new double[configurations * xi_block_size]();
+ vec_exs = new double[configurations * xi_block_size]();
+ vec_albeds = new double[configurations * xi_block_size]();
+ vec_scsrt = new double[configurations * xi_block_size]();
+ vec_absrt = new double[configurations * xi_block_size]();
+ vec_exsrt = new double[configurations * xi_block_size]();
+ vec_fsas = new dcomplex[configurations * xi_block_size]();
+ vec_qschu = new double[configurations * xi_block_size]();
+ vec_pschu = new double[configurations * xi_block_size]();
+ vec_s0mag = new double[configurations * xi_block_size]();
+ vec_cosav = new double[configurations * xi_block_size]();
+ vec_raprs = new double[configurations * xi_block_size]();
+ vec_tqek1 = new double[configurations * xi_block_size]();
+ vec_tqek2 = new double[configurations * xi_block_size]();
+ vec_tqsk1 = new double[configurations * xi_block_size]();
+ vec_tqsk2 = new double[configurations * xi_block_size]();
+ if (nsph == 1) {
+ vec_fsat = NULL;
+ vec_qschut = NULL;
+ vec_pschut = NULL;
+ vec_s0magt = NULL;
+ } else {
+ vec_fsat = new dcomplex[xi_block_size]();
+ vec_qschut = new double[xi_block_size]();
+ vec_pschut = new double[xi_block_size]();
+ vec_s0magt = new double[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]();
+ vec_dir_scand = new double[ndirs]();
+ vec_dir_cfmp = new double[ndirs]();
+ vec_dir_cfsp = new double[ndirs]();
+ vec_dir_sfmp = 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_sas11 = new dcomplex[nsph * ndirs * xi_block_size]();
+ vec_dir_sas21 = new dcomplex[nsph * ndirs * xi_block_size]();
+ vec_dir_sas12 = new dcomplex[nsph * ndirs * xi_block_size]();
+ vec_dir_sas22 = new dcomplex[nsph * ndirs * xi_block_size]();
+ vec_dir_fx = new double[nsph * _num_theta * _num_phi * xi_block_size]();
+ vec_dir_fy = new double[nsph * _num_theta * _num_phi * xi_block_size]();
+ vec_dir_fz = new double[nsph * _num_theta * _num_phi * xi_block_size]();
+ vec_dir_muls = new double[16 * nsph * ndirs * xi_block_size]();
+ vec_dir_mulslr = new double[16 * nsph * ndirs * xi_block_size]();
+ // 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;
+ }
+}
+
+SphereOutputInfo::SphereOutputInfo(const std::string &hdf5_name) {
+ // DA QUI ; TODO: AGGIUSTARE TUTTE LE DIMENSIONI LEGATE A configurations
+ unsigned int flags = H5F_ACC_RDONLY;
+ HDFFile *hdf_file = new HDFFile(hdf5_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("LM", "INT32_(1)", &lm);
+ 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("ISAM", "INT32_(1)", &isam);
+ status = hdf_file->read("JWTM", "INT32_(1)", &jwtm);
+ 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("NUM_CONF", "FLOAT64_(1)", &configurations);
+ 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;
+ lcalc = 0;
+ arg = 0.0 + I * 0.0;
+ 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);
+ str_type = "FLOAT64_(" + to_string(configurations * xi_block_size) + ")";
+ vec_sphere_sizes = new double[xi_block_size];
+ status = hdf_file->read("VEC_SPH_SIZES", str_type, vec_sphere_sizes);
+ str_type = "FLOAT64_(" + to_string(2 * configurations * xi_block_size) + ")";
+ vec_sphere_ref_indices = new dcomplex[xi_block_size];
+ status = hdf_file->read("VEC_SPH_REFRI", str_type, vec_sphere_ref_indices);
+ str_type = "FLOAT64_(" + to_string(configurations * xi_block_size) + ")";
+ vec_scs = new double[xi_block_size];
+ status = hdf_file->read("VEC_SCS", str_type, vec_scs);
+ vec_abs = new double[xi_block_size];
+ status = hdf_file->read("VEC_ABS", str_type, vec_abs);
+ vec_exs = new double[xi_block_size];
+ status = hdf_file->read("VEC_EXS", str_type, vec_exs);
+ vec_albeds = new double[xi_block_size];
+ status = hdf_file->read("VEC_ALBEDS", str_type, vec_albeds);
+ vec_scsrt = new double[xi_block_size];
+ status = hdf_file->read("VEC_SCSRT", str_type, vec_scsrt);
+ vec_absrt = new double[xi_block_size];
+ status = hdf_file->read("VEC_ABSRT", str_type, vec_absrt);
+ vec_exsrt = new double[xi_block_size];
+ status = hdf_file->read("VEC_EXSRT", str_type, vec_exsrt);
+ str_type = "FLOAT64_(" + to_string(2 * configurations * xi_block_size) + ")";
+ vec_fsas = new dcomplex[xi_block_size];
+ status = hdf_file->read("VEC_FSAS", str_type, vec_fsas);
+ str_type = "FLOAT64_(" + to_string(configurations * xi_block_size) + ")";
+ vec_qschu = new double[xi_block_size];
+ status = hdf_file->read("VEC_QSCHU", str_type, vec_qschu);
+ vec_pschu = new double[xi_block_size];
+ status = hdf_file->read("VEC_PSCHU", str_type, vec_pschu);
+ vec_s0mag = new double[xi_block_size];
+ status = hdf_file->read("VEC_S0MAG", str_type, vec_s0mag);
+ vec_cosav = new double[xi_block_size];
+ status = hdf_file->read("VEC_COSAV", str_type, vec_cosav);
+ vec_raprs = new double[xi_block_size];
+ status = hdf_file->read("VEC_RAPRS", str_type, vec_raprs);
+ vec_tqek1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_TQEK1", str_type, vec_tqek1);
+ vec_tqek2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_TQEK2", str_type, vec_tqek2);
+ vec_tqsk1 = new double[xi_block_size];
+ status = hdf_file->read("VEC_TQSK1", str_type, vec_tqsk1);
+ vec_tqsk2 = new double[xi_block_size];
+ status = hdf_file->read("VEC_TQSK2", str_type, vec_tqsk2);
+ if (nsph != 1) {
+ str_type = "FLOAT64_(" + to_string(2 * xi_block_size) + ")";
+ vec_fsat = new dcomplex[xi_block_size];
+ status = hdf_file->read("VEC_FSAT", str_type, vec_fsat);
+ str_type = "FLOAT64_(" + to_string(xi_block_size) + ")";
+ vec_qschut = new double[xi_block_size];
+ status = hdf_file->read("VEC_QSCHUT", str_type, vec_qschut);
+ vec_pschut = new double[xi_block_size];
+ status = hdf_file->read("VEC_PSCHUT", str_type, vec_pschut);
+ vec_s0magt = new double[xi_block_size];
+ status = hdf_file->read("VEC_S0MAGT", str_type, vec_s0magt);
+ } else {
+ vec_fsat = NULL;
+ vec_qschut = NULL;
+ vec_pschut = NULL;
+ vec_s0magt = NULL;
+ }
+ // Initialize directions (they are scale-independent)
+ _num_theta = (thstp == 0.0) ? 1 : 1 + (int)((thlst - th) / thstp);
+ _num_thetas = (thsstp == 0.0) ? 1 : 1 + (int)((thslst - ths) / thsstp);
+ _num_phi = (phstp == 0.0) ? 1 : 1 + (int)((phlst - ph) / phstp);
+ _num_phis = (phsstp == 0.0) ? 1 : 1 + (int)((phslst - phs) / phsstp);
+ vec_dir_tidg = new double[_num_theta];
+ vec_dir_tsdg = new double[_num_thetas];
+ vec_dir_pidg = new double[_num_phi];
+ vec_dir_psdg = new double[_num_phis];
+ 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;
+ }
+ 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(2 * nsph * ndirs * xi_block_size) + ")";
+ vec_dir_sas11 = new dcomplex[nsph * ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SAS11", str_type, vec_dir_sas11);
+ vec_dir_sas21 = new dcomplex[nsph * ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SAS21", str_type, vec_dir_sas21);
+ vec_dir_sas12 = new dcomplex[nsph * ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SAS12", str_type, vec_dir_sas12);
+ vec_dir_sas22 = new dcomplex[nsph * ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_SAS22", str_type, vec_dir_sas22);
+ str_type = "FLOAT64_(" + to_string(nsph * _num_theta * _num_phi * xi_block_size) + ")";
+ vec_dir_fx = new double[nsph * _num_theta * _num_phi * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FX", str_type, vec_dir_fx);
+ vec_dir_fy = new double[nsph * _num_theta * _num_phi * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FY", str_type, vec_dir_fy);
+ vec_dir_fz = new double[nsph * _num_theta * _num_phi * xi_block_size];
+ status = hdf_file->read("VEC_DIR_FZ", str_type, vec_dir_fz);
+ str_type = "FLOAT64_(" + to_string(16 * nsph * ndirs * xi_block_size) + ")";
+ vec_dir_muls = new double[16 * nsph * ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_MULS", str_type, vec_dir_muls);
+ vec_dir_mulslr = new double[16 * nsph * ndirs * xi_block_size];
+ status = hdf_file->read("VEC_DIR_MULSLR", str_type, vec_dir_mulslr);
+ status = hdf_file->close();
+ delete hdf_file;
+ } else {
+ if (hdf_file != NULL) delete hdf_file;
+ UnrecognizedFormatException ex("Error: " + hdf5_name + " not recognized as a valid HDF5 file!");
+ throw ex;
+ }
+}
+
+SphereOutputInfo::~SphereOutputInfo() {
+ delete[] vec_jxi;
+ delete[] vec_ier;
+ delete[] vec_vk;
+ delete[] vec_xi;
+ delete[] vec_sphere_sizes;
+ delete[] vec_sphere_ref_indices;
+ delete[] vec_scs;
+ delete[] vec_abs;
+ delete[] vec_exs;
+ delete[] vec_albeds;
+ delete[] vec_scsrt;
+ delete[] vec_absrt;
+ delete[] vec_exsrt;
+ delete[] vec_fsas;
+ delete[] vec_qschu;
+ delete[] vec_pschu;
+ delete[] vec_s0mag;
+ delete[] vec_cosav;
+ delete[] vec_raprs;
+ delete[] vec_tqek1;
+ delete[] vec_tqek2;
+ delete[] vec_tqsk1;
+ delete[] vec_tqsk2;
+ if (nsph != 1) {
+ delete[] vec_fsat;
+ delete[] vec_qschut;
+ delete[] vec_pschut;
+ delete[] vec_s0magt;
+ }
+ 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_cfsp;
+ delete[] vec_dir_sfmp;
+ delete[] vec_dir_sfsp;
+ delete[] vec_dir_un;
+ delete[] vec_dir_uns;
+ delete[] vec_dir_sas11;
+ delete[] vec_dir_sas21;
+ delete[] vec_dir_sas12;
+ delete[] vec_dir_sas22;
+ delete[] vec_dir_fx;
+ delete[] vec_dir_fy;
+ delete[] vec_dir_fz;
+ delete[] vec_dir_muls;
+ delete[] vec_dir_mulslr;
+}
+
+long SphereOutputInfo::compute_size(
+ ScattererConfiguration *sc, GeometryConfiguration *gc,
+ int first_xi, int xi_length
+) {
+ // Size of the configuration set
+ long result = 18 * sizeof(int);
+ result += 12 * sizeof(double);
+ result += 47 * sizeof(long);
+ result += sizeof(dcomplex);
+ // Get configuration parameters
+ int _nsph = gc->number_of_spheres;
+ double _th = gc->in_theta_start;
+ double _thstp = gc->in_theta_step;
+ double _thlst = gc->in_theta_end;
+ int num_theta = (_thstp == 0.0) ? 1 : 1 + (int)((_thlst - _th) / _thstp);
+ double _ths = gc->sc_theta_start;
+ double _thsstp = gc->sc_theta_step;
+ double _thslst = gc->sc_theta_end;
+ int num_thetas = (_thsstp == 0.0) ? 1 : 1 + (int)((_thslst - _ths) / _thsstp);
+ double _ph = gc->in_phi_start;
+ double _phstp = gc->in_phi_step;
+ double _phlst = gc->in_phi_end;
+ int num_phi = (_phstp == 0.0) ? 1 : 1 + (int)((_phlst - _ph) / _phstp);
+ double _phs = gc->sc_phi_start;
+ double _phsstp = gc->sc_phi_step;
+ double _phslst = gc->sc_phi_end;
+ int num_phis = (_phsstp == 0.0) ? 1 : 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 _nconf = sc->configurations;
+ // Size of the data set
+ result += _xi_block_size * (sizeof(short) + sizeof(int));
+ result += 2 * _xi_block_size * sizeof(double);
+ result += 16 * _nconf * _xi_block_size * sizeof(double);
+ result += 2 * _nconf * _xi_block_size * sizeof(dcomplex);
+ result += (num_theta + num_thetas + num_phi + num_phis) * sizeof(double);
+ result += 11 * _ndirs * sizeof(double);
+ result += 4 * _nsph * _ndirs * _xi_block_size * sizeof(dcomplex);
+ result += 3 * _nsph * num_theta * num_phi * sizeof(double);
+ result += 32 * _nsph * _ndirs * _xi_block_size * sizeof(double);
+ if (_nsph != 1) {
+ result += _xi_block_size * sizeof(dcomplex);
+ result += 3 * _xi_block_size * sizeof(double);
+ }
+ return result;
+}
+
+long SphereOutputInfo::compute_size() {
+ // Size of the configuration set
+ long result = 18 * sizeof(int);
+ result += 12 * sizeof(double);
+ result += 47 * sizeof(long);
+ result += sizeof(dcomplex);
+ // Size of the data set
+ result += xi_block_size * (sizeof(short) + sizeof(int));
+ result += 2 * xi_block_size * sizeof(double);
+ result += 16 * configurations * xi_block_size * sizeof(double);
+ result += 2 * configurations * xi_block_size * sizeof(dcomplex);
+ result += (_num_theta + _num_thetas + _num_phi + _num_phis) * sizeof(double);
+ result += 11 * ndirs * sizeof(double);
+ result += 4 * nsph * ndirs * xi_block_size * sizeof(dcomplex);
+ result += 3 * nsph * _num_theta * _num_phi * sizeof(double);
+ result += 32 * nsph * ndirs * xi_block_size * sizeof(double);
+ if (nsph != 1) {
+ result += xi_block_size * sizeof(dcomplex);
+ result += 3 * xi_block_size * sizeof(double);
+ }
+ return result;
+}
+
+int SphereOutputInfo::insert(const SphereOutputInfo &rhs) {
+ int result = 0;
+ result += (rhs.nsph == nsph) ? 0 : 1;
+ result += (rhs.inpol == inpol) ? 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 sizes depend 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));
+ if (nsph != 1) {
+ memcpy(vec_fsat + offset, rhs.vec_fsat, chunk_size * sizeof(dcomplex));
+ memcpy(vec_qschut + offset, rhs.vec_qschut, chunk_size * sizeof(double));
+ memcpy(vec_pschut + offset, rhs.vec_pschut, chunk_size * sizeof(double));
+ memcpy(vec_s0magt + offset, rhs.vec_s0magt, chunk_size * sizeof(double));
+ }
+
+ // Insert vectors whose sizes depend on configurations and wavelengths
+ offset = (xi1 - _first_xi) * configurations;
+ chunk_size = rhs.xi_block_size * configurations;
+ 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));
+ memcpy(vec_scs + offset, rhs.vec_scs, chunk_size * sizeof(double));
+ memcpy(vec_abs + offset, rhs.vec_abs, chunk_size * sizeof(double));
+ memcpy(vec_exs + offset, rhs.vec_exs, chunk_size * sizeof(double));
+ memcpy(vec_albeds + offset, rhs.vec_albeds, chunk_size * sizeof(double));
+ memcpy(vec_scsrt + offset, rhs.vec_scsrt, chunk_size * sizeof(double));
+ memcpy(vec_absrt + offset, rhs.vec_absrt, chunk_size * sizeof(double));
+ memcpy(vec_exsrt + offset, rhs.vec_exsrt, chunk_size * sizeof(double));
+ memcpy(vec_fsas + offset, rhs.vec_fsas, chunk_size * sizeof(dcomplex));
+ memcpy(vec_qschu + offset, rhs.vec_qschu, chunk_size * sizeof(double));
+ memcpy(vec_pschu + offset, rhs.vec_pschu, chunk_size * sizeof(double));
+ memcpy(vec_s0mag + offset, rhs.vec_s0mag, chunk_size * sizeof(double));
+ memcpy(vec_cosav + offset, rhs.vec_cosav, chunk_size * sizeof(double));
+ memcpy(vec_raprs + offset, rhs.vec_raprs, chunk_size * sizeof(double));
+ memcpy(vec_tqek1 + offset, rhs.vec_tqek1, chunk_size * sizeof(double));
+ memcpy(vec_tqek2 + offset, rhs.vec_tqek2, chunk_size * sizeof(double));
+ memcpy(vec_tqsk1 + offset, rhs.vec_tqsk1, chunk_size * sizeof(double));
+ memcpy(vec_tqsk2 + offset, rhs.vec_tqsk2, chunk_size * sizeof(double));
+
+ // Insert vectors whose sizes depend on NSPH, directions and wavelengths
+ offset = (xi1 - _first_xi) * nsph * ndirs;
+ chunk_size = rhs.xi_block_size * nsph * ndirs;
+ 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_muls + 16 * offset, rhs.vec_dir_muls, 16 * chunk_size * sizeof(double));
+ memcpy(vec_dir_mulslr + 16 * offset, rhs.vec_dir_mulslr, 16 * chunk_size * sizeof(double));
+
+ // Insert vectors whose sizes depend on NSPH, incidence directions and wavelengths
+ offset = (xi1 - _first_xi) * nsph * _num_theta * _num_phi;
+ chunk_size = rhs.xi_block_size * nsph * _num_theta * _num_phi;
+ memcpy(vec_dir_fx + offset, rhs.vec_dir_fx, chunk_size * sizeof(double));
+ memcpy(vec_dir_fy + offset, rhs.vec_dir_fy, chunk_size * sizeof(double));
+ memcpy(vec_dir_fz + offset, rhs.vec_dir_fz, chunk_size * sizeof(double));
+ // TODO: fix the vector sizes in HDF5 writer and MPI communicators
+ }
+ return result;
+}
+
+int SphereOutputInfo::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 SphereOutputInfo::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("LM");
+ rec_type_list->append("INT32_(1)");
+ rec_ptr_list->append(&lm);
+ 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("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("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->set(0, "NUM_CONF");
+ rec_type_list->set(0, "INT32_(1)");
+ rec_ptr_list->set(0, &configurations);
+ 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("VEC_SPH_SIZES");
+ rec_type_list->append("FLOAT64_(" + to_string(configurations * xi_block_size) + ")");
+ rec_ptr_list->append(vec_sphere_sizes);
+ rec_name_list->append("VEC_SPH_REFRI");
+ rec_type_list->append("FLOAT64_(" + to_string(2 * configurations * xi_block_size) + ")");
+ rec_ptr_list->append(vec_sphere_ref_indices);
+ str_type = "FLOAT64_(" + to_string(configurations * xi_block_size) + ")";
+ rec_name_list->append("VEC_SCS");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_scs);
+ rec_name_list->append("VEC_ABS");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_abs);
+ rec_name_list->append("VEC_EXS");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_exs);
+ rec_name_list->append("VEC_ALBEDS");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_albeds);
+ rec_name_list->append("VEC_SCSRT");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_scsrt);
+ rec_name_list->append("VEC_ABSRT");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_absrt);
+ rec_name_list->append("VEC_EXSRT");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_exsrt);
+ str_type = "FLOAT64_(" + to_string(2 * configurations * xi_block_size) + ")";
+ rec_name_list->append("VEC_FSAS");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_fsas);
+ str_type = "FLOAT64_(" + to_string(configurations * xi_block_size) + ")";
+ rec_name_list->append("VEC_QSCHU");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_qschu);
+ rec_name_list->append("VEC_PSCHU");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_pschu);
+ rec_name_list->append("VEC_S0MAG");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_s0mag);
+ rec_name_list->append("VEC_COSAV");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_cosav);
+ rec_name_list->append("VEC_RAPRS");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_raprs);
+ rec_name_list->append("VEC_TQEK1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_tqek1);
+ rec_name_list->append("VEC_TQEK2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_tqek2);
+ rec_name_list->append("VEC_TQSK1");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_tqsk1);
+ rec_name_list->append("VEC_TQSK2");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_tqsk2);
+ if (nsph != 1) {
+ str_type = "FLOAT64_(" + to_string(2 * xi_block_size) + ")";
+ rec_name_list->append("VEC_FSAT");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_fsat);
+ str_type = "FLOAT64_(" + to_string(xi_block_size) + ")";
+ rec_name_list->append("VEC_QSCHUT");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_qschut);
+ rec_name_list->append("VEC_PSCHUT");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_pschut);
+ rec_name_list->append("VEC_S0MAGT");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_s0magt);
+ }
+ 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(2 * nsph * ndirs * xi_block_size) + ")";
+ 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(nsph * _num_theta * _num_phi * xi_block_size) + ")";
+ rec_name_list->append("VEC_DIR_FX");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fx);
+ rec_name_list->append("VEC_DIR_FY");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fy);
+ rec_name_list->append("VEC_DIR_FZ");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_fz);
+ str_type = "FLOAT64_(" + to_string(16 * nsph * ndirs * xi_block_size) + ")";
+ rec_name_list->append("VEC_DIR_MULS");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_muls);
+ rec_name_list->append("VEC_DIR_MULSLR");
+ rec_type_list->append(str_type);
+ rec_ptr_list->append(vec_dir_mulslr);
+
+ // 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 SphereOutputInfo::write_legacy(const std::string &file_name) {
+ const dcomplex cc0 = 0.0 + I * 0.0;
+ int result = 0;
+ FILE *p_outfile = fopen(file_name.c_str(), "w");
+ if (p_outfile != NULL) {
+ if (vec_jxi[0] == 1) {
+ fprintf(p_outfile, " READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM\n");
+ fprintf(
+ p_outfile, " %5d%5d%5d%5d%5d%5d\n",
+ nsph, lm, inpol, npnt, npntts, isam
+ );
+ fprintf(p_outfile, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
+ fprintf(
+ p_outfile, " %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",
+ th, thstp, thlst, ths, thsstp, thslst
+ );
+ fprintf(p_outfile, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
+ fprintf(
+ p_outfile, " %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.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 (inpol == 0) fprintf(p_outfile, " LIN\n \n");
+ else fprintf(p_outfile, " CIRC\n \n");
+ 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;
+ fprintf(p_outfile, "========== JXI =%3d ====================\n", jxi + 1);
+ if (idfc >= 0) {
+ fprintf(p_outfile, " VK=%15.7lE, XI=%15.7lE\n", vec_xi[jxi], vec_vk[jxi]);
+ } else { // IDFC < 0
+ fprintf(p_outfile, " XI=%15.7lE\n", vec_xi[jxi]);
+ }
+ if (vec_ier[jxi] == 1) {
+ fprintf(p_outfile, " STOP IN DME\n");
+ fprintf(
+ p_outfile, " AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n",
+ (int)vec_ier[jxi], lcalc, real(arg), imag(arg)
+ );
+ }
+ for (int ci = 0; ci < configurations; ci++) {
+ int cindex = jxi * configurations + ci;
+ fprintf(p_outfile, " SPHERE %2d\n", ci + 1);
+ if (vec_sphere_ref_indices[cindex] == cc0) {
+ fprintf(p_outfile, " SIZE=%15.7lE\n", vec_sphere_sizes[cindex]);
+ } else {
+ fprintf(
+ p_outfile, " 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, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+ fprintf(
+ p_outfile, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_scs[cindex], vec_abs[cindex], vec_exs[cindex], vec_albeds[cindex]
+ );
+ fprintf(p_outfile, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
+ fprintf(
+ p_outfile, " %14.7lE%15.7lE%15.7lE\n",
+ vec_scsrt[cindex], vec_absrt[cindex], vec_exsrt[cindex]
+ );
+ fprintf(
+ p_outfile, " FSAS=%15.7lE%15.7lE\n",
+ real(vec_fsas[cindex]), imag(vec_fsas[cindex])
+ );
+ fprintf(
+ p_outfile, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ vec_qschu[cindex], vec_pschu[cindex], vec_s0mag[cindex]
+ );
+ fprintf(
+ p_outfile, " COSAV=%15.7lE, RAPRS=%15.7lE\n",
+ vec_cosav[cindex], vec_raprs[cindex]
+ );
+ fprintf(
+ p_outfile, " IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n",
+ vec_tqek1[cindex], vec_tqsk1[cindex]
+ );
+ fprintf(
+ p_outfile, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n",
+ vec_tqek2[cindex], vec_tqsk2[cindex]
+ );
+ } // ci configuration loop
+ if (nsph != 1) {
+ fprintf(
+ p_outfile, " FSAT=(%15.7lE,%15.7lE)\n",
+ real(vec_fsat[jxi]), imag(vec_fsat[jxi])
+ );
+ fprintf(
+ p_outfile, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ vec_qschut[jxi], vec_pschut[jxi], vec_s0magt[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++) {
+ int dir_index = ndirs * jxi + done_dirs;
+ bool goto190 = (done_dirs == 0) && ((jxi > 0) || (jth > 0) || (jph > 0));
+ 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[3 * done_dirs],
+ vec_dir_un[3 * done_dirs + 1],
+ vec_dir_un[3 * done_dirs + 2]
+ );
+ fprintf(
+ p_outfile, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n",
+ vec_dir_uns[3 * done_dirs],
+ vec_dir_uns[3 * done_dirs + 1],
+ vec_dir_uns[3 * done_dirs + 2]
+ );
+ } else {
+ fprintf(
+ p_outfile, " UN=(%12.5lE,%12.5lE,%12.5lE)\n",
+ vec_dir_un[3 * done_dirs],
+ vec_dir_un[3 * done_dirs + 1],
+ vec_dir_un[3 * done_dirs + 2]
+ );
+ }
+ for (int i = 0; i < nsph; i++) {
+ int cindex = jxi * nsph * ndirs + nsph * done_dirs + i;
+ fprintf(p_outfile, " SPHERE %2d\n", i + 1);
+ fprintf(
+ p_outfile, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
+ real(vec_dir_sas11[cindex]),
+ imag(vec_dir_sas11[cindex]),
+ real(vec_dir_sas21[cindex]),
+ imag(vec_dir_sas21[cindex])
+ );
+ fprintf(
+ p_outfile, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
+ real(vec_dir_sas12[cindex]),
+ imag(vec_dir_sas12[cindex]),
+ real(vec_dir_sas22[cindex]),
+ imag(vec_dir_sas22[cindex])
+ );
+ if (jths == 0 && jphs == 0) {
+ fprintf(
+ p_outfile, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
+ vec_dir_fx[jxi * nsph * _num_theta * _num_phi + jth * nsph * _num_phi + jph * nsph + i],
+ vec_dir_fy[jxi * nsph * _num_theta * _num_phi + jth * nsph * _num_phi + jph * nsph + i],
+ vec_dir_fz[jxi * nsph * _num_theta * _num_phi + jth * nsph * _num_phi + jph * nsph + i]
+ );
+ }
+ fprintf(p_outfile, " MULS\n");
+ for (int j = 0; j < 4; j++) {
+ int muls_index = 16 * cindex + 4 * j;
+ fprintf(
+ p_outfile, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_dir_muls[muls_index],
+ vec_dir_muls[muls_index + 1],
+ vec_dir_muls[muls_index + 2],
+ vec_dir_muls[muls_index + 3]
+ );
+ } // j muls loop
+ fprintf(p_outfile, " MULSLR\n");
+ for (int j = 0; j < 4; j++) {
+ int muls_index = 16 * cindex + 4 * j;
+ fprintf(
+ p_outfile, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ vec_dir_mulslr[muls_index],
+ vec_dir_mulslr[muls_index + 1],
+ vec_dir_mulslr[muls_index + 2],
+ vec_dir_mulslr[muls_index + 3]
+ );
+ } // j mulslr loop
+ } // i sphere loop
+ 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 SphereOutputInfo::mpireceive(const mixMPI *mpidata, int pid) {
+ int result = 0;
+ int chk_nsph, chk_inpol, 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_nsph, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(&chk_inpol, 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_nsph == nsph) ? 0 : 1;
+ result += (chk_inpol == inpol) ? 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);
+ if (nsph != 1) {
+ MPI_Recv(vec_fsat + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_qschut + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_pschut + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_s0magt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ }
+
+ // Receive vectors whose sizes depend on configurations and wavelengths
+ MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ offset = (xi1 - _first_xi) * configurations;
+ 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);
+ MPI_Recv(vec_scs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_abs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_exs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_albeds + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_scsrt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_absrt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_exsrt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_fsas + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_qschu + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_pschu + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_s0mag + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_cosav + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_raprs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_tqek1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_tqek2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_tqsk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_tqsk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+ // Receive vectors whose sizes depend on NSPH, directions and wavelengths
+ MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ offset = (xi1 - _first_xi) * nsph * ndirs;
+ 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_muls + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_mulslr + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+ // Receive vectors whose sizes depend on NSPH, incidence directions and wavelengths
+ MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ offset = (xi1 - _first_xi) * nsph * _num_theta * _num_phi;
+ MPI_Recv(vec_dir_fx + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fy + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+ MPI_Recv(vec_dir_fz + offset, 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 SphereOutputInfo::mpisend(const mixMPI *mpidata) {
+ int result = 0;
+ int chunk_size;
+ // Send output metadata for configuration cross-check
+ MPI_Send(&nsph, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(&inpol, 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);
+ if (nsph != 1) {
+ MPI_Send(vec_fsat, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_qschut, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_pschut, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_s0magt, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ }
+
+ // Send vectors whose sizes depend on configurations and scales
+ chunk_size = xi_block_size * configurations;
+ 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);
+ MPI_Send(vec_scs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_abs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_exs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_albeds, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_scsrt, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_absrt, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_exsrt, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_fsas, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_qschu, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_pschu, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_s0mag, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_cosav, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_raprs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_tqek1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_tqek2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_tqsk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_tqsk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+
+ // Send vectors whose sizes depend on NSPH, directions and wavelengths
+ chunk_size = xi_block_size * nsph * ndirs;
+ MPI_Send(&chunk_size, 1, MPI_INT32_T, 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_muls, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_mulslr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+
+ // Send vectors whose sizes depend on NSPH, incidence directions and wavelengths
+ chunk_size = xi_block_size * nsph * _num_theta * _num_phi;
+ MPI_Send(vec_dir_fx, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fy, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ MPI_Send(vec_dir_fz, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+ }
+ return result;
+}
+#endif // MPI_VERSION
+// >>> END OF SphereOutputInfo CLASS IMPLEMENTATION <<<
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 9f8c661698188683e17ea41db885f12f1fd658b3..dafa02def09757897e92c1b8ff336f04388ad1b4 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -51,6 +51,10 @@
#include "../include/TransitionMatrix.h"
#endif
+#ifndef INCLUDE_OUTPUTS_H_
+#include "../include/outputs.h"
+#endif
+
using namespace std;
/*! \brief C++ implementation of SPH
@@ -90,6 +94,7 @@ void sphere(const string& config_file, const string& data_file, const string& ou
}
int s_nsph = sconf->number_of_spheres;
int nsph = gconf->number_of_spheres;
+ int configurations = sconf->configurations;
if (s_nsph == nsph) {
int isq, ibf;
double *argi, *args, *gaps;
@@ -127,6 +132,7 @@ void sphere(const string& config_file, const string& data_file, const string& ou
}
double frx = 0.0, fry = 0.0, frz = 0.0;
double cfmp, cfsp, sfmp, sfsp;
+ const dcomplex cc0 = 0.0 * I + 0.0;
int jw;
int l_max = gconf->l_max;
ParticleDescriptor *c1 = new ParticleDescriptorSphere(gconf, sconf);
@@ -150,48 +156,9 @@ void sphere(const string& config_file, const string& data_file, const string& ou
argi = new double[1];
args = new double[1];
gaps = new double[2];
- FILE *output = fopen((output_path + "/c_OSPH").c_str(), "w");
- fprintf(output, " READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM\n");
- fprintf(
- output,
- " %5d%5d%5d%5d%5d%5d\n",
- nsph,
- l_max,
- in_pol,
- npnt,
- npntts,
- meridional_type
- );
- fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
- fprintf(
- output,
- " %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",
- in_theta_start,
- in_theta_step,
- in_theta_end,
- sc_theta_start,
- sc_theta_step,
- sc_theta_end
- );
- fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
- fprintf(
- output,
- " %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",
- in_phi_start,
- in_phi_step,
- in_phi_end,
- sc_phi_start,
- sc_phi_step,
- sc_phi_end
- );
- fprintf(output, " READ(IR,*)JWTM\n");
- fprintf(output, " %5d\n", jwtm);
- fprintf(output, " READ(ITIN)NSPHT\n");
- fprintf(output, " READ(ITIN)(IOG(I),I=1,NSPH)\n");
- fprintf(output, " READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
- fprintf(output, " READ(ITIN)(XIV(I),I=1,NXI)\n");
- fprintf(output, " READ(ITIN)NSHL(I),ROS(I)\n");
- fprintf(output, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n \n");
+ // FILE *output = fopen((output_path + "/c_OSPH").c_str(), "w");
+ mixMPI *mpidata = new mixMPI(); // Fake MPI configuration
+ SphereOutputInfo *p_output = new SphereOutputInfo(sconf, gconf, mpidata);
double sml = 1.0e-3;
int nth = 0, nph = 0;
if (in_theta_step != 0.0)
@@ -253,7 +220,6 @@ void sphere(const string& config_file, const string& data_file, const string& ou
thdps(l_max, zpv);
double exdc = sconf->exdc;
double exri = sqrt(exdc);
- fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
fstream tppoan;
string tppoan_name = output_path + "/c_TPPOAN";
tppoan.open(tppoan_name.c_str(), ios::binary|ios::out);
@@ -271,9 +237,6 @@ void sphere(const string& config_file, const string& data_file, const string& ou
for (int nsi = 0; nsi < nsph; nsi++)
tppoan.write(reinterpret_cast<char *>(&(c1->iog[nsi])), sizeof(int));
- if (in_pol == 0) fprintf(output, " LIN\n");
- else fprintf(output, " CIRC\n");
- fprintf(output, " \n");
double wp = sconf->wp;
double xip = sconf->xip;
double wn = wp / 3.0e8;
@@ -283,22 +246,24 @@ void sphere(const string& config_file, const string& data_file, const string& ou
int nxi = sconf->number_of_scales;
if (idfc < 0) {
vk = xip * wn;
- fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
- fprintf(output, " \n");
+ p_output->vec_vk[0] = vk;
}
+ int ndirs = nkks;
for (int jxi488 = 0; jxi488 < nxi; jxi488++) {
+ int oindex = 0;
int jxi = jxi488 + 1;
logger->log("INFO: running scale iteration " + to_string(jxi) + " of " + to_string(nxi) + ".\n");
- fprintf(output, "========== JXI =%3d ====================\n", jxi);
double xi = sconf->get_scale(jxi488);
if (idfc >= 0) {
vk = xi * wn;
vkarg = vk;
- fprintf(output, " VK=%15.7lE, XI=%15.7lE\n", xi, vk);
+ p_output->vec_vk[jxi488] = vk;
+ p_output->vec_xi[jxi488] = xi;
} else { // IDFC < 0
vkarg = xi * vk;
sqsfi = 1.0 / (xi * xi);
- fprintf(output, " XI=%15.7lE\n", xi);
+ p_output->vec_vk[jxi488] = vk;
+ p_output->vec_xi[jxi488] = xi;
}
tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
for (int i132 = 0; i132 < nsph; i132++) {
@@ -329,10 +294,10 @@ void sphere(const string& config_file, const string& data_file, const string& ou
int lcalc = 0;
dme(l_max, i, npnt, npntts, vkarg, exdc, exri, c1, jer, lcalc, arg);
if (jer != 0) {
- fprintf(output, " STOP IN DME\n");
- fprintf(output, " AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n", jer, lcalc, real(arg), imag(arg));
+ p_output->vec_ier[jxi] = 1;
+ p_output->lcalc = lcalc;
+ p_output->arg = arg;
tppoan.close();
- fclose(output);
delete sconf;
delete gconf;
delete c1;
@@ -398,55 +363,47 @@ void sphere(const string& config_file, const string& data_file, const string& ou
double sqk = vk * vk * exdc;
aps(zpv, l_max, nsph, c1, sqk, gaps);
rabas(in_pol, l_max, nsph, c1, tqse, tqspe, tqss, tqsps);
+ int last_configuration = 0;
for (int i170 = 0; i170 < nsph; i170++) {
int i = i170 + 1;
if (c1->iog[i170] >= i) {
+ last_configuration++;
+ oindex = jxi488 * configurations + last_configuration - 1;
double albeds = c1->sscs[i170] / c1->sexs[i170];
c1->sqscs[i170] *= sqsfi;
c1->sqabs[i170] *= sqsfi;
c1->sqexs[i170] *= sqsfi;
- fprintf(output, " SPHERE %2d\n", i);
if (c1->nshl[i170] != 1) {
- fprintf(output, " SIZE=%15.7lE\n", c1->vsz[i170]);
+ p_output->vec_sphere_ref_indices[oindex] = cc0;
+ p_output->vec_sphere_sizes[oindex] = c1->vsz[i170];
} else {
- fprintf(
- output,
- " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n",
- c1->vsz[i170],
- real(c1->vkt[i170]),
- imag(c1->vkt[i170])
- );
+ p_output->vec_sphere_ref_indices[oindex] = c1->vkt[i170];
+ p_output->vec_sphere_sizes[oindex] = c1->vsz[i170];
}
- fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
- fprintf(
- output,
- " %14.7lE%15.7lE%15.7lE%15.7lE\n",
- c1->sscs[i170], c1->sabs[i170],
- c1->sexs[i170], albeds
- );
- fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
- fprintf(
- output,
- " %14.7lE%15.7lE%15.7lE\n",
- c1->sqscs[i170], c1->sqabs[i170],
- c1->sqexs[i170]
- );
- fprintf(output, " FSAS=%15.7lE%15.7lE\n", real(c1->fsas[i170]), imag(c1->fsas[i170]));
+ p_output->vec_scs[oindex] = c1->sscs[i170];
+ p_output->vec_abs[oindex] = c1->sabs[i170];
+ p_output->vec_exs[oindex] = c1->sexs[i170];
+ p_output->vec_albeds[oindex] = albeds;
+ p_output->vec_scsrt[oindex] = c1->sqscs[i170];
+ p_output->vec_absrt[oindex] = c1->sqabs[i170];
+ p_output->vec_exsrt[oindex] = c1->sqexs[i170];
+ p_output->vec_fsas[oindex] = c1->fsas[i170];
double csch = 2.0 * vk * sqsfi / c1->gcsv[i170];
s0 = c1->fsas[i170] * exri;
double qschu = csch * imag(s0);
double pschu = csch * real(s0);
double s0mag = cs0 * cabs(s0);
- fprintf(
- output,
- " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
- qschu, pschu, s0mag
- );
+ p_output->vec_qschu[oindex] = qschu;
+ p_output->vec_pschu[oindex] = pschu;
+ p_output->vec_s0mag[oindex] = s0mag;
double rapr = c1->sexs[i170] - gaps[i170];
double cosav = gaps[i170] / c1->sscs[i170];
- fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
- fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 1, tqse[0][i170], tqss[0][i170]);
- fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170], tqss[1][i170]);
+ p_output->vec_cosav[oindex] = cosav;
+ p_output->vec_raprs[oindex] = rapr;
+ p_output->vec_tqek1[oindex] = tqse[0][i170];
+ p_output->vec_tqsk1[oindex] = tqss[0][i170];
+ p_output->vec_tqek2[oindex] = tqse[1][i170];
+ p_output->vec_tqsk2[oindex] = tqss[1][i170];
tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170])), sizeof(double));
tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170])), sizeof(double));
double val = real(tqspe[0][i170]);
@@ -470,19 +427,18 @@ void sphere(const string& config_file, const string& data_file, const string& ou
} // End if iog[i170] >= i
} // i170 loop
if (nsph != 1) {
- fprintf(output, " FSAT=(%15.7lE,%15.7lE)\n", real(tfsas), imag(tfsas));
+ p_output->vec_fsat[jxi488] = tfsas;
double csch = 2.0 * vk * sqsfi / gcs;
s0 = tfsas * exri;
double qschu = csch * imag(s0);
double pschu = csch * real(s0);
double s0mag = cs0 * cabs(s0);
- fprintf(
- output,
- " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
- qschu, pschu, s0mag
- );
+ p_output->vec_qschut[jxi488] = qschu;
+ p_output->vec_pschut[jxi488] = pschu;
+ p_output->vec_s0magt[jxi488] = s0mag;
}
th = th1;
+ int done_dirs = 0;
for (int jth486 = 0; jth486 < nth; jth486++) { // OpenMP parallelizable section
int jth = jth486 + 1;
ph = ph1;
@@ -556,61 +512,51 @@ void sphere(const string& config_file, const string& data_file, const string& ou
tppoan.write(reinterpret_cast<char *>(&(u[1])), sizeof(double));
tppoan.write(reinterpret_cast<char *>(&(u[2])), sizeof(double));
}
- fprintf(
- output,
- "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n",
- jth, jph, jths, jphs
- );
- fprintf(
- output,
- " TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n",
- th, ph, ths, phs
- );
- fprintf(output, " SCAND=%10.3lE\n", scan);
- fprintf(output, " CFMP=%15.7lE, SFMP=%15.7lE\n", cfmp, sfmp);
- fprintf(output, " CFSP=%15.7lE, SFSP=%15.7lE\n", cfsp, sfsp);
+ p_output->vec_dir_scand[done_dirs] = scan;
+ p_output->vec_dir_cfmp[done_dirs] = cfmp;
+ p_output->vec_dir_cfsp[done_dirs] = cfsp;
+ p_output->vec_dir_sfmp[done_dirs] = sfmp;
+ p_output->vec_dir_sfsp[done_dirs] = sfsp;
if (meridional_type >= 0) {
- fprintf(output, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
- fprintf(output, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n", uns[0], uns[1], uns[2]);
+ p_output->vec_dir_un[3 * done_dirs] = un[0];
+ p_output->vec_dir_un[3 * done_dirs + 1] = un[1];
+ p_output->vec_dir_un[3 * done_dirs + 2] = un[2];
+ p_output->vec_dir_uns[3 * done_dirs] = uns[0];
+ p_output->vec_dir_uns[3 * done_dirs + 1] = uns[1];
+ p_output->vec_dir_uns[3 * done_dirs + 2] = uns[2];
} else {
- fprintf(output, " UN=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
+ p_output->vec_dir_un[3 * done_dirs] = un[0];
+ p_output->vec_dir_un[3 * done_dirs + 1] = un[1];
+ p_output->vec_dir_un[3 * done_dirs + 2] = un[2];
+ p_output->vec_dir_uns[3 * done_dirs] = 0.0;
+ p_output->vec_dir_uns[3 * done_dirs + 1] = 0.0;
+ p_output->vec_dir_uns[3 * done_dirs + 2] = 0.0;
}
sscr2(nsph, l_max, vk, exri, c1);
+ last_configuration = 0;
for (int ns226 = 0; ns226 < nsph; ns226++) {
int ns = ns226 + 1;
- fprintf(output, " SPHERE %2d\n", ns);
- fprintf(
- output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
- real(c1->sas[ns226][0][0]), imag(c1->sas[ns226][0][0]),
- real(c1->sas[ns226][1][0]), imag(c1->sas[ns226][1][0])
- );
- fprintf(
- output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
- real(c1->sas[ns226][0][1]), imag(c1->sas[ns226][0][1]),
- real(c1->sas[ns226][1][1]), imag(c1->sas[ns226][1][1])
- );
- if (jths == 1 && jphs == 1)
- fprintf(
- output, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
- frx, fry, frz
- );
+ oindex = jxi488 * nsph * ndirs + nsph * done_dirs + ns226;
+ p_output->vec_dir_sas11[oindex] = c1->sas[ns226][0][0];
+ p_output->vec_dir_sas21[oindex] = c1->sas[ns226][1][0];
+ p_output->vec_dir_sas12[oindex] = c1->sas[ns226][0][1];
+ p_output->vec_dir_sas22[oindex] = c1->sas[ns226][1][1];
+ p_output->vec_dir_fx[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frx;
+ p_output->vec_dir_fy[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = fry;
+ p_output->vec_dir_fz[jxi488 * nsph * nth * nph + nsph * nph * (jth - 1) + nsph * (jph - 1) + ns226] = frz;
for (int i225 = 0; i225 < 16; i225++) c1->vint[i225] = c1->vints[ns226][i225];
mmulc(c1->vint, cmullr, cmul);
- fprintf(output, " MULS\n ");
for (int imul = 0; imul < 4; imul++) {
+ int muls_index = 16 * jxi488 * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;
for (int jmul = 0; jmul < 4; jmul++) {
- fprintf(output, "%15.7lE", cmul[imul][jmul]);
+ p_output->vec_dir_muls[muls_index + jmul] = cmul[imul][jmul];
}
- if (imul < 3) fprintf(output, "\n ");
- else fprintf(output, "\n");
}
- fprintf(output, " MULSLR\n ");
for (int imul = 0; imul < 4; imul++) {
+ int muls_index = 16 * jxi488 * nsph * ndirs + 16 * nsph * done_dirs + 4 * imul;
for (int jmul = 0; jmul < 4; jmul++) {
- fprintf(output, "%15.7lE", cmullr[imul][jmul]);
+ p_output->vec_dir_mulslr[muls_index + jmul] = cmullr[imul][jmul];
}
- if (imul < 3) fprintf(output, "\n ");
- else fprintf(output, "\n");
}
for (int vi = 0; vi < 16; vi++) {
double value = real(c1->vint[vi]);
@@ -625,6 +571,7 @@ void sphere(const string& config_file, const string& data_file, const string& ou
}
} // ns226 loop
if (meridional_type < 1) phs += sc_phi_step;
+ done_dirs++;
} // jphs480 loop
if (meridional_type <= 1) thsl += sc_theta_step;
} // jths482 loop
@@ -632,13 +579,18 @@ void sphere(const string& config_file, const string& data_file, const string& ou
} // jph484 loop on elevation
th += in_theta_step;
} // jth486 loop on azimuth
+ p_output->vec_jxi[jxi488] = jxi488 + 1;
logger->log("INFO: finished scale iteration " + to_string(jxi) + " of " + to_string(nxi) + ".\n");
} //jxi488 loop on scales
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
logger->err("ERROR: failed to open TPPOAN file.\n");
}
- fclose(output);
+ // fclose(output);
+ p_output->write(output_path + "/c_OSPH.hd5", "HDF5");
+ p_output->write(output_path + "/c_OSPH", "LEGACY");
+ delete p_output;
+ delete mpidata;
delete c1;
for (int zi = l_max - 1; zi > -1; zi--) {
for (int zj = 0; zj < 3; zj++) {