diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp index a093e0c9a5eda818f12701d3fc355e5b7d32cb79..b685453d6aaf15996a52419508ee5525cfc7a115 100644 --- a/src/cluster/cluster.cpp +++ b/src/cluster/cluster.cpp @@ -16,16 +16,19 @@ using namespace std; * */ -//! \brief C++ implementation of CLU -void cluster() { - printf("INFO: making legacy configuration ...\n"); - ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/cluster/DEDFB"); - conf->write_formatted("c_OEDFB_clu"); - conf->write_binary("c_TEDF_clu"); - delete conf; - printf("INFO: reading binary configuration ...\n"); - ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF_clu"); - GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/cluster/DCLU"); +/*! \brief C++ implementation of CLU + * + * \param config_file: `string` Name of the configuration file. + * \param data_file: `string` Name of the input data file. + * \param output_path: `string` Directory to write the output files in. + */ +void cluster(string config_file, string data_file, string output_path) { + printf("INFO: making legacy configuration..."); + ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(config_file); + sconf->write_formatted(output_path + "/c_OEDFB_clu"); + sconf->write_binary(output_path + "/c_TEDF_clu"); + GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(data_file); + printf(" done.\n"); if (sconf->number_of_spheres == gconf->number_of_spheres) { // Shortcuts to variables stored in configuration objects int nsph = gconf->number_of_spheres; @@ -78,7 +81,7 @@ void cluster() { C1_AddOns *c1ao = new C1_AddOns(c4); // End of add-ons initialization C6 *c6 = new C6(c4->lmtpo); - FILE *output = fopen("c_OCLU", "w"); + FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w"); int jer = 0, lcalc = 0; complex<double> arg(0.0, 0.0), ccsam(0.0, 0.0); int max_ici = 0; @@ -226,7 +229,8 @@ void cluster() { double vk = 0.0; // NOTE: Not really sure it should be initialized at 0 fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri); fstream tppoan; - tppoan.open("c_TPPOAN", ios::out | ios::binary); + string tppoan_name = output_path + "/c_TPPOAN"; + tppoan.open(tppoan_name.c_str(), ios::out | ios::binary); if (tppoan.is_open()) { tppoan.write(reinterpret_cast<char *>(&iavm), sizeof(int)); tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int)); @@ -244,6 +248,7 @@ void cluster() { fprintf(output, " \n"); } for (int jxi488 = 1; jxi488 <= nxi; jxi488++) { + printf("INFO: running scale iteration...\n"); int jaw = 1; fprintf(output, "========== JXI =%3d ====================\n", jxi488); double xi = sconf->scale_vec[jxi488 - 1]; @@ -293,20 +298,25 @@ void cluster() { } if (jer != 0) break; } // i132 loop + printf("INFO: initializing matrix..."); cms(am, c1, c1ao, c4, c6); - ccsam = summat(am, 960, 960); - printf("DEBUG: after CMS CCSAM = (%lE,%lE)\n", ccsam.real(), ccsam.imag()); + printf(" done.\n"); + //ccsam = summat(am, 960, 960); + //printf("DEBUG: after CMS CCSAM = (%lE,%lE)\n", ccsam.real(), ccsam.imag()); int ndit = 2 * nsph * c4->nlim; + printf("INFO: inverting matrix..."); lucin(am, mxndm, ndit, jer); - ccsam = summat(am, 960, 960); - printf("DEBUG: after LUCIN CCSAM = (%lE,%lE)\n", ccsam.real(), ccsam.imag()); + printf(" done.\n"); + //ccsam = summat(am, 960, 960); + //printf("DEBUG: after LUCIN CCSAM = (%lE,%lE)\n", ccsam.real(), ccsam.imag()); if (jer != 0) break; // jxi488 loop: goes to memory clean ztm(am, c1, c1ao, c4, c6, c9); if (sconf->idfc >= 0) { if (jxi488 == gconf->jwtm) { int is = 1; fstream ttms_file; - ttms_file.open("c_TTMS", ios::out | ios::binary); + string ttms_name = output_path + "/c_TTMS"; + ttms_file.open(ttms_name.c_str(), ios::out | ios::binary); if (ttms_file.is_open()) { ttms_file.write(reinterpret_cast<char *>(&is), sizeof(int)); ttms_file.write(reinterpret_cast<char *>(&lm), sizeof(int)); @@ -340,7 +350,7 @@ void cluster() { double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0; std::complex<double> s0(0.0, 0.0); scr0(vk, exri, c1, c1ao, c3, c4); - printf("DEBUG: after SCR0 TFSAS = (%lE, %lE)\n", c3->tfsas.real(), c3->tfsas.imag()); + //printf("DEBUG: after SCR0 TFSAS = (%lE, %lE)\n", c3->tfsas.real(), c3->tfsas.imag()); double sqk = vk * vk * sconf->exdc; aps(zpv, c4->li, nsph, c1, sqk, gaps); rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps); @@ -533,7 +543,7 @@ void cluster() { complex<double> s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri; double qschum = s0m.imag() * csch; double pschum = s0m.real() * csch; - double s0magm = sqrt((s0m.real() + s0m.imag()) * (s0m.real() - s0m.imag())) * cs0; + double s0magm = abs(s0m) * cs0; double rfinrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].real() / c3->tfsas.real(); double extcrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag() / c3->tfsas.imag(); if (inpol == 0) { @@ -871,7 +881,7 @@ void cluster() { } // jph484 loop th += thstp; } // jth486 loop - printf("INFO: done jxi488 iteration.\n"); + printf("INFO: done scale.\n"); } // jxi488 loop tppoan.close(); } else { // In case TPPOAN could not be opened. Should never happen. @@ -958,5 +968,5 @@ void cluster() { } delete sconf; delete gconf; - printf("Done.\n"); + printf("Finished: output written to %s.\n", (output_path + "/c_OCLU").c_str()); } diff --git a/src/cluster/np_cluster.cpp b/src/cluster/np_cluster.cpp index 7b0e087a71e00ec4b63f7bb8446d5d0629793447..e9a85e79d16d82e359654456a5d0345396c08fa0 100644 --- a/src/cluster/np_cluster.cpp +++ b/src/cluster/np_cluster.cpp @@ -4,7 +4,7 @@ #include <cstdio> #include <string> #ifndef INCLUDE_CONFIGURATION_H_ -#include "include/Configuration.h" +#include "../include/Configuration.h" #endif using namespace std; diff --git a/src/include/Configuration.h b/src/include/Configuration.h index 12e3d779f562e8fe8116dd6fb549c24f724479c6..e784fa38bf0bd7b2f681102f0cbba5bbb7628028 100644 --- a/src/include/Configuration.h +++ b/src/include/Configuration.h @@ -65,8 +65,8 @@ public: */ class GeometryConfiguration { //! Temporary work-around to allow cluster() and sphere() peeking in. - friend void cluster(); - friend void sphere(); + friend void cluster(std::string, std::string, std::string); + friend void sphere(std::string, std::string, std::string); protected: //! \brief Number of spherical components. int number_of_spheres; @@ -191,8 +191,8 @@ public: */ class ScattererConfiguration { //! Temporary work-around to allow cluster() and sphere() peeking in. - friend void cluster(); - friend void sphere(); + friend void cluster(std::string, std::string, std::string); + friend void sphere(std::string, std::string, std::string); protected: //! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x LAYERS]. std::complex<double> ***dc0_matrix; diff --git a/src/sphere/np_sphere.cpp b/src/sphere/np_sphere.cpp index 4522e0282b5c2b3c19cbdd01a6112bc87a290c59..cec49e8e78ad43edc34426830c41ecf7b42c100f 100644 --- a/src/sphere/np_sphere.cpp +++ b/src/sphere/np_sphere.cpp @@ -4,7 +4,7 @@ #include <cstdio> #include <string> #ifndef INCLUDE_CONFIGURATION_H_ -#include "include/Configuration.h" +#include "../include/Configuration.h" #endif using namespace std; diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp index 7565ed594e2f9caca942a4d04d350cf620b4bab8..e5adf0421b10dbbbe737d4da361e834fee6766d6 100644 --- a/src/sphere/sphere.cpp +++ b/src/sphere/sphere.cpp @@ -11,18 +11,18 @@ using namespace std; -//! \brief C++ implementation of SPH -void sphere() { +/*! \brief C++ implementation of SPH + * + * \param config_file: `string` Name of the configuration file. + * \param data_file: `string` Name of the input data file. + * \param output_path: `string` Directory to write the output files in. + */ +void sphere(string config_file, string data_file, string output_path) { complex<double> arg, s0, tfsas; double th, ph; printf("INFO: making legacy configuration ...\n"); - ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/sphere/DEDFB"); - conf->write_formatted("c_OEDFB_sph"); - conf->write_binary("c_TEDF_sph"); - delete conf; - printf("INFO: reading binary configuration ...\n"); - ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF_sph"); - GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/sphere/DSPH"); + ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(config_file); + GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(data_file); if (sconf->number_of_spheres == gconf->number_of_spheres) { int isq, ibf; double *argi, *args, *gaps; @@ -71,7 +71,7 @@ void sphere() { argi = new double[1]; args = new double[1]; gaps = new double[2]; - FILE *output = fopen("c_OSPH", "w"); + FILE *output = fopen((output_path + "/c_OSPH").c_str(), "w"); fprintf(output, " READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM\n"); fprintf( output, @@ -179,7 +179,8 @@ void sphere() { double exri = sqrt(sconf->exdc); fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri); fstream tppoan; - tppoan.open("c_TPPOAN_sph", ios::binary|ios::out); + string tppoan_name = output_path + "/c_TPPOAN_sph"; + tppoan.open(tppoan_name.c_str(), ios::binary|ios::out); if (tppoan.is_open()) { int imode = 10; tppoan.write(reinterpret_cast<char *>(&imode), sizeof(int)); @@ -205,6 +206,7 @@ void sphere() { fprintf(output, " \n"); } for (int jxi488 = 0; jxi488 < sconf->number_of_scales; jxi488++) { + printf("INFO: running scale iteration...\n"); int jxi = jxi488 + 1; fprintf(output, "========== JXI =%3d ====================\n", jxi); double xi = sconf->scale_vec[jxi488]; @@ -259,7 +261,8 @@ void sphere() { // This is the condition that writes the transition matrix to output. int is = 1111; fstream ttms; - ttms.open("c_TTMS_sph", ios::binary | ios::out); + string ttms_name = output_path + "/c_TTMS_sph"; + ttms.open(ttms_name.c_str(), ios::binary | ios::out); if (ttms.is_open()) { ttms.write(reinterpret_cast<char *>(&is), sizeof(int)); ttms.write(reinterpret_cast<char *>(&(gconf->l_max)), sizeof(int)); @@ -281,7 +284,7 @@ void sphere() { double cs0 = 0.25 * vk * vk * vk / half_pi; //printf("DEBUG: cs0 = %lE\n", cs0); sscr0(tfsas, nsph, gconf->l_max, vk, exri, c1); - printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag()); + //printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag()); double sqk = vk * vk * sconf->exdc; aps(zpv, gconf->l_max, nsph, c1, sqk, gaps); rabas(gconf->in_pol, gconf->l_max, nsph, c1, tqse, tqspe, tqss, tqsps); @@ -519,6 +522,7 @@ void sphere() { } // jph484 loop on elevation th += gconf->in_theta_step; } // jth486 loop on azimuth + printf("INFO: done scale.\n"); } //jxi488 loop on scales tppoan.close(); } else { // In case TPPOAN could not be opened. Should never happen. @@ -568,7 +572,7 @@ void sphere() { delete[] tqss; delete[] tqspe; delete[] tqsps; - printf("Done.\n"); + printf("Finished: output written to %s.\n", (output_path + "/c_OSPH").c_str()); } else { // NSPH mismatch between geometry and scatterer configurations. throw UnrecognizedConfigurationException( "Inconsistent geometry and scatterer configurations."