Extend conversion of CLU to C++ (685d3e1b) · Commits · Giacomo Mulas / NP_TMcode

src/cluster/cluster.cpp

+75 −6

Original line number	Diff line number	Diff line
		@@ -19,13 +19,13 @@ using namespace std;
		//! \brief C++ implementation of CLU
		void cluster() {
		printf("INFO: making legacy configuration ...\n");
		ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("DEDFB_clu");
		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("DCLU");
		GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/cluster/DCLU");
		if (sconf->number_of_spheres == gconf->number_of_spheres) {
		// Shortcuts to variables stored in configuration objects
		int nsph = gconf->number_of_spheres;
		@@ -93,6 +93,19 @@ void cluster() {
		for (int ai = 0; ai < mxndm; ai++) am[ai] = new complex<double>[mxndm];
		const int ndi = c4->nsph * c4->nlim;
		C9 c9 = new C9(ndi, c4->nlem, 2 ndi, 2 * c4->nlem);
		double *gaps = new double[nsph]();
		double tqse, tqss;
		complex<double> tqspe, tqsps;
		tqse = new double*[2];
		tqspe = new complex<double>*[2];
		tqss = new double*[2];
		tqsps = new complex<double>*[2];
		for (int ti = 0; ti < 2; ti++) {
		tqse[ti] = new double[nsph]();
		tqspe[ti] = new complex<double>[nsph]();
		tqss[ti] = new double[nsph]();
		tqsps[ti] = new complex<double>[nsph]();
		}
		// End of global variables for CLU
		fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
		fprintf(output, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
		@@ -155,9 +168,7 @@ void cluster() {
		for (int zj = 0; zj < 3; zj++) {
		zpv[zi][zj] = new double*[2];
		for (int zk = 0; zk < 2; zk++) {
		zpv[zi][zj][zk] = new double[2];
		zpv[zi][zj][zk][0] = 0.0;
		zpv[zi][zj][zk][1] = 0.0;
		zpv[zi][zj][zk] = new double[2]();
		}
		}
		}
		@@ -271,7 +282,54 @@ void cluster() {
		}
		}
		// label 156: continue from here

		if (inpol == 0) {
		fprintf(output, " LIN\n");
		} else { // label 158
		fprintf(output, " CIRC\n");
		}
		// label 160
		double cs0 = 0.25 * vk * vk * vk / acos(0.0);
		double csch;
		scr0(vk, exri, c1, c1ao, c3, c4);
		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);
		if (c4->li != c4->le) fprintf(output, " SPHERES; LMX=LI\n");
		for (int i170 = 1; i170 <= nsph; i170++) {
		if (c1->iog[i170 - 1] >= i170) {
		int i = i170 - 1;
		double albeds = c1->sscs[i] / c1->sexs[i];
		c1->sqscs[i] *= sqsfi;
		c1->sqabs[i] *= sqsfi;
		c1->sqexs[i] *= sqsfi;
		fprintf(output, " SPHERE %2d\n", i170);
		if (c1->nshl[i] != 1) {
		fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i]);
		} else { // label 162
		fprintf(output, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE,%15.7lE\n", c2->vsz[i], c2->vkt[i].real(), c2->vkt[i].imag());
		}
		// label 164
		fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
		fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", c1->sscs[i], c1->sabs[i], c1->sexs[i], albeds);
		fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
		fprintf(output, " %14.7lE%15.7lE%15.7lE\n", c1->sqscs[i], c1->sqabs[i], c1->sqexs[i]);
		fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i].real(), c1->fsas[i].imag());
		csch = 2.0 * vk * sqsfi / c1->gcsv[i];
		std::complex<double> s0 = c1->fsas[i] * exri;
		double qschu = s0.imag() * csch;
		double pschu = s0.real() * csch;
		double s0mag = sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag())) * cs0;
		fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
		double rapr = c1->sexs[i] - gaps[i];
		double cosav = gaps[i] / c1->sscs[i];
		fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
		fprintf(output, " IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[0][i], tqss[0][i]);
		fprintf(output, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[1][i], tqss[1][i]);
		}
		} // i170 loop
		fprintf(output, " FSAT=%15.7lE%15.7lE\n", c3->tfsas.real(), c3->tfsas.imag());
		printf("INFO: done jxi488 iteration.\n");
		} // jxi488 loop
		tppoan.close();
		} else { // In case TPPOAN could not be opened. Should never happen.
		@@ -296,6 +354,17 @@ void cluster() {
		delete[] zpv;
		for (int ai = mxndm - 1; ai > -1; ai--) delete[] am[ai];
		delete[] am;
		delete[] gaps;
		for (int ti = 1; ti > -1; ti--) {
		delete[] tqse[ti];
		delete[] tqss[ti];
		delete[] tqspe[ti];
		delete[] tqsps[ti];
		}
		delete[] tqse;
		delete[] tqss;
		delete[] tqspe;
		delete[] tqsps;
		} else { // NSPH mismatch between geometry and scatterer configurations.
		throw UnrecognizedConfigurationException(
		"Inconsistent geometry and scatterer configurations."

src/include/clu_subs.h

+58 −1

Original line number	Diff line number	Diff line
		@@ -21,14 +21,19 @@
		#include <complex>

		// >>> DECLARATION OF SPH_SUBS <<<
		extern void aps(double ***zpv, int li, int nsph, C1 c1, double sqk, double *gaps);
		extern std::complex<double> dconjg(std::complex<double> value);
		extern void dme(
		int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
		C1 c1, C2 c2, int &jer, int &lcalc, std::complex<double> &arg
		);
		extern void sphar(double cth, double sth, double cph, double sph, int lm, std::complex<double> *ylm);
		extern void rabas(
		int inpol, int li, int nsph, C1 c1, double tqse, std::complex<double> *tqspe,
		double tqss, std::complex<double> tqsps
		);
		extern void rbf(int n, double x, int &nm, double sj[]);
		extern void rnf(int n, double x, int &nm, double sy[]);
		extern void sphar(double cth, double sth, double cph, double sph, int lm, std::complex<double> *ylm);
		extern void thdps(int lm, double ****zpv);
		// >>> END OF SPH_SUBS DECLARATION <<<
		void logmat(std::complex<double> **mat, int rows, int cols);
		@@ -910,6 +915,58 @@ void r3j000(int j2, int j3, C6 *c6) {
		}
		}

		/*! \brief C++ porting of SCR0
		*
		* \param vk: `double` QUESTION: definition?
		* \param exri: `double` External medium refractive index. QUESTION: correct?
		* \param c1: `C1 *` Pointer to a C1 instance.
		* \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.
		* \param c3: `C3 *` Pointer to a C3 instance.
		* \param c4: `C4 *` Pointer to a C4 structure.
		*/
		void scr0(double vk, double exri, C1 c1, C1_AddOns c1ao, C3 c3, C4 c4) {
		const std::complex<double> cc0(0.0, 0.0);
		double exdc = exri * exri;
		double ccs = 4.0 * acos(0.0) / (vk * vk);
		double cccs = ccs / exdc;
		std::complex<double> sum21, rm, re, csam;
		csam = -(ccs / (exri * vk)) * std::complex<double>(0.0, 0.5);
		double scs = 0.0, ecs = 0.0, acs = 0.0;
		c3->tfsas = cc0;
		for (int i14 = 1; i14 <= c4->nsph; i14++) {
		int iogi = c1->iog[i14 - 1];
		if (iogi >= i14) {
		double sums = 0.0;
		sum21 = cc0;
		for (int l10 = 1; l10 <= c4->li; l10++) {
		double fl = 1.0 * (l10 + l10 + 1);
		rm = 1.0 / c1->rmi[l10 - 1][i14 - 1];
		re = 1.0 / c1->rei[l10 - 1][i14 - 1];
		double rvalue = (dconjg(rm) * rm + dconjg(re) * re).real() * fl;
		sums += rvalue;
		sum21 += ((rm + re) * fl);
		} // l10 loop
		sum21 *= -1.0;
		double scasec = cccs * sums;
		double extsec = -cccs * sum21.real();
		double abssec = extsec - scasec;
		c1->sscs[i14 - 1] = scasec;
		c1->sexs[i14 - 1] = extsec;
		c1->sabs[i14 - 1] = abssec;
		double gcss = c1->gcsv[i14 - 1];
		c1->sqscs[i14 - 1] = scasec / gcss;
		c1->sqexs[i14 - 1] = extsec / gcss;
		c1->sqabs[i14 - 1] = abssec / gcss;
		c1->fsas[i14 - 1] = sum21 * csam;
		}
		// label 12
		scs += c1->sscs[iogi - 1];
		ecs += c1->sexs[iogi - 1];
		acs += c1->sabs[iogi - 1];
		c3->tfsas += c1->fsas[iogi - 1];
		} // i14 loop
		}

		/*! \brief C++ porting of STR
		*
		* \param rcf: `double **` Matrix of double.