Improve format matching with legacy output (8d84efcd) · Commits · Giacomo Mulas / NP_TMcode

src/cluster/cluster.cpp

+178 −6

Original line number	Diff line number	Diff line
		@@ -136,6 +136,18 @@ void cluster() {
		double *unsmp = new double[3]();
		double *upmp = new double[3]();
		double *upsmp = new double[3]();
		double *argi = new double[1]();
		double *args = new double[1]();
		double *duk = new double[3]();
		double cextlr, cext;
		cextlr = new double*[4];
		cext = new double*[4];
		for (int ci = 0; ci < 4; ci++) {
		cextlr[ci] = new double[4]();
		cext[ci] = new double[4]();
		}
		int isq, ibf;
		double scan, cfmp, sfmp, cfsp, sfsp;
		// 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",
		@@ -338,7 +350,7 @@ void cluster() {
		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());
		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");
		@@ -373,14 +385,11 @@ void cluster() {
		for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP portable?
		ph = ph1;
		double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
		// argi[NSPEF], with NSPEF=1 if IDOT=0, else NSPEF=NSPH
		double *argi;
		for (int jph484 = 1; jph484 <= nph; jph484++) {
		int jw = 0;
		if (nk != 1 \|\| jxi488 <= 1) {
		upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
		if (isam >= 0) {
		argi = new double[1];
		wmamp(
		0, cost, sint, cosp, sinp, inpol, c4->le, 0,
		nsph, argi, u, upmp, unmp, c1
		@@ -400,6 +409,160 @@ void cluster() {
		}
		// label 184
		double thsl = ths1;
		double phsph = 0.0;
		for (int jths = 1; jths <= nths; jths++) {
		ths = thsl;
		int icspnv = 0;
		if (isam > 1) ths += thsca;
		if (isam >= 1) {
		phsph = 0.0;
		if (ths < 0.0 \|\| ths > 180.0) phsph = 180.0;
		if (ths < 0.0) ths *= -1.0;
		if (ths > 180.0) ths = 360.0 - ths;
		if (phsph != 0.0) icspnv = 1;
		}
		// label 186
		phs = phs1;
		for (int jphs = 1; jphs <= nphs; jphs++) {
		double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0;
		if (isam >= 1) {
		phs = ph + phsph;
		if (phs > 360.0) phs -= 360.0;
		}
		// label 188
		if (!((nks == 1 && jxi488 == 1) \|\| jth486 > 1 \|\| jph484 > 1)) {
		upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
		if (isam >= 0)
		wmamp(
		2, costs, sints, cosps, sinps, inpol, c4->le,
		0, nsph, args, us, upsmp, unsmp, c1
		);
		}
		// label 190
		if (nkks != 1 \|\| jxi488 <= 1) {
		upvsp(
		u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns,
		duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp
		);
		if (isam < 0) {
		wmasp(
		cost, sint, cosp, sinp, costs, sints, cosps, sinps,
		u, up, un, us, ups, uns, isq, ibf, inpol, c4->le,
		0, nsph, argi, args, c1
		);
		} else { // label 192
		for (int i193 = 0; i193 < 3; i193++) {
		up[i193] = upmp[i193];
		un[i193] = unmp[i193];
		ups[i193] = upsmp[i193];
		uns[i193] = unsmp[i193];
		}
		}
		}
		// label 194
		if (iavm == 1) crsm1(vk, exri, c1, c1ao, c4, c6);
		if (isam < 0) {
		apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
		raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
		jw = 1;
		}
		// label 196
		tppoan.write(reinterpret_cast<char *>(&th), sizeof(double));
		tppoan.write(reinterpret_cast<char *>(&ph), sizeof(double));
		tppoan.write(reinterpret_cast<char *>(&ths), sizeof(double));
		tppoan.write(reinterpret_cast<char *>(&phs), sizeof(double));
		tppoan.write(reinterpret_cast<char *>(&scan), sizeof(double));
		if (jaw != 0) {
		jaw = 0;
		mextc(vk, exri, c1ao->fsacm, cextlr, cext);
		// We now have 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 < 3; i++) {
		double value = c1ao->scscm[i];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->scscpm[i].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->scscpm[i].imag();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->ecscm[i];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->ecscpm[i].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = c1ao->ecscpm[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 = gapm[i][j];
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = gappm[i][j].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = gappm[i][j].imag();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		}
		fprintf(output, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
		int jlr = 2;
		for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
		int ipol = (ilr210 % 2 == 0) ? 1 : -1;
		if (ilr210 == 2) jlr = 1;
		double extsm = c1ao->ecscm[ilr210 - 1];
		double qextm = extsm * sqsfi / c3->gcs;
		double extrm = extsm / c3->ecs;
		double scasm = c1ao->scscm[ilr210 - 1];
		double albdm = scasm / extsm;
		double qscam = scasm *sqsfi / c3->gcs;
		double scarm = scasm / c3->scs;
		double abssm = extsm - scasm;
		double qabsm = abssm * sqsfi / c3->gcs;
		double absrm = abssm / c3->acs;
		double acsecs = c3->acs / c3->ecs;
		if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;
		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 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) {
		fprintf(output, " LIN %2d\n", ipol);
		} else { // label 206
		fprintf(output, " CIRC %2d\n", ipol);
		}
		// label 208
		fprintf(output, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");
		fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", scasm, abssm, extsm, albdm);
		fprintf(output, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");
		fprintf(output, " %14.7lE%15.7lE%15.7lE\n", qscam, qabsm, qextm);
		fprintf(output, " ---- SCCRT --- ABCRT --- EXCRT ----\n");
		fprintf(output, " %14.7lE%15.7lE%15.7lE\n", scarm, absrm, extrm);
		fprintf(
		output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
		ilr210, ilr210, c1ao->fsacm[ilr210 - 1][ilr210 - 1].real(),
		c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag(), jlr, ilr210,
		c1ao->fsacm[jlr - 1][ilr210 - 1].real(), c1ao->fsacm[jlr - 1][ilr210 - 1].imag()
		);
		fprintf(
		output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
		ilr210, ilr210, rfinrm, ilr210, ilr210, extcrm
		);
		fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschum, pschum, s0magm);
		double rapr = c1ao->ecscm[ilr210 - 1] - gapm[2][ilr210 - 1];
		double cosav = gapm[2][ilr210 - 1] / c1ao->scscm[ilr210 - 1];
		double fz = rapr;
		fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
		fprintf(output, " Fk=%15.7lE\n", fz);
		} // ilr210 loop
		// from RMBRIF
		}
		// label 212
		} // jphs loop
		} // jths loop
		} // jph484 loop
		} // jth486 loop
		printf("INFO: done jxi488 iteration.\n");
		@@ -466,6 +629,15 @@ void cluster() {
		delete[] unsmp;
		delete[] upmp;
		delete[] upsmp;
		delete[] argi;
		delete[] args;
		delete[] duk;
		for (int ci = 3; ci > -1; ci--) {
		delete[] cextlr[ci];
		delete[] cext[ci];
		}
		delete[] cextlr;
		delete[] cext;
		} else { // NSPH mismatch between geometry and scatterer configurations.
		throw UnrecognizedConfigurationException(
		"Inconsistent geometry and scatterer configurations."

src/include/Commons.h

+2 −0

Original line number	Diff line number	Diff line
		@@ -228,6 +228,8 @@ public:
		//! \brief QUESTION: definition?
		std::complex<double> *vintm;
		//! \brief QUESTION: definition?
		std::complex<double> **vints;
		//! \brief QUESTION: definition?
		std::complex<double> *vintt;
		//! \brief QUESTION: definition?
		std::complex<double> **fsac;

src/include/clu_subs.h

+309 −20

Original line number	Diff line number	Diff line
		@@ -40,11 +40,22 @@ extern void upvmp(
		double thd, double phd, int icspnv, double &cost, double &sint,
		double &cosp, double &sinp, double u, double up, double *un
		);
		extern void upvsp(
		double u, double upmp, double unmp, double us, double upsmp, double unsmp,
		double up, double un, double ups, double uns, double *duk, int &isq,
		int &ibf, double &scand, double &cfmp, double &sfmp, double &cfsp, double &sfsp
		);
		extern void wmamp(
		int iis, double cost, double sint, double cosp, double sinp, int inpol,
		int lm, int idot, int nsph, double arg, double u, double *up,
		double un, C1 c1
		);
		extern void wmasp(
		double cost, double sint, double cosp, double sinp, double costs, double sints,
		double cosps, double sinps, double u, double up, double un, double us,
		double ups, double uns, int isq, int ibf, int inpol, int lm, int idot,
		int nsph, double argi, double args, C1 *c1
		);
		// >>> END OF SPH_SUBS DECLARATION <<<
		void logmat(std::complex<double> **mat, int rows, int cols);

		@@ -245,6 +256,109 @@ void r3jjr(int j2, int j3, int m2, int m3, C6 *c6) {
		}
		}

		/*! \brief C++ porting of R3JMR
		*
		* \param j1: `int`
		* \param j2: `int`
		* \param j3: `int`
		* \param m1: `int`
		* \param c6: `C6 *`
		*/
		void r3jmr(int j1, int j2, int j3, int m1, C6 *c6) {
		int mmx = (j2 < j3 - m1) ? j2 : j3 - m1;
		int mmn = (-j2 > -(j3 + m1)) ? -j2 : -(j3 + m1);
		int nmmo = mmx - mmn;
		int j1po = j1 + 1;
		int j1tpo = j1po + j1;
		int isn = 1;
		if ((j2 - j3 - m1) % 2 != 0) isn = -1;
		if (nmmo <= 0) {
		double sj = 1.0 * j1tpo;
		double cnr = (1.0 / sqrt(sj)) * isn;
		c6->rac3j[0] = cnr;
		// returns
		} else { // label 15
		int j1s = j1 * j1po;
		int j2po = j2 + 1;
		int j2s = j2 * j2po;
		int j3po = j3 + 1;
		int j3s = j3 * j3po;
		int id = j1s - j2s - j3s;
		int m2 = mmx;
		int m3 = m1 + m2;
		double cm = sqrt(1.0 * (j2po - m2) * (j2 + m2) * (j3po - m3) * (j3 + m3));
		double dm = 1.0 * (id + m2 * m3 * 2);
		if (nmmo <= 1) {
		c6->rac3j[0] = dm / cm;
		double sj = (1.0 + c6->rac3j[0] * c6->rac3j[0]) * j1tpo;
		double cnr = 1.0 / sqrt(sj) * isn;
		c6->rac3j[1] = cnr;
		c6->rac3j[0] *= cnr;
		// returns
		} else { // label 20
		int nm = nmmo + 1;
		int nmat = (nm + 1) / 2;
		c6->rac3j[nm - 1] = 1.0;
		c6->rac3j[nmmo - 1] = dm / cm;
		double sjt = 1.0;
		double sjr = 1.0;
		if (nmat != nmmo) {
		int nbr = nmmo - nmat;
		for (int ibr45 = 1; ibr45 <= nbr; ibr45++) {
		int irr = nm - ibr45;
		m2--;
		m3 = m1 + m2;
		double cmp = cm;
		cm = sqrt(1.0 * (j2po - m2) * (j2 + m2) * (j3po - m3) * (j3 + m3));
		sjt = c6->rac3j[irr - 1] * c6->rac3j[irr - 1];
		dm = 1.0 * (id + m2 * m3 * 2);
		c6->rac3j[irr - 1] = ((dm - c6->rac3j[irr] cmp) / cm);
		sjr += sjt;
		} // ibr45 loop
		}
		// label 50
		double osjt = sjt;
		sjt = c6->rac3j[nmat - 1] * c6->rac3j[nmat - 1];
		if (sjt >= osjt) {
		sjr += sjt;
		} else { // label 55
		nmat++;
		}
		// label 60
		double racmat = c6->rac3j[nmat - 1];
		c6->rac3j[0] = 1.0;
		m2 = mmn;
		m3 = m1 + m2;
		double cmp = sqrt(1.0 * (j2 - m2) * (j2po + m2) * (j3 - m3) * (j3po + m3));
		dm = 1.0 * (id + m2 * m3 * 2);
		c6->rac3j[1] = dm / cmp;
		double sjl = 1.0;
		int nmatmo = nmat - 1;
		if (nmatmo > 1) {
		for (int irl70 = 2; irl70 <= nmatmo; irl70++) {
		m2++;
		m3 = m1 + m2;
		cm = cmp;
		cmp = sqrt(1.0 * (j2 - m2) * (j2po + m2) * (j3 - m3) * (j3po + m3));
		sjt = c6->rac3j[irl70 - 1] * c6->rac3j[irl70 - 1];
		dm = 1.0 * (id + m2 * m3 * 2);
		c6->rac3j[irl70] = (c6->rac3j[irl70 - 1] * dm - c6->rac3j[irl70 - 2] * cm) / cmp;
		sjl += sjt;
		}
		}// label 75
		double ratrac = racmat / c6->rac3j[nmat - 1];
		double rats = ratrac * ratrac;
		double sj = (sjr + sjl * rats) * j1tpo;
		c6->rac3j[nmat - 1] = racmat;
		double cnr = 1.0 / sqrt(sj) * isn;
		for (int irr80 = nmat; irr80 <= nm; irr80++) c6->rac3j[irr80 - 1] *= cnr;
		double cnl = cnr * ratrac;
		for (int irl85 = 1; irl85 <= nmatmo; irl85++) c6->rac3j[irl85 - 1] *= cnl;
		// returns
		}
		}
		}

		/*! \brief C++ porting of GHIT
		*
		* \param ihi: `int`
		@@ -663,7 +777,7 @@ void apcra(
		// label 26
		} // ilmp loop
		} // l28 loop
		for (int l30 = 1; l30 <= le; l30++) { // This can be omitted, since svs was initialized at 0
		for (int l30 = 1; l30 <= le; l30++) { // 0-init: can be omitted
		for (int ilmp = 1; ilmp <= 3; ilmp++) {
		for (int ipa = 1; ipa <= 2; ipa++) {
		for (int ipamp = 1; ipamp <= 2; ipamp++) {
		@@ -757,9 +871,9 @@ void apcra(
		sum1 = cc0;
		sum2 = cc0;
		for (int l68 = 1; l68 <= le; l68++) {
		int lpo = l68 + 1;
		int ltpo = l68 + lpo;
		int imm = l68 * lpo;
		//int lpo = l68 + 1;
		//int ltpo = l68 + lpo;
		//int imm = l68 * lpo;
		for (int ilmp = 1; ilmp <= 3; ilmp++) {
		if ((l68 == 1 && ilmp == 1) \|\| (l68 == le && ilmp == 3)) continue; // ilmp loop
		if (inpol == 0) {
		@@ -921,6 +1035,135 @@ void cms(std::complex<double> *am, C1 c1, C1_AddOns c1ao, C4 c4, C6 *c6) {
		} // n1 loop
		}

		/*! \brief C++ porting of CRSM1
		*
		* \param vk: `double`
		* \param exri: `double`
		* \param c1: `C1 *`
		* \param c1ao: `C1_AddOns *`
		* \param c4: `C4 *`
		* \param c6: `C6 *`
		*/
		void crsm1(double vk, double exri, C1 c1, C1_AddOns c1ao, C4 c4, C6 c6) {
		std::complex<double> *svf, svw, *svs;
		const std::complex<double> cc0(0.0, 0.0);
		std::complex<double> cam(0.0, 0.0);
		const int le4po = 4 * c4->le + 1;
		svf = new std::complex<double>**[le4po];
		svw = new std::complex<double>**[le4po];
		svs = new std::complex<double>*[le4po];
		for (int si = 0; si < le4po; si++) {
		svf[si] = new std::complex<double>*[le4po];
		svw[si] = new std::complex<double>*[4];
		svs[si] = new std::complex<double>[4]();
		for (int sj = 0; sj < le4po; sj++) svf[si][sj] = new std::complex<double>[4]();
		for (int sj = 0; sj < 4; sj++) svw[si][sj] = new std::complex<double>[4]();
		}
		double exdc = exri * exri;
		double ccs = 1.0 / (vk * vk);
		const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
		double cint = ccs / (pi4sq * exdc);
		int letpo = c4->le + c4->le + 1;
		for (int i20 = 0; i20 < 16; i20++) c1ao->vintm[i20] = cc0; // 0-init: can be omitted
		for (int lpo40 = 1; lpo40 <= letpo; lpo40++) {
		int l = lpo40 - 1;
		int ltpo = lpo40 + l;
		int immn = letpo - l;
		int immx = letpo + l;
		for (int imf = immn; imf <= immx; imf++) { // 0-init: can be omitted
		for (int ims = immn; ims <= immx; ims++) {
		for (int ipo = 1; ipo <= 4; ipo++) {
		svf[imf - 1][ims - 1][ipo - 1] = cc0;
		} // ipo loop
		} // ims loop
		} // imf loop
		for (int l1 = 1; l1 <= c4->le; l1++) {
		int il1 = l1 * (l1 + 1);
		for (int l2 = 1; l2 <= c4->le; l2++) {
		int abs_l2ml1 = (l2 > l1) ? l2 - l1 : l1 - l2;
		if (l < abs_l2ml1 \|\| l > l2 + l1) continue; // l2 loop
		int il2 = l2 * (l2 + 1);
		for (int im = immn; im >= immx; im++) { // 0-init: can be omitted
		for (int ipa = 1; ipa <= 4; ipa++) {
		svs[im - 1][ipa - 1] = cc0;
		for (int ipo = 1; ipo <= 4; ipo++) {
		svw[im - 1][ipa - 1][ipo - 1] = cc0;
		} // ipo loop
		} // ipa loop
		} // im loop
		for (int im = immn; im <= immx; im++) {
		int m = im - letpo;
		r3jmr(l, l1, l2, m, c6);
		int m1mnmo = (-l1 > -l2 - m) ? -(l1 + 1) : -(l2 + m + 1);
		int nm1 = (l1 < l2 - m) ? (l1 - m1mnmo) : (l2 - m - m1mnmo);
		for (int im1 = 1; im1 <= nm1; im1++) {
		int m1 = -im1 - m1mnmo;
		int isn = 1;
		if (m1 % 2 != 0) isn = -1;
		double cg3j = c6->rac3j[im1 - 1] * isn;
		int ilm1 = il1 + m1;
		int ilm2 = il2 + m1 - m;
		int ipa = 0;
		for (int ipa1 = 1; ipa1 <= 2; ipa1++) {
		int i1 = ilm1;
		if (ipa1 == 2) i1 = ilm1 + c4->nlem;
		for (int ipa2 = 1; ipa2 <= 2; ipa2++) {
		int i2 = ilm2;
		if (ipa2 == 2) i2 = ilm2 + c4->nlem;
		ipa++;
		svs[im - 1][ipa - 1] += (c1ao->am0m[i1 - 1][i2 - 1] * cg3j);
		int ipo = 0;
		for (int ipo2 = 1; ipo2 <= 2; ipo2++) {
		for (int ipo1 = 3; ipo1 <= 4; ipo1++) {
		ipo++;
		svw[im - 1][ipa - 1][ipo - 1] += (c1->w[i1 - 1][ipo1 - 1] * c1->w[i2 - 1][ipo2 - 1] * cg3j);
		} // ipo1 loop
		} // ipo2 loop
		} // ipa2 loop
		} // ipa1 loop
		} // im1 loop
		// label 32 loops
		for (int imf = immn; imf <= immx; imf++) {
		for (int ims = immn; ims <= immx; ims++) {
		for (int ipo = 1; ipo <= 4; ipo++) {
		for (int ipa = 1; ipa <= 4; ipa++) {
		svf[imf - 1][ims - 1][ipo - 1] += (svw[imf - 1][ipa - 1][ipo - 1] * svs[ims - 1][ipa - 1]);
		} // ipa loop
		} // ipo loop
		} // ims loop
		} // imf loop
		// ends loop level 34, which are l2 loop and l1 loop
		} // im loop
		} // l2 loop
		} // l1 loop
		for (int imf = immn; imf <= immx; imf++) {
		for (int ims = immn; ims <= immx; ims++) {
		int i = 0;
		for (int ipo1 = 1; ipo1 <= 4; ipo1++) {
		cam = dconjg(svf[imf - 1][ims - 1][ipo1 - 1]);
		for (int ipo2 = 1; ipo2 <= 4; ipo2++) {
		i++;
		c1ao->vintm[i - 1] += (svf[imf - 1][ims - 1][ipo2 - 1] * cam * (1.0 * ltpo));
		}
		} // ipo1 loop
		} // ims loop
		} // imf loop
		} // lpo40 loop
		for (int i42 = 0; i42 < 16; i42++) c1ao->vintm[i42] *= cint;

		// Clean memory
		for (int si = le4po - 1; si > -1; si--) {
		for (int sj = le4po - 1; sj > -1; sj--) delete[] svf[si][sj];
		for (int sj = 3; sj > -1; sj--) delete[] svw[si][sj];
		delete[] svf[si];
		delete[] svw[si];
		delete[] svs[si];
		}
		delete[] svf;
		delete[] svw;
		delete[] svs;
		}

		/*! \brief C++ porting of HJV
		*
		* \param exri: `double`
		@@ -1114,6 +1357,64 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
		delete[] v;
		}

		/*! \brief C++ porting of MEXTC
		*
		* \param vk: `double`
		* \param exri: `double`
		* \param fsac: Matrix of complex
		* \param cextlr: `double **`
		* \param cext: `double **`
		*/
		void mextc(double vk, double exri, std::complex<double> fsac, double cextlr, double **cext) {
		double fa11r = fsac[0][0].real();
		double fa11i = fsac[0][0].imag();
		double fa21r = fsac[1][0].real();
		double fa21i = fsac[1][0].imag();
		double fa12r = fsac[0][1].real();
		double fa12i = fsac[0][1].imag();
		double fa22r = fsac[1][1].real();
		double fa22i = fsac[1][1].imag();
		cextlr[0][0] = fa11i * 2.0;
		cextlr[0][1] = 0.0;
		cextlr[0][2] = -fa12i;
		cextlr[0][3] = -fa12r;
		cextlr[1][0] = 0.0;
		cextlr[1][1] = fa22i * 2.0;
		cextlr[1][2] = -fa21i;
		cextlr[1][3] = fa21r;
		cextlr[2][0] = -fa21i * 2.0;
		cextlr[2][1] = -fa12i * 2.0;
		cextlr[2][2] = fa11i + fa22i;
		cextlr[2][3] = fa22r - fa11r;
		cextlr[3][0] = fa21r * 2.0;
		cextlr[3][1] = -fa12r * 2.0;
		cextlr[3][2] = fa11r - fa22r;
		cextlr[3][3] = cextlr[2][2];
		cext[0][0] = cextlr[3][3];
		cext[1][1] = cextlr[3][3];
		cext[2][2] = cextlr[3][3];
		cext[2][3] = cextlr[2][3];
		cext[3][2] = cextlr[3][2];
		cext[3][3] = cextlr[3][3];
		cext[0][1] = fa11i - fa22i;
		cext[0][2] = -fa12i - fa21i;
		cext[0][3] = fa21r - fa12r;
		cext[1][0] = cext[0][1];
		cext[1][2] = fa21i - fa12i;
		cext[3][1] = fa12r + fa21r;
		cext[1][3] = -cext[3][1];
		cext[2][0] = cext[0][2];
		cext[2][1] = -cext[1][2];
		cext[3][0] = cext[1][3];
		double ckm = vk / exri;
		for (int i10 = 0; i10 < 4; i10++) {
		for (int j10 = 0; j10 < 4; j10++) {
		cextlr[i10][j10] *= ckm;
		cext[i10][j10] *= ckm;
		}
		}
		}

		/*! \brief C++ porting of PCRSM0
		*
		* \param vk: `double`
		@@ -1515,7 +1816,7 @@ void scr0(double vk, double exri, C1 c1, C1_AddOns c1ao, C3 c3, C4 c4) {
		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;
		//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];
		@@ -1544,9 +1845,9 @@ void scr0(double vk, double exri, C1 c1, C1_AddOns c1ao, C3 c3, C4 c4) {
		c1->fsas[i14 - 1] = sum21 * csam;
		}
		// label 12
		scs += c1->sscs[iogi - 1];
		ecs += c1->sexs[iogi - 1];
		acs += c1->sabs[iogi - 1];
		c3->scs += c1->sscs[iogi - 1];
		c3->ecs += c1->sexs[iogi - 1];
		c3->acs += c1->sabs[iogi - 1];
		c3->tfsas += c1->fsas[iogi - 1];
		} // i14 loop
		}
		@@ -1670,18 +1971,6 @@ void ztm(std::complex<double> *am, C1 c1, C1_AddOns c1ao, C4 c4, C6 *c6, C9
		} // im2 loop
		} // l2 loop
		} // n2 loop
		/* // DEBUG CODE
		std::complex<double> dbgtst(0.0, 0.0);
		for (int di = 0; di < ndi; di++) {
		for (int dj = 0; dj < c4->nlem; dj++) dbgtst += c9->gis[di][dj];
		}
		printf("DEBUG: in ZTM init sgis = (%lE, %lE)\n", dbgtst.real(), dbgtst.imag());
		dbgtst = std::complex<double>(0.0, 0.0);
		for (int di = 0; di < ndi; di++) {
		for (int dj = 0; dj < c4->nlem; dj++) dbgtst += c9->gls[di][dj];
		}
		printf("DEBUG: in ZTM init sgls = (%lE, %lE)\n", dbgtst.real(), dbgtst.imag());
		// END DEBUG */
		for (int i1 = 1; i1 <= ndi; i1++) { // GPU portable?
		int i1e = i1 + ndi;
		for (int i3 = 1; i3 <= c4->nlem; i3++) {

src/libnptm/Commons.cpp

+5 −1

Original line number	Diff line number	Diff line
		@@ -106,9 +106,11 @@ C1_AddOns::C1_AddOns(C4 *c4) {
		for (int ai = 0; ai < nlemt; ai++) {
		am0m[ai] = new complex<double>[nlemt]();
		}
		vint = new complex<double>[16]();
		vint = new complex<double>[16](); // QUESTION: is dimension 16 generally fixed?
		vintm = new complex<double>[16]();
		vintt = new complex<double>[16]();
		vints = new complex<double>*[nsph];
		for (int vi = 0; vi < nsph; vi++) vints[vi] = new complex<double>[16]();
		fsac = new complex<double>*[2];
		sac = new complex<double>*[2];
		fsacm = new complex<double>*[2];
		@@ -143,6 +145,8 @@ C1_AddOns::~C1_AddOns() {
		delete[] vint;
		delete[] vintm;
		delete[] vintt;
		for (int vi = nsph - 1; vi > -1; vi--) delete[] vints[vi];
		delete[] vints;
		for (int fi = 1; fi > -1; fi--) {
		delete[] fsac[fi];
		delete[] sac[fi];