Implement dynamic memory management for spherical harmonics

src/include/sph_subs.h

@@ -954,7 +954,9 @@ void sphar(
 		double cosrth, double sinrth, double cosrph, double sinrph,
 		int ll, std::complex<double> *ylm
 ) {
-	double sinrmp[40], cosrmp[40], plegn[861];
+	int rmp_size = ll;
+	int plegn_size = (ll + 1) * ll / 2 + ll + 1;
+	double sinrmp[rmp_size], cosrmp[rmp_size], plegn[plegn_size];
 	double four_pi = 8.0 * acos(0.0);
 	double pi4irs = 1.0 / sqrt(four_pi);
 	double x = cosrth;
@@ -1240,7 +1242,8 @@ void wmamp(
 		int lm, int idot, int nsph, double *arg, double *u, double *up,
 		double *un, C1 *c1
 ) {
-	std::complex<double> *ylm = new std::complex<double>[1682];
+	int ylm_size = (lm + 1) * (lm + 1) + 1;
+	std::complex<double> *ylm = new std::complex<double>[ylm_size];
 	int nlmp = lm * (lm + 2) + 2;
 	ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
 	if (idot != 0) {
@@ -1295,7 +1298,8 @@ void wmasp(
 		double *ups, double *uns, int isq, int ibf, int inpol, int lm, int idot,
 		int nsph, double *argi, double *args, C1 *c1
 ) {
-	std::complex<double> *ylm = new std::complex<double>[1682];
+	int ylm_size = (lm + 1) * (lm + 1) + 1;
+	std::complex<double> *ylm = new std::complex<double>[ylm_size];
 	int nlmp = lm * (lm + 2) + 2;
 	ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
 	if (idot != 0) {

src/sphere/sphere.cpp

@@ -2,8 +2,12 @@
 #include <fstream>
 #include <string>
 #include <complex>
+#ifndef INCLUDE_CONFIGURATION_H_
 #include "../include/Configuration.h"
+#endif
+#ifndef INCLUDE_SPH_SUBS_H_
 #include "../include/sph_subs.h"
+#endif
 
 using namespace std;
 
@@ -15,13 +19,13 @@ void sphere() {
 	double *argi, *args, *gaps;
 	double th, ph;
 	printf("INFO: making legacy configuration ...\n");
-	ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/sphere/DEDFB");
-	conf->write_formatted("c_OEDFB");
-	conf->write_binary("c_TEDF");
+	ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("DEDFB_sph");
+	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");
-	GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/sphere/DSPH");
+	ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF_sph");
+	GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("DSPH");
 	if (sconf->number_of_spheres == gconf->number_of_spheres) {
 		int isq, ibf;
 		double cost, sint, cosp, sinp;
@@ -49,18 +53,10 @@ void sphere() {
 		complex<double> *vint = new complex<double>[16];
 		int jw;
 		int nsph = gconf->number_of_spheres;
-		C1 *c1 = new C1(nsph, gconf->l_max);
+		C1 *c1 = new C1(nsph, gconf->l_max, sconf->nshl_vec, sconf->iog_vec);
 		for (int i = 0; i < nsph; i++) {
-			c1->iog[i] = sconf->iog_vec[i];
-			c1->nshl[i] = sconf->nshl_vec[i];
 			c1->ros[i] = sconf->radii_of_spheres[i];
 		}
-		for (int i = 0; i < gconf->l_max; i++) {
-			c1->rmi[i][0] = complex<double>(0.0, 0.0);
-			c1->rmi[i][1] = complex<double>(0.0, 0.0);
-			c1->rei[i][0] = complex<double>(0.0, 0.0);
-			c1->rei[i][1] = complex<double>(0.0, 0.0);
-		}
 		C2 *c2 = new C2(nsph, 5, gconf->npnt, gconf->npntts);
 		argi = new double[1];
 		args = new double[1];
@@ -172,7 +168,7 @@ void sphere() {
 		double exri = sqrt(sconf->exdc);
 		fprintf(output, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
 		fstream tppoan;
-		tppoan.open("c_TPPOAN", ios::binary|ios::out);
+		tppoan.open("c_TPPOAN_sph", ios::binary|ios::out);
 		if (tppoan.is_open()) {
 			int imode = 10;
 			tppoan.write(reinterpret_cast<char *>(&imode), sizeof(int));
@@ -250,7 +246,7 @@ void sphere() {
 					// This is the condition that writes the transition matrix to output.
 					int is = 1111;
 					fstream ttms;
-					ttms.open("c_TTMS", ios::binary | ios::out);
+					ttms.open("c_TTMS_sph", 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));
@@ -550,9 +546,12 @@ void sphere() {
 		}
 		delete[] cmul;
 		delete[] cmullr;
+		printf("Done.\n");
 	} else { // NSPH mismatch between geometry and scatterer configurations.
 		throw UnrecognizedConfigurationException(
 				"Inconsistent geometry and scatterer configurations."
 		);
 	}
+	delete sconf;
+	delete gconf;
 }