Implement command-line options for input and output

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());
 }

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;

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;

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;

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."