Extend sph work-flow migration to C++ (2571ccb6) · Commits · Giacomo Mulas / NP_TMcode

src/include/sph_subs.h

+520 −0

Original line number	Diff line number	Diff line
		@@ -362,6 +362,188 @@ void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
		}
		}

		/*! \brief C++ porting of ORUNVE
		*
		* \param u1: `double *`
		* \param u2: `double *`
		* \param u3: `double *`
		* \param iorth: `int`
		* \param torth: `double`
		*/
		void orunve( double u1, double u2, double *u3, int iorth, double torth) {
		if (iorth <= 0) {
		double cp = u1[0] * u2[0] + u1[1] * u2[1] + u1[2] * u2[2];
		double abs_cp = cp;
		if (abs_cp < 0.0) abs_cp *= -1.0;
		if (iorth == 0 \|\| abs_cp >= torth) {
		double fn = 1.0 / sqrt(1.0 - cp * cp);
		u3[0] = (u1[1] * u2[2] - u1[2] * u2[1]) * fn;
		u3[1] = (u1[2] * u2[0] - u1[0] * u2[2]) * fn;
		u3[2] = (u1[0] * u2[1] - u1[1] * u2[0]) * fn;
		return;
		}
		}
		u3[0] = u1[1] * u2[2] - u1[2] * u2[1];
		u3[1] = u1[2] * u2[0] - u1[0] * u2[2];
		u3[2] = u1[0] * u2[1] - u1[1] * u2[0];
		}

		/*! \brief C++ porting of PWMA
		*
		* \param up: `double *`
		* \param un: `double *`
		* \param ylm: Vector of complex
		* \param inpol: `int`
		* \param lw: `int`
		* \param isq: `int`
		* \param c1: `C1 *`
		*/
		void pwma(
		double up, double un, std::complex<double> *ylm, int inpol, int lw,
		int isq, C1 *c1
		) {
		//std::complex<double> cp1, cm1, cc1, cp2, c2, cc2, uim, cl;
		double four_pi = 8.0 * acos(0.0);
		int is = isq;
		if (isq == -1) is = 0;
		int ispo = is + 1;
		int ispt = is + 2;
		int nlwm = lw * (lw + 2);
		int nlwmt = nlwm + nlwm;
		double sqrtwi = 1.0 / sqrt(2.0);
		std::complex<double> uim(0.0, 1.0);
		std::complex<double> cm1 = 0.5 * std::complex<double>(up[0], up[1]);
		std::complex<double> cp1 = 0.5 * std::complex<double>(up[0], -up[1]);
		double cz1 = up[2];
		std::complex<double> cm2 = 0.5 * std::complex<double>(un[0], un[1]);
		std::complex<double> cp2 = 0.5 * std::complex<double>(un[0], -un[1]);
		double cz2 = un[2];
		for (int l20 = 1; l20 <= lw; l20++) {
		int lf = l20 + 1;
		int lftl = lf * l20;
		double x = 1.0 * lftl;
		std::complex<double> cl = (four_pi / sqrt(x)) * std::pow(uim, l20);
		int mv = l20 + lf;
		int m = -lf;
		for (int mf20 = 1; mf20 <= mv; mf20++) {
		m += 1;
		int k = lftl + m;
		x = 1.0 * (lftl - m * (m + 1));
		double cp = sqrt(x);
		x = 1.0 * (lftl - m * (m -1));
		double cm = sqrt(x);
		double cz = 1.0 * m;
		c1->w[k - 1][ispo - 1] = dconjg(
		cp1 * cp * ylm[k + 1] +
		cm1 * cm * ylm[k - 1] +
		cz1 * cz * ylm[k]
		) * cl;
		c1->w[k - 1][ispt - 1] = dconjg(
		cp2 * cp * ylm[k + 1] +
		cm2 * cm * ylm[k - 1] +
		cz2 * cz * ylm[k]
		) * cl;
		}
		for (int k30 = 1; k30 <= nlwm; k30++) {
		int i = k30 + nlwm;
		c1->w[i - 1][ispo - 1] = uim * c1->w[k30 - 1][ispt - 1];
		c1->w[i - 1][ispt - 1] = -uim * c1->w[k30 - 1][ispo - 1];
		}
		if (inpol != 0) {
		for (int k40 = 1; k40 <= nlwm; k40++) {
		int i = k40 + nlwm;
		std::complex<double> cc1 = sqrtwi * (c1->w[k40 - 1][ispo - 1] + uim * c1->w[k40 - 1][ispt - 1]);
		std::complex<double> cc2 = sqrtwi * (c1->w[k40 - 1][ispo - 1] - uim * c1->w[k40 - 1][ispt - 1]);
		c1->w[k40 - 1][ispo - 1] = cc2;
		c1->w[i - 1][ispo - 1] = -cc2;
		c1->w[k40 - 1][ispt - 1] = cc1;
		c1->w[i - 1][ispt - 1] = cc1;
		}
		}
		if (isq != 0) {
		for (int i50 = 1; i50 <= 2; i50++) {
		int ipt = i50 + 2;
		int ipis = i50 + is;
		for (int k50 = 1; k50 <= nlwmt; k50++)
		c1->w[k50 - 1][ipt - 1] = dconjg(c1->w[k50 - 1][ipis - 1]);
		}
		}
		}
		}

		/*! \brief C++ porting of RABAS
		*
		* \param inpol: `int`
		* \param li: `int`
		* \param nsph: `int`
		* \param c1: `C1 *`
		* \param TQSE: Matrix of complex.
		* \param TQSPE: Matrix of complex.
		* \param TQSS: Matrix of complex.
		* \param TQSPS: Matrix of complex.
		*/
		void rabas(
		int inpol, int li, int nsph, C1 c1, double tqse, std::complex<double> *tqspe,
		double tqss, std::complex<double> tqsps
		) {
		std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
		std::complex<double> uim = std::complex<double>(0.0, 1.0);
		double two_pi = 4.0 * acos(0.0);
		for (int i80 = 1; i80 <= nsph; i80++) {
		if(c1->iog[i80 - 1] >= i80) {
		tqse[0][i80 - 1] = 0.0;
		tqse[1][i80 - 1] = 0.0;
		tqspe[0][i80 - 1] = cc0;
		tqspe[1][i80 - 1] = cc0;
		tqss[0][i80 - 1] = 0.0;
		tqss[1][i80 - 1] = 0.0;
		tqsps[0][i80 - 1] = cc0;
		tqsps[1][i80 - 1] = cc0;
		for (int l70 = 1; l70 <= li; l70++) {
		double fl = 1.0 * l70 + l70 + 1;
		std::complex<double> rm = 1.0 / c1->rmi[l70 - 1][i80 - 1];
		double rmm = (rm * dconjg(rm)).real();
		std::complex<double> re = 1.0 / c1->rei[l70 - 1][i80 - 1];
		double rem = (re * dconjg(re)).real();
		if (inpol == 0) {
		std::complex<double> pce = ((rm + re) * uim) * fl;
		std::complex<double> pcs = ((rmm + rem) * fl) * uim;
		tqspe[0][i80 - 1] -= pce;
		tqspe[1][i80 - 1] += pce;
		tqsps[0][i80 - 1] -= pcs;
		tqsps[1][i80 - 1] += pcs;
		} else {
		double ce = (rm + re).real() * fl;
		double cs = (rmm + rem) * fl;
		tqse[0][i80 - 1] -= ce;
		tqse[1][i80 - 1] += ce;
		tqss[0][i80 - 1] -= cs;
		tqss[1][i80 - 1] += cs;
		}
		}
		if (inpol == 0) {
		tqspe[0][i80 - 1] *= two_pi;
		tqspe[1][i80 - 1] *= two_pi;
		tqsps[0][i80 - 1] *= two_pi;
		tqsps[1][i80 - 1] *= two_pi;
		printf("DEBUG: TQSPE(1, %d ) = ( %lE, %lE)\n", i80, tqspe[0][i80 - 1].real(), tqspe[0][i80 - 1].imag());
		printf("DEBUG: TQSPE(2, %d ) = ( %lE, %lE)\n", i80, tqspe[1][i80 - 1].real(), tqspe[1][i80 - 1].imag());
		printf("DEBUG: TQSPS(1, %d ) = ( %lE, %lE)\n", i80, tqsps[0][i80 - 1].real(), tqsps[0][i80 - 1].imag());
		printf("DEBUG: TQSPS(2, %d ) = ( %lE, %lE)\n", i80, tqsps[1][i80 - 1].real(), tqsps[1][i80 - 1].imag());
		} else {
		tqse[0][i80 - 1] *= two_pi;
		tqse[1][i80 - 1] *= two_pi;
		tqss[0][i80 - 1] *= two_pi;
		tqss[1][i80 - 1] *= two_pi;
		printf("DEBUG: TQSE(1, %d ) = %lE)\n", i80, tqse[0][i80 - 1]);
		printf("DEBUG: TQSE(2, %d ) = %lE)\n", i80, tqse[1][i80 - 1]);
		printf("DEBUG: TQSS(1, %d ) = %lE)\n", i80, tqss[0][i80 - 1]);
		printf("DEBUG: TQSS(2, %d ) = %lE)\n", i80, tqss[1][i80 - 1]);
		}
		}
		}
		}

		/*! \brief C++ porting of RBF
		*
		* This is the Real Bessel Function.
		@@ -639,6 +821,95 @@ void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri
		}
		}

		/*! \brief C++ porting of SPHAR
		*
		* \param cosrth: `double`
		* \param sinrth: `double`
		* \param cosrph: `double`
		* \param sinrph: `double`
		* \param ll: `int`
		* \param ylm: Vector of complex.
		*/
		void sphar(
		double cosrth, double sinrth, double cosrph, double sinrph,
		int ll, std::complex<double> *ylm
		) {
		double sinrmp[40], cosrmp[40], plegn[861];
		double four_pi = 8.0 * acos(0.0);
		double pi4irs = 1.0 / sqrt(four_pi);
		double x = cosrth;
		double y = sinrth;
		//printf("DEBUG: X = %lE\n", x);
		if (y < 0.0) y *= -1.0;
		//printf("DEBUG: Y = %lE\n", y);
		double cllmo = 3.0;
		double cll = 1.5;
		double ytol = y;
		plegn[0] = 1.0;
		//printf("DEBUG: PLEGN( %d ) = %lE\n", 1, plegn[0]);
		plegn[1] = x * sqrt(cllmo);
		//printf("DEBUG: PLEGN( %d ) = %lE\n", 2, plegn[1]);
		plegn[2] = ytol * sqrt(cll);
		//printf("DEBUG: PLEGN( %d ) = %lE\n", 3, plegn[2]);
		sinrmp[0] = sinrph;
		cosrmp[0] = cosrph;
		int k;
		if (ll >= 2) {
		k = 3;
		for (int l20 = 2; l20 <= ll; l20++) {
		int lmo = l20 - 1;
		int ltpo = l20 + l20 + 1;
		int ltmo = ltpo - 2;
		int lts = ltpo * ltmo;
		double cn = 1.0 * lts;
		for (int mpo10 = 1; mpo10 <= lmo; mpo10++) {
		int m = mpo10 - 1;
		int mpopk = mpo10 + k;
		int ls = (l20 + m) * (l20 - m);
		double cd = 1.0 * ls;
		double cnm = 1.0 * ltpo * (lmo + m) * (l20 - mpo10);
		double cdm = 1.0 * ls * (ltmo - 2);
		plegn[mpopk - 1] = plegn[mpopk - l20 - 1] * x * sqrt(cn / cd) -
		plegn[mpopk - ltmo - 1] * sqrt(cnm / cdm);
		//printf("DEBUG: PLEGN( %d ) = %lE\n", mpopk, plegn[mpopk - 1]);
		}
		int lpk = l20 + k;
		double cltpo = 1.0 * ltpo;
		plegn[lpk - 1] = plegn[k - 1] * x * sqrt(cltpo);
		//printf("DEBUG: PLEGN( %d ) = %lE\n", lpk, plegn[lpk - 1]);
		k = lpk + 1;
		double clt = 1.0 * (ltpo - 1);
		cll *= (cltpo / clt);
		ytol *= y;
		plegn[k - 1] = ytol * sqrt(cll);
		//printf("DEBUG: PLEGN( %d ) = %lE\n", k, plegn[k - 1]);
		sinrmp[l20 - 1] = sinrph * cosrmp[lmo - 1] + cosrph * sinrmp[lmo - 1];
		cosrmp[l20 - 1] = cosrph * cosrmp[lmo - 1] - sinrph * sinrmp[lmo - 1];
		} // end l20 loop
		}
		//l = 0;
		//bool goto40 = true;
		for (int l = 0; l < ll; l++) {
		//int m = 0;
		k = l * (l + 1);
		int l0y = k + 1;
		int l0p = k / 2 + 1;
		ylm[l0y - 1] = pi4irs * plegn[l0p - 1];
		//printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", l0y, ylm[l0y - 1].real(), ylm[l0y - 1].imag());
		for (int m = 0; m < l; m++) {
		int lmp = l0p + m;
		double save = pi4irs * plegn[lmp - 1];
		int lmy = l0y + m;
		ylm[lmy - 1] = save * std::complex<double>(cosrmp[m - 1], sinrmp[m - 1]);
		if (m % 2 != 0) ylm[lmy - 1] *= -1.0;
		//printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", lmy, ylm[lmy - 1].real(), ylm[lmy - 1].imag());
		lmy = l0y - m;
		ylm[lmy - 1] = save * std::complex<double>(cosrmp[m - 1], -sinrmp[m - 1]);
		//printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", lmy, ylm[lmy - 1].real(), ylm[lmy - 1].imag());
		}
		}
		}

		/*! \brief C++ porting of THDPS
		*
		* \param lm: `int`
		@@ -680,6 +951,255 @@ void thdps(int lm, double ****zpv) {
		}
		}

		/*! \brief C++ porting of UPVMP
		*
		* \param thd: `double`
		* \param phd: `double`
		* \param icspnv: `int`
		* \param cost: `double`
		* \param sint: `double`
		* \param cosp: `double`
		* \param sinp: `double`
		* \param u: `double *`
		* \param up: `double *`
		* \param un: `double *`
		*/
		void upvmp(
		double thd, double phd, int icspnv, double &cost, double &sint,
		double &cosp, double &sinp, double u, double up, double *un
		) {
		double half_pi = acos(0.0);
		double rdr = half_pi / 90.0;
		double th = thd * rdr;
		double ph = phd * rdr;
		cost = cos(th);
		sint = sin(th);
		cosp = cos(ph);
		sinp = sin(ph);
		//printf("DEBUG: cost = %lE\n", cost);
		//printf("DEBUG: sint = %lE\n", sint);
		//printf("DEBUG: cosp = %lE\n", cosp);
		//printf("DEBUG: sinp = %lE\n", sinp);
		u[0] = cosp * sint;
		u[1] = sinp * sint;
		u[2] = cost;
		up[0] = cosp * cost;
		up[1] = sinp * cost;
		up[2] = -sint;
		un[0] = -sinp;
		un[1] = cosp;
		un[2] = 0.0;
		if (icspnv != 0) {
		up[0] *= -1.0;
		up[1] *= -1.0;
		up[2] *= -1.0;
		un[0] *= -1.0;
		un[1] *= -1.0;
		}
		}

		/*! \brief C++ porting of UPVSP
		*
		* \param u: `double *`
		* \param upmp: `double *`
		* \param unmp: `double *`
		* \param us: `double *`
		* \param upsmp: `double *`
		* \param unsmp: `double *`
		* \param up: `double *`
		* \param un: `double *`
		* \param ups: `double *`
		* \param uns: `double *`
		* \param duk: `double *`
		* \param isq: `int &`
		* \param ibf: `int &`
		* \param scand: `double &`
		* \param cfmp: `double &`
		* \param sfmp: `double &`
		* \param cfsp: `double &`
		* \param sfsp: `double &`
		*/
		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
		) {
		double rdr = acos(0.0) / 90.0;
		double small = 1.0e-6;
		isq = 0;
		scand = u[0] * us[0] + u[1] * us[1] + u[2] * us[2];
		double abs_scand = scand - 1.0;
		if (abs_scand < 0.0) abs_scand *= -1.0;
		if (abs_scand >= small) {
		abs_scand = scand + 1.0;
		if (abs_scand < 0.0) abs_scand *= -1.0;
		if (abs_scand >= small) {
		scand = acos(scand) / rdr;
		duk[0] = u[0] - us[0];
		duk[1] = u[1] - us[1];
		duk[2] = u[2] - us[2];
		ibf = 0;
		} else {
		// label 15
		scand = 180.0;
		duk[0] = 2.0 * u[0];
		duk[1] = 2.0 * u[1];
		duk[2] = 2.0 * u[2];
		ibf = 1;
		ups[0] = -upsmp[0];
		ups[1] = -upsmp[1];
		ups[2] = -upsmp[2];
		uns[0] = -unsmp[0];
		uns[1] = -unsmp[1];
		uns[2] = -unsmp[2];
		}
		} else {
		// label 10
		scand = 0.0;
		duk[0] = 0.0;
		duk[1] = 0.0;
		duk[2] = 0.0;
		ibf = -1;
		isq = -1;
		ups[0] = upsmp[0];
		ups[1] = upsmp[1];
		ups[2] = upsmp[2];
		uns[0] = unsmp[0];
		uns[1] = unsmp[1];
		uns[2] = unsmp[2];
		}
		if (ibf == -1 \|\| ibf == 1) { // label 20
		up[0] = upmp[0];
		up[1] = upmp[1];
		up[2] = upmp[2];
		un[0] = unmp[0];
		un[1] = unmp[1];
		un[2] = unmp[2];
		} else { // label 25
		orunve(u, us, un, -1, small);
		uns[0] = un[0];
		uns[1] = un[1];
		uns[2] = un[2];
		orunve(un, u, up, 1, small);
		orunve(uns, us, ups, 1, small);
		}
		// label 85
		cfmp = upmp[0] * up[0] + upmp[1] * up[1] + upmp[2] * up[2];
		sfmp = unmp[0] * up[0] + unmp[1] * up[1] + unmp[2] * up[2];
		cfsp = ups[0] * upsmp[0] + ups[1] * upsmp[1] + ups[2] * upsmp[2];
		sfsp = uns[0] * upsmp[0] + uns[1] * upsmp[1] + uns[2] * upsmp[2];
		}

		/*! \brief C++ porting of WMAMP
		*
		* \param iis: `int`
		* \param cost: `double`
		* \param sint: `double`
		* \param cosp: `double`
		* \param sinp: `double`
		* \param inpol: `int`
		* \param lm: `int`
		* \param idot: `int`
		* \param nsph: `int`
		* \param arg: `double *`
		* \param u: `double *`
		* \param up: `double *`
		* \param un: `double *`
		* \param c1: `C1 *`
		*/
		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
		) {
		std::complex<double> *ylm = new std::complex<double>[1682];
		int nlmp = lm * (lm + 2) + 2;
		ylm[nlmp] = std::complex<double>(0.0, 0.0);
		if (idot != 0) {
		if (idot != 1) {
		for (int n40 = 0; n40 < nsph; n40++) {
		arg[n40] = u[0] * c1->rxx[n40] + u[1] * c1->ryy[n40] + u[2] * c1->rzz[n40];
		}
		} else {
		for (int n50 = 0; n50 < nsph; n50++) {
		arg[n50] = c1->rzz[n50];
		}
		}
		if (iis == 2) {
		for (int n60 = 0; n60 < nsph; n60++) arg[n60] *= -1;
		}
		}
		sphar(cost, sint, cosp, sinp, lm, ylm);
		pwma(up, un, ylm, inpol, lm, iis, c1);
		delete[] ylm;
		}

		/*! \brief C++ porting of WMASP
		*
		* \param cost: `double`
		* \param sint: `double`
		* \param cosp: `double`
		* \param sinp: `double`
		* \param costs: `double`
		* \param sints: `double`
		* \param cosps: `double`
		* \param sinps: `double`
		* \param u: `double *`
		* \param up: `double *`
		* \param un: `double *`
		* \param us: `double *`
		* \param ups: `double *`
		* \param uns: `double *`
		* \param isq: `int`
		* \param ibf: `int`
		* \param inpol: `int`
		* \param lm: `int`
		* \param idot: `int`
		* \param nsph: `int`
		* \param argi: `double *`
		* \param args: `double *`
		* \param c1: `C1 *`
		*/
		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
		) {
		std::complex<double> *ylm = new std::complex<double>[1682];
		int nlmp = lm * (lm + 2) + 2;
		ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
		if (idot != 0) {
		if (idot != 1) {
		for (int n40 = 1; n40 <= nsph; n40++) {
		argi[n40 - 1] = u[0] * c1->rxx[n40 - 1] + u[1] * c1->ryy[n40 - 1] + u[2] * c1->rzz[n40 - 1];
		if (ibf != 0) {
		args[n40 - 1] = argi[n40 - 1] * ibf;
		} else {
		args[n40 - 1] = -1.0 * (us[0] * c1->rxx[n40 - 1] + us[1] * c1->ryy[n40 - 1] + us[2] * c1->rzz[n40 - 1]);
		}
		}
		} else { // label 50
		for (int n60 = 1; n60 <= nsph; n60++) {
		argi[n60 - 1] = cost * c1->rzz[n60 - 1];
		if (ibf != 0) {
		args[n60 - 1] = argi[n60 - 1] * ibf;
		} else {
		args[n60 - 1] = -costs * c1->rzz[n60 - 1];
		}
		}
		}
		}
		sphar(cost, sint, cosp, sinp, lm, ylm);
		pwma(up, un, ylm, inpol, lm, isq, c1);
		if (ibf >= 0) {
		sphar(costs, sints, cosps, sinps, lm, ylm);
		pwma(ups, uns, ylm, inpol, lm, 2, c1);
		}
		delete[] ylm;
		}


		/*! \brief C++ porting of DME
		*
		* \param li: `int`

src/sphere/sphere.cpp

+218 −3

File changed.

Preview size limit exceeded, changes collapsed.