Finish porting of CLU.F to C++ (5473322a) · Commits · Giacomo Mulas / NP_TMcode

src/cluster/cluster.cpp

+317 −3

Original line number	Diff line number	Diff line
		@@ -94,6 +94,8 @@ void cluster() {
		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 *tqev = new double[3]();
		double *tqsv = new double[3]();
		double tqse, tqss, tqce, tqcs;
		complex<double> tqspe, tqsps, tqcpe, tqcps;
		tqse = new double*[2];
		@@ -114,6 +116,7 @@ void cluster() {
		tqcs[ti] = new double[3]();
		tqcps[ti] = new complex<double>[3]();
		}
		double *gapv = new double[3]();
		complex<double> gapp, gappm;
		double gap, gapm;
		gapp = new complex<double>*[3];
		@@ -140,11 +143,16 @@ void cluster() {
		double *args = new double[1]();
		double *duk = new double[3]();
		double cextlr, cext;
		double cmullr, cmul;
		cextlr = new double*[4];
		cext = new double*[4];
		cmullr = new double*[4];;
		cmul = new double*[4];
		for (int ci = 0; ci < 4; ci++) {
		cextlr[ci] = new double[4]();
		cext[ci] = new double[4]();
		cmullr[ci] = new double[4]();
		cmul[ci] = new double[4]();
		}
		int isq, ibf;
		double scan, cfmp, sfmp, cfsp, sfsp;
		@@ -558,12 +566,311 @@ void cluster() {
		fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
		fprintf(output, " Fk=%15.7lE\n", fz);
		} // ilr210 loop
		// from RMBRIF
		double rmbrif = (c1ao->fsacm[0][0].real() - c1ao->fsacm[1][1].real()) / c1ao->fsacm[0][0].real();
		double rmdchr = (c1ao->fsacm[0][0].imag() - c1ao->fsacm[1][1].imag()) / c1ao->fsacm[0][0].imag();
		fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rmbrif);
		fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rmdchr);
		}
		// label 212
		} // jphs loop
		} // jths loop
		fprintf(output, "******** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ******************\n", jth486, jph484, 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);
		if (isam >= 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]);
		} else { // label 214
		fprintf(output, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", un[0], un[1], un[2]);
		}
		// label 220
		if (inpol == 0) {
		fprintf(output, " LIN\n");
		} else { // label 222
		fprintf(output, " CIRC\n");
		}
		// label 224
		scr2(vk, vkarg, exri, duk, c1, c1ao, c3, c4);
		if (c4->li != c4->le) fprintf(output, " SPHERES; LMX=MIN0(LI,LE)\n");
		for (int i226 = 1; i226 <= nsph; i226++) {
		if (c1->iog[i226 - 1] >= i226) {
		fprintf(output, " SPHERE %2d\n", i226);
		fprintf(
		output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
		c1->sas[i226 - 1][0][0].real(), c1->sas[i226 - 1][0][0].imag(),
		c1->sas[i226 - 1][1][0].real(), c1->sas[i226 - 1][1][0].imag()
		);
		fprintf(
		output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
		c1->sas[i226 - 1][0][1].real(), c1->sas[i226 - 1][0][1].imag(),
		c1->sas[i226 - 1][1][1].real(), c1->sas[i226 - 1][1][1].imag()
		);
		for (int j225 = 0; j225 < 16; j225++) { // QUESTION: check that 16 is a fixed dimension
		c1ao->vint[j225] = c1ao->vints[i226 - 1][j225];
		} // j225 loop
		mmulc(c1ao->vint, cmullr, cmul);
		fprintf(output, " MULS\n");
		for (int i1 = 0; i1 < 4; i1++) {
		fprintf(
		output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
		cmul[i1][0], cmul[i1][1], cmul[i1][2], cmul[i1][3]
		);
		} // i1 loop
		fprintf(output, " MULSLR\n");
		for (int i1 = 0; i1 < 4; i1++) {
		fprintf(
		output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
		cmullr[i1][0], cmullr[i1][1], cmullr[i1][2], cmullr[i1][3]
		);
		} // i1 loop
		}
		} // i226 loop
		fprintf(
		output, " SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n",
		c3->tsas[0][0].real(), c3->tsas[0][0].imag(),
		c3->tsas[1][0].real(), c3->tsas[1][0].imag()
		);
		fprintf(
		output, " SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n",
		c3->tsas[0][1].real(), c3->tsas[0][1].imag(),
		c3->tsas[1][1].real(), c3->tsas[1][1].imag()
		);
		fprintf(output, " CLUSTER\n");
		pcros(vk, exri, c1, c1ao, c4);
		mextc(vk, exri, c1ao->fsac, cextlr, cext);
		mmulc(c1ao->vint, cmullr, cmul);
		if (jw != 0) {
		jw = 0;
		// Some implicit loops writing to binary.
		for (int i = 0; i < 4; i++) {
		for (int j = 0; j < 4; j++) {
		double value = cext[i][j];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		}
		for (int i = 0; i < 2; i++) {
		double value = c1ao->scsc[i];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->scscp[i].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->scscp[i].imag();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->ecsc[i];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->ecscp[i].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->ecscp[i].imag();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		for (int i = 0; i < 3; i++) {
		for (int j = 0; j < 2; j++) {
		double value = gap[i][j];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = gapp[i][j].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = gapp[i][j].imag();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		}
		for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 3; j++) {
		double value = tqce[i][j];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = tqcpe[i][j].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = tqcpe[i][j].imag();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		}
		for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 3; j++) {
		double value = tqcs[i][j];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = tqcps[i][j].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = tqcps[i][j].imag();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		}
		for (int i = 0; i < 3; i++) {
		double value = u[i];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = up[i];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = un[i];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		}
		// label 254
		for (int i = 0; i < 16; i++) {
		double value = c1ao->vint[i].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->vint[i].imag();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		for (int i = 0; i < 4; i++) {
		for (int j = 0; j < 4; j++) {
		double value = cmul[i][j];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		}
		int jlr = 2;
		for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
		int ipol = (ilr290 % 2 == 0) ? 1 : -1;
		if (ilr290 == 2) jlr = 1;
		double extsec = c1ao->ecsc[ilr290 - 1];
		double qext = extsec * sqsfi / c3->gcs;
		double extrat = extsec / c3->ecs;
		double scasec = c1ao->scsc[ilr290 - 1];
		double albedc = scasec / extsec;
		double qsca = scasec * sqsfi / c3->gcs;
		double scarat = scasec / c3->scs;
		double abssec = extsec - scasec;
		double qabs = abssec * sqsfi / c3->gcs;
		double absrat = 1.0;
		double ratio = c3->acs / c3->ecs;
		if (ratio < -1.0e-6 \|\| ratio > 1.0e-6) absrat = abssec / c3->acs;
		s0 = c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
		double qschu = s0.imag() * csch;
		double pschu = s0.real() * csch;
		s0mag = sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag())) * cs0;
		double refinr = c1ao->fsac[ilr290 - 1][ilr290 - 1].real() / c3->tfsas.real();
		double extcor = c1ao->fsac[ilr290 - 1][ilr290 - 1].imag() / c3->tfsas.imag();
		if (inpol == 0) {
		fprintf(output, " LIN %2d\n", ipol);
		} else { // label 273
		fprintf(output, " CIRC %2d\n", ipol);
		}
		// label 275
		fprintf(output, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");
		fprintf(
		output, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
		scasec, abssec, extsec, albedc
		);
		fprintf(output, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");
		fprintf(
		output, " %14.7lE%15.7lE%15.7lE\n",
		qsca, qabs, qext
		);
		fprintf(output, " ---- SCCRT --- ABCRT --- EXCRT ----\n");
		fprintf(
		output, " %14.7lE%15.7lE%15.7lE\n",
		scarat, absrat, extrat
		);
		fprintf(
		output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
		ilr290, ilr290, c1ao->fsac[ilr290 - 1][ilr290 - 1].real(), c1ao->fsac[ilr290 - 1][ilr290 - 1].imag(),
		jlr, ilr290, c1ao->fsac[jlr - 1][ilr290 - 1].real(), c1ao->fsac[jlr - 1][ilr290 - 1].imag()
		);
		fprintf(
		output, " SAC(%1d,%1d)=%15.7lE%15.7lE SAC(%1d,%1d)=%15.7lE%15.7lE\n",
		ilr290, ilr290, c1ao->sac[ilr290 - 1][ilr290 - 1].real(), c1ao->sac[ilr290 - 1][ilr290 - 1].imag(),
		jlr, ilr290, c1ao->sac[jlr - 1][ilr290 - 1].real(), c1ao->sac[jlr - 1][ilr290 - 1].imag()
		);
		fprintf(
		output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
		ilr290, ilr290, refinr, ilr290, ilr290, extcor
		);
		fprintf(
		output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
		qschu, pschu, s0mag
		);
		bool goto190 = isam >= 0 && (jths > 1 \|\| jphs > 1);
		if (!goto190) {
		gapv[0] = gap[0][ilr290 - 1];
		gapv[1] = gap[1][ilr290 - 1];
		gapv[2] = gap[2][ilr290 - 1];
		double extins = c1ao->ecsc[ilr290 - 1];
		double scatts = c1ao->scsc[ilr290 - 1];
		double rapr, cosav, fp, fn, fk, fx, fy, fz;
		rftr(u, up, un, gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz);
		fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
		fprintf(output, " Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n", fp, fn, fk);
		fprintf(output, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n", fx, fy, fz);
		tqev[0] = tqce[ilr290 - 1][0];
		tqev[1] = tqce[ilr290 - 1][1];
		tqev[2] = tqce[ilr290 - 1][2];
		tqsv[0] = tqcs[ilr290 - 1][0];
		tqsv[1] = tqcs[ilr290 - 1][1];
		tqsv[2] = tqcs[ilr290 - 1][2];
		double tep, ten, tek, tsp, tsn, tsk;
		tqr(u, up, un, tqev, tqsv, tep, ten, tek, tsp, tsn, tsk);
		fprintf(output, " TQEl=%15.7lE, TQEr=%15.7lE, TQEk=%15.7lE\n", tep, ten, tek);
		fprintf(output, " TQSl=%15.7lE, TQSr=%15.7lE, TQSk=%15.7lE\n", tsp, tsn, tsk);
		fprintf(
		output, " TQEx=%15.7lE, TQEy=%15.7lE, TQEz=%15.7lE\n",
		tqce[ilr290 - 1][0], tqce[ilr290 - 1][1], tqce[ilr290 - 1][2]
		);
		fprintf(
		output, " TQSx=%15.7lE, TQSy=%15.7lE, TQSz=%15.7lE\n",
		tqcs[ilr290 - 1][0], tqcs[ilr290 - 1][1], tqcs[ilr290 - 1][2]
		);
		}
		} //ilr290 loop
		double rbirif = (c1ao->fsac[0][0].real() - c1ao->fsac[1][1].real()) / c1ao->fsac[0][0].real();
		double rdichr = (c1ao->fsac[0][0].imag() - c1ao->fsac[1][1].imag()) / c1ao->fsac[0][0].imag();
		fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rbirif);
		fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rdichr);
		fprintf(output, " MULC\n");
		for (int i = 0; i < 4; i++) {
		fprintf(
		output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
		cmul[i][0], cmul[i][1], cmul[i][2], cmul[i][3]
		);
		}
		fprintf(output, " MULCLR\n");
		for (int i = 0; i < 4; i++) {
		fprintf(
		output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
		cmullr[i][0], cmullr[i][1], cmullr[i][2], cmullr[i][3]
		);
		}
		if (iavm != 0) {
		mmulc(c1ao->vintm, cmullr, cmul);
		// Some implicit loops writing to binary.
		for (int i = 0; i < 16; i++) {
		double value;
		value = c1ao->vintm[i].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->vintm[i].imag();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		for (int i = 0; i < 4; i++) {
		for (int j = 0; j < 4; j++) {
		double value = cmul[i][j];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		}
		fprintf(output, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
		if (inpol == 0) {
		fprintf(output, " LIN\n");
		} else { // label 316
		fprintf(output, " CIRC\n");
		}
		// label 318
		fprintf(output, " MULC\n");
		for (int i = 0; i < 4; i++) {
		fprintf(
		output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
		cmul[i][0], cmul[i][1], cmul[i][2], cmul[i][3]
		);
		}
		fprintf(output, " MULCLR\n");
		for (int i = 0; i < 4; i++) {
		fprintf(
		output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
		cmullr[i][0], cmullr[i][1], cmullr[i][2], cmullr[i][3]
		);
		}
		}
		// label 420, continues jphs loop
		if (isam < 1) phs += phsstp;
		} // jphs loop, labeled 480
		if (isam <= 1) thsl += thsstp;
		} // jths loop, labeled 482
		ph += phstp;
		} // jph484 loop
		th += thstp;
		} // jth486 loop
		printf("INFO: done jxi488 iteration.\n");
		} // jxi488 loop
		@@ -609,6 +916,8 @@ void cluster() {
		delete[] tqcpe;
		delete[] tqcs;
		delete[] tqcps;
		delete[] tqev;
		delete[] tqsv;
		for (int gi = 2; gi > -1; gi--) {
		delete[] gapp[gi];
		delete[] gappm[gi];
		@@ -619,6 +928,7 @@ void cluster() {
		delete[] gappm;
		delete[] gap;
		delete[] gapm;
		delete[] gapv;
		delete[] u;
		delete[] us;
		delete[] un;
		@@ -635,9 +945,13 @@ void cluster() {
		for (int ci = 3; ci > -1; ci--) {
		delete[] cextlr[ci];
		delete[] cext[ci];
		delete[] cmullr[ci];
		delete[] cmul[ci];
		}
		delete[] cextlr;
		delete[] cext;
		delete[] cmullr;
		delete[] cmul;
		} else { // NSPH mismatch between geometry and scatterer configurations.
		throw UnrecognizedConfigurationException(
		"Inconsistent geometry and scatterer configurations."

src/include/Commons.h

+4 −0

Original line number	Diff line number	Diff line
		@@ -238,8 +238,12 @@ public:
		//! \brief QUESTION: definition?
		std::complex<double> **fsacm;
		//! \brief QUESTION: definition?
		double *scsc;
		//! \brief QUESTION: definition?
		std::complex<double> *scscp;
		//! \brief QUESTION: definition?
		double *ecsc;
		//! \brief QUESTION: definition?
		double *ecscm;
		//! \brief QUESTION: definition?
		double *scscm;

src/include/clu_subs.h

+221 −18

Original line number	Diff line number	Diff line
		@@ -34,6 +34,7 @@ extern void rabas(
		);
		extern void rbf(int n, double x, int &nm, double sj[]);
		extern void rnf(int n, double x, int &nm, double sy[]);
		extern void mmulc(std::complex<double> vint, double cmullr, double *cmul);
		extern void sphar(double cth, double sth, double cph, double sph, int lm, std::complex<double> *ylm);
		extern void thdps(int lm, double ****zpv);
		extern void upvmp(
		@@ -57,7 +58,6 @@ extern void wmasp(
		int nsph, double argi, double args, C1 *c1
		);
		// >>> END OF SPH_SUBS DECLARATION <<<
		void logmat(std::complex<double> **mat, int rows, int cols);

		/*! \brief C++ porting of CDTP
		*
		@@ -1415,6 +1415,74 @@ void mextc(double vk, double exri, std::complex<double> fsac, double cextlr,
		}
		}

		/*! \brief C++ porting of PCROS
		*
		* This function is intended to evaluate the particle cross-section. QUESTIUON: correct?
		* \param vk: `double`
		* \param exri: `double`
		* \param c1: `C1 *`
		* \param c1ao: `C1_AddOns *`
		* \param c4: `C4 *`
		*/
		void pcros(double vk, double exri, C1 c1, C1_AddOns c1ao, C4 *c4) {
		const std::complex<double> cc0(0.0, 0.0);
		std::complex<double> sump, sum1, sum2, sum3, sum4, am, amp, cc, csam;
		const double exdc = exri * exri;
		double ccs = 1.0 / (vk * vk);
		double cccs = ccs / exdc;
		csam = -(ccs / (exri * vk)) * std::complex<double>(0.0, 0.5);
		const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
		double cfsq =4.0 / (pi4sq ccs ccs);
		const int nlemt = c4->nlem + c4->nlem;
		int jpo = 2;
		for (int ipo18 = 1; ipo18 <= 2; ipo18++) {
		if (ipo18 == 2) jpo = 1;
		int ipopt = ipo18 + 2;
		int jpopt = jpo + 2;
		double sum = 0.0;
		sump = cc0;
		sum1 = cc0;
		sum2 = cc0;
		sum3 = cc0;
		sum4 = cc0;
		for (int i12 = 1; i12 <= nlemt; i12++) {
		int i = i12 - 1;
		am = cc0;
		amp = cc0;
		for (int j10 = 1; j10 <= nlemt; j10++) {
		int j = j10 - 1;
		am += (c1ao->am0m[i][j] * c1->w[j][ipo18 - 1]);
		amp += (c1ao->am0m[i][j] * c1->w[j][jpo - 1]);
		} // j10 loop
		sum += (dconjg(am) * am).real();
		sump += (dconjg(amp) * amp);
		sum1 += (dconjg(c1->w[i][ipo18 - 1]) * am);
		sum2 += (dconjg(c1->w[i][jpo - 1]) * am);
		sum3 += (c1->w[i][ipopt - 1] * am);
		sum4 += (c1->w[i][jpopt - 1] * am);
		} // i12 loop
		c1ao->scsc[ipo18 - 1] = cccs * sum;
		c1ao->scscp[ipo18 - 1] = cccs * sump;
		c1ao->ecsc[ipo18 - 1] = -cccs * sum1.real();
		c1ao->ecscp[ipo18 - 1] = -cccs * sum2;
		c1ao->fsac[ipo18 - 1][ipo18 - 1] = csam * sum1;
		c1ao->fsac[jpo - 1][ipo18 - 1] = csam * sum2;
		c1ao->sac[ipo18 - 1][ipo18 - 1] = csam * sum3;
		c1ao->sac[jpo - 1][ipo18 - 1] = csam * sum4;
		} // ipo18 loop
		int i = 0;
		for (int ipo1 = 1; ipo1 <= 2; ipo1++) {
		for (int jpo1 = 1; jpo1 <= 2; jpo1++) {
		cc = dconjg(c1ao->sac[jpo1 - 1][ipo1 - 1]);
		for (int ipo2 = 1; ipo2 <= 2; ipo2 ++) {
		for (int jpo2 = 1; jpo2 <= 2; jpo2++) {
		c1ao->vint[i++] = c1ao->sac[jpo2 - 1][ipo2 - 1] * cc * cfsq;
		} // jpo2 loop
		} // ipo2 loop
		} // jpo1 loop
		} // ipo1 loop
		}

		/*! \brief C++ porting of PCRSM0
		*
		* \param vk: `double`
		@@ -1800,6 +1868,39 @@ void raba(
		delete[] ctqcs;
		}

		/*! \brief C++ porting of RFTR
		*
		* \param u: `double *`
		* \param up: `double *`
		* \param un: `double *`
		* \param gapv: `double *`
		* \param extins: `double`
		* \param scatts: `double`
		* \param rapr: `double &`
		* \param cosav: `double &`
		* \param fp: `double &`
		* \param fn: `double &`
		* \param fk: `double &`
		* \param fx: `double &`
		* \param fy: `double &`
		* \param fz: `double &`
		*/
		void rftr(
		double u, double up, double un, double gapv, double extins, double scatts,
		double &rapr, double &cosav, double &fp, double &fn, double &fk, double &fx,
		double &fy, double &fz
		) {
		fk = u[0] * gapv[0] + u[1] * gapv[1] + u[2] * gapv[2];
		rapr = extins - fk;
		cosav = fk / scatts;
		fp = -(up[0] * gapv[0] + up[1] * gapv[1] + up[2] * gapv[2]);
		fn = -(un[0] * gapv[0] + un[1] * gapv[1] + un[2] * gapv[2]);
		fk = rapr;
		fx = u[0] * extins - gapv[0];
		fy = u[1] * extins - gapv[1];
		fz = u[2] * extins - gapv[2];
		}

		/*! \brief C++ porting of SCR0
		*
		* \param vk: `double` QUESTION: definition?
		@@ -1852,6 +1953,99 @@ void scr0(double vk, double exri, C1 c1, C1_AddOns c1ao, C3 c3, C4 c4) {
		} // i14 loop
		}

		/*! \brief C++ porting of SCR2
		*
		* \param vk: `double` QUESTION: definition?
		* \param vkarg: `double` QUESTION: definition?
		* \param exri: `double` External medium refractive index
		* \param duk: `double *` QUESTION: definition?
		* \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 scr2(
		double vk, double vkarg, double exri, double duk, C1 c1, C1_AddOns *c1ao,
		C3 c3, C4 c4) {
		const std::complex<double> cc0(0.0, 0.0);
		const std::complex<double> uim(0.0, 1.0);
		std::complex<double> s11, s21, s12, s22, rm, re, csam, cph, phas, cc;
		double ccs = 1.0 / (vk * vk);
		csam = -(ccs / (exri * vk)) * std::complex<double>(0.0, 0.5);
		const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
		double cfsq = 4.0 / (pi4sq * ccs * ccs);
		cph = uim * exri * vkarg;
		int ls = (c4->li < c4->le) ? c4->li : c4->le;
		c3->tsas[0][0] = cc0;
		c3->tsas[1][0] = cc0;
		c3->tsas[0][1] = cc0;
		c3->tsas[1][1] = cc0;
		for (int i14 = 1; i14 <= c4->nsph; i14++) {
		int i = i14 - 1;
		int iogi = c1->iog[i14 - 1];
		if (iogi >= i14) {
		int k = 0;
		s11 = cc0;
		s21 = cc0;
		s12 = cc0;
		s22 = cc0;
		for (int l10 = 1; l10 <= ls; l10++) {
		int l = l10 - 1;
		rm = 1.0 / c1->rmi[l][i];
		re = 1.0 / c1->rei[l][i];
		int ltpo = l10 + l10 + 1;
		for (int im10 = 1; im10 <= ltpo; im10++) {
		k++;
		int ke = k + c4->nlem;
		s11 -= (c1->w[k - 1][2] * c1->w[k - 1][0] * rm + c1->w[ke - 1][2] * c1->w[ke - 1][0] * re);
		s21 -= (c1->w[k - 1][3] * c1->w[k - 1][0] * rm + c1->w[ke - 1][3] * c1->w[ke - 1][0] * re);
		s12 -= (c1->w[k - 1][2] * c1->w[k - 1][1] * rm + c1->w[ke - 1][2] * c1->w[ke - 1][1] * re);
		s22 -= (c1->w[k - 1][3] * c1->w[k - 1][1] * rm + c1->w[ke - 1][3] * c1->w[ke - 1][1] * re);
		} // im10 loop
		} // l10 loop
		c1->sas[i][0][0] = s11 * csam;
		c1->sas[i][1][0] = s21 * csam;
		c1->sas[i][0][1] = s12 * csam;
		c1->sas[i][1][1] = s22 * csam;
		}
		// label 12
		phas = exp(cph * (duk[0] * c1->rxx[i] + duk[1] * c1->ryy[i] + duk[2] * c1->rzz[i]));
		c3->tsas[0][0] += (c1->sas[iogi - 1][0][0] * phas);
		c3->tsas[1][0] += (c1->sas[iogi - 1][1][0] * phas);
		c3->tsas[0][1] += (c1->sas[iogi - 1][0][1] * phas);
		c3->tsas[1][1] += (c1->sas[iogi - 1][1][1] * phas);
		} // i14 loop
		for (int i24 = 1; i24 <= c4->nsph; i24++) {
		int iogi = c1->iog[i24 - 1];
		if (iogi >= i24) {
		int j = 0;
		for (int ipo1 = 1; ipo1 <=2; ipo1++) {
		for (int jpo1 = 1; jpo1 <= 2; jpo1++) {
		cc = dconjg(c1->sas[i24 - 1][jpo1 - 1][ipo1 - 1]);
		for (int ipo2 = 1; ipo2 <= 2; ipo2++) {
		for (int jpo2 = 1; jpo2 <= 2; jpo2++) {
		j++;
		c1ao->vints[i24 - 1][j - 1] = c1->sas[i24 - 1][jpo2 - 1][ipo2 - 1] * cc * cfsq;
		} // jpo2 loop
		} // ipo2 loop
		} // jpo1 loop
		} // ipo1 loop
		}
		} // i24 loop
		int j = 0;
		for (int ipo1 = 1; ipo1 <=2; ipo1++) {
		for (int jpo1 = 1; jpo1 <= 2; jpo1++) {
		cc = dconjg(c3->tsas[jpo1 - 1][ipo1 - 1]);
		for (int ipo2 = 1; ipo2 <= 2; ipo2++) {
		for (int jpo2 = 1; jpo2 <= 2; jpo2++) {
		j++;
		c1ao->vintt[j - 1] = c3->tsas[jpo2 - 1][ipo2 - 1] * cc * cfsq;
		} // jpo2 loop
		} // ipo2 loop
		} // jpo1 loop
		} // ipo1 loop
		}

		/*! \brief C++ porting of STR
		*
		* \param rcf: `double **` Matrix of double.
		@@ -1936,6 +2130,32 @@ void str(double *rcf, C1 c1, C1_AddOns c1ao, C3 c3, C4 c4, C6 c6) {
		delete[] ylm;
		}

		/*! \brief C++ porting of TQR
		*
		* \param u: `double *`
		* \param up: `double *`
		* \param un: `double *`
		* \param tqev: `double *`
		* \param tqsv: `double *`
		* \param tep: `double &`
		* \param ten: `double &`
		* \param tek: `double &`
		* \param tsp: `double &`
		* \param tsn: `double &`
		* \param tsk: `double &`
		*/
		void tqr(
		double u, double up, double un, double tqev, double *tqsv, double &tep,
		double &ten, double &tek, double &tsp, double &tsn, double &tsk
		) {
		tep = up[0] * tqev[0] + up[1] * tqev[1] + up[2] * tqev[2];
		ten = un[0] * tqev[0] + un[1] * tqev[1] + un[2] * tqev[2];
		tek = u[0] * tqev[0] + u[1] * tqev[1] + u[2] * tqev[2];
		tsp = up[0] * tqsv[0] + up[1] * tqsv[1] + up[2] * tqsv[2];
		tsn = un[0] * tqsv[0] + un[1] * tqsv[1] + un[2] * tqsv[2];
		tsk = u[0] * tqsv[0] + u[1] * tqsv[1] + u[2] * tqsv[2];
		}

		/*! \brief C++ porting of ZTM
		*
		* \param am: Matrix of complex.
		@@ -2025,23 +2245,6 @@ void ztm(std::complex<double> *am, C1 c1, C1_AddOns c1ao, C4 c4, C6 *c6, C9
		} // i0 loop
		}

		/*! \brief Write a matrix to a log file (debug function).
		*
		* \param mat: Matrix of complex.
		* \param rows: `int`
		* \param cols: `int`
		*/
		void logmat(std::complex<double> **mat, int rows, int cols) {
		FILE* logfile = fopen("c_matrix.log", "w");
		for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
		fprintf(logfile, "R:%3d,%3d,%15.7lE\n", i + 1, j + 1, mat[i][j].real());
		fprintf(logfile, "I:%3d,%3d,%15.7lE\n", i + 1, j + 1, mat[i][j].imag());
		}
		}
		fclose(logfile);
		}

		/*! \brief Sum all the elements of a matrix (debug function).
		*
		* \param mat: Matrix of complex.

src/libnptm/Commons.cpp

+4 −0

File changed.

Preview size limit exceeded, changes collapsed.