Complete the migration of sph.f to C++ (432f5809) · Commits · Giacomo Mulas / NP_TMcode

src/include/sph_subs.h

+206 −47

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

		/*! \brief C++ porting of MMULC
		*
		* \param vint: Vector of complex.
		* \param cmullr: `double **`
		* \param cmul: `double **`
		*/
		void mmulc(std::complex<double> vint, double cmullr, double *cmul) {
		double sm2 = vint[0].real();
		double s24 = vint[1].real();
		double d24 = vint[1].imag();
		double sm4 = vint[5].real();
		double s23 = vint[8].real();
		double d32 = vint[8].imag();
		double s34 = vint[9].real();
		double d34 = vint[9].imag();
		double sm3 = vint[10].real();
		double s31 = vint[11].real();
		double d31 = vint[11].imag();
		double s21 = vint[12].real();
		double d12 = vint[12].imag();
		double s41 = vint[13].real();
		double d14 = vint[13].imag();
		double sm1 = vint[15].real();
		cmullr[0][0] = sm2;
		cmullr[0][1] = sm3;
		cmullr[0][2] = -s23;
		cmullr[0][3] = -d32;
		cmullr[1][0] = sm4;
		cmullr[1][1] = sm1;
		cmullr[1][2] = -s41;
		cmullr[1][3] = -d14;
		cmullr[2][0] = -s24 * 2.0;
		cmullr[2][1] = -s31 * 2.0;
		cmullr[2][2] = s21 + s34;
		cmullr[2][3] = d34 + d12;
		cmullr[3][0] = -d24 * 2.0;
		cmullr[3][1] = -d31 * 2.0;
		cmullr[3][2] = d34 - d12;
		cmullr[3][3] = s21 - s34;
		cmul[0][0] = (sm2 + sm3 + sm4 + sm1) * 0.5;
		cmul[0][1] = (sm2 - sm3 + sm4 - sm1) * 0.5;
		cmul[0][2] = -s23 - s41;
		cmul[0][3] = -d32 - d14;
		cmul[1][0] = (sm2 + sm3 - sm4 - sm1) * 0.5;
		cmul[1][1] = (sm2 - sm3 - sm4 + sm1) * 0.5;
		cmul[1][2] = -s23 + s41;
		cmul[1][3] = -d32 + d14;
		cmul[2][0] = -s24 - s31;
		cmul[2][1] = -s24 + s31;
		cmul[2][2] = s21 + s34;
		cmul[2][3] = d34 + d12;
		cmul[3][0] = -d24 - d31;
		cmul[3][1] = -d24 + d31;
		cmul[3][2] = d34 - d12;
		cmul[3][3] = s21 - s34;
		}

		/*! \brief C++ porting of ORUNVE
		*
		* \param u1: `double *`
		@@ -418,11 +475,14 @@ void pwma(
		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];
		//printf("DEBUG: in PWMA YLM(2) = (%lE, %lE)\n", ylm[1].real(), ylm[1].imag());
		for (int l20 = 1; l20 <= lw; l20++) {
		//if (ispo == 1) printf("DEBUG: l20 = %d\n", 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);
		std::complex<double> cl = std::complex<double>(four_pi / sqrt(x), 0.0);
		for (int powi = 1; powi <= l20; powi++) cl *= uim;
		int mv = l20 + lf;
		int m = -lf;
		for (int mf20 = 1; mf20 <= mv; mf20++) {
		@@ -438,16 +498,21 @@ void pwma(
		cm1 * cm * ylm[k - 1] +
		cz1 * cz * ylm[k]
		) * cl;
		//if (ispo == 1) printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k, ispo, c1->w[k - 1][ispo - 1].real(), c1->w[k - 1][ispo - 1].imag());
		c1->w[k - 1][ispt - 1] = dconjg(
		cp2 * cp * ylm[k + 1] +
		cm2 * cm * ylm[k - 1] +
		cz2 * cz * ylm[k]
		) * cl;
		//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k, ispt, c1->w[k - 1][ispt - 1].real(), c1->w[k - 1][ispt - 1].imag());
		}
		}
		for (int k30 = 1; k30 <= nlwm; k30++) {
		int i = k30 + nlwm;
		c1->w[i - 1][ispo - 1] = uim * c1->w[k30 - 1][ispt - 1];
		//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispo, c1->w[i - 1][ispo - 1].real(), c1->w[i - 1][ispo - 1].imag());
		c1->w[i - 1][ispt - 1] = -uim * c1->w[k30 - 1][ispo - 1];
		//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispt, c1->w[i - 1][ispt - 1].real(), c1->w[i - 1][ispt - 1].imag());
		}
		if (inpol != 0) {
		for (int k40 = 1; k40 <= nlwm; k40++) {
		@@ -455,20 +520,30 @@ void pwma(
		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;
		//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k40, ispo, c1->w[k40 - 1][ispo - 1].real(), c1->w[k40 - 1][ispo - 1].imag());
		c1->w[i - 1][ispo - 1] = -cc2;
		//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispo, c1->w[i - 1][ispo - 1].real(), c1->w[i - 1][ispo - 1].imag());
		c1->w[k40 - 1][ispt - 1] = cc1;
		//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k40, ispt, c1->w[k40 - 1][ispt - 1].real(), c1->w[k40 - 1][ispt - 1].imag());
		c1->w[i - 1][ispt - 1] = cc1;
		//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispt, c1->w[i - 1][ispt - 1].real(), c1->w[i - 1][ispt - 1].imag());
		}
		} else {
		if (isq == 0) {
		return;
		}
		}
		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++)
		for (int k50 = 1; k50 <= nlwmt; k50++) {
		c1->w[k50 - 1][ipt - 1] = dconjg(c1->w[k50 - 1][ipis - 1]);
		//printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k50, ipt, c1->w[k50 - 1][ipt - 1].real(), c1->w[k50 - 1][ipt - 1].imag());
		}
		}
		}
		//printf("DEBUG: out PWMA W(1,1) = (%lE, %lE)\n", c1->w[0][0].real(), c1->w[0][0].imag());
		}

		/*! \brief C++ porting of RABAS
		@@ -783,8 +858,8 @@ void rnf(int n, double x, int &nm, double sy[]) {
		* \param c1: `C1 *`
		*/
		void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
		std::complex<double> sum21, rm, re, cc0, csam;
		cc0 = std::complex<double>(0.0, 0.0);
		std::complex<double> sum21, rm, re, csam;
		std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
		double exdc = exri * exri;
		double ccs = 4.0 * acos(0.0) / (vk * vk);
		double cccs = ccs / exdc;
		@@ -821,6 +896,79 @@ void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri
		}
		}

		/*! \brief C++ porting of SSCR2
		*
		* \param nsph: `int`
		* \param lm: `int`
		* \param vk: `double`
		* \param exri: `double`
		* \param c1: `C1 *`
		*/
		void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
		std::complex<double> s11, s21, s12, s22, rm, re, csam, cc;
		std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
		double ccs = 1.0 / (vk * vk);
		csam = -(ccs / (exri * vk)) * std::complex<double>(0.0, 0.5);
		double pigfsq = 64.0 * acos(0.0) * acos(0.0);
		double cfsq = 4.0 / (pigfsq * ccs * ccs);
		int nlmm = lm * (lm + 2);
		//printf("DEBUG: in SSCR2 W(1,1) = (%lE,%lE)\n", c1->w[0][0].real(), c1->w[0][0].imag());
		for (int i14 = 1; i14 <= nsph; i14++) {
		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 <= lm; l10++) {
		rm = 1.0 / c1->rmi[l10 - 1][i14 - 1];
		re = 1.0 / c1->rei[l10 - 1][i14 - 1];
		//printf("DEBUG: rm = (%lE,%lE)\n", rm.real(), rm.imag());
		//printf("DEBUG: re = (%lE,%lE)\n", re.real(), re.imag());
		int ltpo = l10 + l10 + 1;
		for (int im10 = 1; im10 <= ltpo; im10++) {
		k += 1;
		int ke = k + nlmm;
		//printf("DEBUG: W( %d, 3) = (%lE,%lE)\n", k, c1->w[k - 1][2].real(), c1->w[k - 1][2].imag());
		//printf("DEBUG: W( %d, 1) = (%lE,%lE)\n", k, c1->w[k - 1][0].real(), c1->w[k - 1][0].imag());
		//printf("DEBUG: W( %d, 3) = (%lE,%lE)\n", ke, c1->w[ke - 1][2].real(), c1->w[ke - 1][2].imag());
		//printf("DEBUG: W( %d, 1) = (%lE,%lE)\n", ke, c1->w[ke - 1][0].real(), c1->w[ke - 1][0].imag());
		s11 = s11 - c1->w[k - 1][2] * c1->w[k - 1][0] * rm - c1->w[ke - 1][2] * c1->w[ke - 1][0] * re;
		//printf("DEBUG: W( %d, 4) = (%lE,%lE)\n", k, c1->w[k - 1][3].real(), c1->w[k - 1][3].imag());
		//printf("DEBUG: W( %d, 4) = (%lE,%lE)\n", ke, c1->w[ke - 1][3].real(), c1->w[ke - 1][3].imag());
		s21 = s21 - c1->w[k - 1][3] * c1->w[k - 1][0] * rm - c1->w[ke - 1][3] * c1->w[ke - 1][0] * re;
		//printf("DEBUG: W( %d, 2) = (%lE,%lE)\n", k, c1->w[k - 1][1].real(), c1->w[k - 1][1].imag());
		//printf("DEBUG: W( %d, 2) = (%lE,%lE)\n", ke, c1->w[ke - 1][1].real(), c1->w[ke - 1][1].imag());
		s12 = s12 - c1->w[k - 1][2] * c1->w[k - 1][1] * rm - c1->w[ke - 1][2] * c1->w[ke - 1][1] * re;
		s22 = s22 - c1->w[k - 1][3] * c1->w[k - 1][1] * rm - c1->w[ke - 1][3] * c1->w[ke - 1][1] * re;
		}
		}
		c1->sas[i14 - 1][0][0] = s11 * csam;
		c1->sas[i14 - 1][1][0] = s21 * csam;
		c1->sas[i14 - 1][0][1] = s12 * csam;
		c1->sas[i14 - 1][1][1] = s22 * csam;
		}
		} // loop i14
		for (int i24 = 1; i24 <= nsph; i24++) {
		int iogi = c1->iog[i24 - 1];
		if (iogi >= i24) {
		int j = 0;
		for (int ipo1 = 0; ipo1 < 2; ipo1++) {
		for (int jpo1 = 0; jpo1 < 2; jpo1++) {
		std::complex<double> cc = dconjg(c1->sas[i24 - 1][jpo1][ipo1]);
		for (int ipo2 = 0; ipo2 < 2; ipo2++) {
		for (int jpo2 = 0; jpo2 < 2; jpo2++) {
		j += 1;
		c1->vints[i24 - 1][j - 1] = c1->sas[i24 - 1][jpo2][ipo2] * cc * cfsq;
		}
		}
		}
		}
		}
		}
		}

		/*! \brief C++ porting of SPHAR
		*
		* \param cosrth: `double`
		@@ -853,9 +1001,8 @@ void sphar(
		//printf("DEBUG: PLEGN( %d ) = %lE\n", 3, plegn[2]);
		sinrmp[0] = sinrph;
		cosrmp[0] = cosrph;
		int k;
		if (ll >= 2) {
		k = 3;
		int k = 3;
		for (int l20 = 2; l20 <= ll; l20++) {
		int lmo = l20 - 1;
		int ltpo = l20 + l20 + 1;
		@@ -887,27 +1034,36 @@ void sphar(
		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;
		// label 30
		int l = 0;
		int m, k, l0y, l0p, lmy, lmp;
		double save;
		label40:
		m = 0;
		k = l * (l + 1);
		int l0y = k + 1;
		int l0p = k / 2 + 1;
		l0y = k + 1;
		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;
		goto label45;
		label44:
		lmp = l0p + m;
		save = pi4irs * plegn[lmp - 1];
		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());
		}
		}
		label45:
		if (m >= l) goto label47;
		m += 1;
		goto label44;
		label47:
		if (l >= ll) return;
		l += 1;
		goto label40;
		}

		/*! \brief C++ porting of THDPS
		@@ -1130,6 +1286,7 @@ void wmamp(
		}
		}
		sphar(cost, sint, cosp, sinp, lm, ylm);
		//printf("DEBUG: in WMAMP and calling PWMA with lm = %d and iis = %d\n", lm, iis);
		pwma(up, un, ylm, inpol, lm, iis, c1);
		delete[] ylm;
		}
		@@ -1191,9 +1348,11 @@ void wmasp(
		}
		}
		sphar(cost, sint, cosp, sinp, lm, ylm);
		//printf("DEBUG: in WMASP and calling PWMA with isq = %d\n", isq);
		pwma(up, un, ylm, inpol, lm, isq, c1);
		if (ibf >= 0) {
		sphar(costs, sints, cosps, sinps, lm, ylm);
		//printf("DEBUG: in WMASP and calling PWMA with isq = 2 and ibf = %d\n", ibf);
		pwma(ups, uns, ylm, inpol, lm, 2, c1);
		}
		delete[] ylm;

src/sphere/sph.f

+61 −7

Original line number	Diff line number	Diff line
		@@ -456,11 +456,19 @@ CCC 1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
		TQSPE(2,I)=TQSPE(2,I)*PIG2
		TQSPS(1,I)=TQSPS(1,I)*PIG2
		TQSPS(2,I)=TQSPS(2,I)*PIG2
		PRINT *,"DEBUG: TQSPE(1,",I,") =",TQSPE(1,I)
		PRINT *,"DEBUG: TQSPE(2,",I,") =",TQSPE(2,I)
		PRINT *,"DEBUG: TQSPS(1,",I,") =",TQSPS(1,I)
		PRINT *,"DEBUG: TQSPS(2,",I,") =",TQSPS(2,I)
		GO TO 80
		75 TQSE(1,I)=TQSE(1,I)*PIG2
		TQSE(2,I)=TQSE(2,I)*PIG2
		TQSS(1,I)=TQSS(1,I)*PIG2
		TQSS(2,I)=TQSS(2,I)*PIG2
		PRINT *,"DEBUG: TQSE(1,",I,") =",TQSE(1,I)
		PRINT *,"DEBUG: TQSE(2,",I,") =",TQSE(2,I)
		PRINT *,"DEBUG: TQSS(1,",I,") =",TQSS(1,I)
		PRINT *,"DEBUG: TQSS(2,",I,") =",TQSS(2,I)
		80 CONTINUE
		RETURN
		END
		@@ -576,6 +584,7 @@ CCC 1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
		PIGFSQ=6.4D+1DACOS(C0)*2
		CFSQ=4.0D0/(PIGFSQCCSCCS)
		NLMM=LM*(LM+2)
		C PRINT *,"DEBUG: in SSCR2 W(1,1) =",W(1,1)
		DO 14 I=1,NSPH
		IOGI=IOG(I)
		IF(IOGI.LT.I)GO TO 14
		@@ -587,14 +596,30 @@ CCC 1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
		DO 10 L=1,LM
		RM=1.0D0/RMI(L,I)
		RE=1.0D0/REI(L,I)
		C PRINT *,"DEBUG: RM =",RM
		C PRINT *,"DEBUG: RE =",RE
		LTPO=L+L+1
		DO 10 IM=1,LTPO
		K=K+1
		KE=K+NLMM
		C PRINT *,"DEBUG: W(",K,",3) =",W(K,3)
		C PRINT *,"DEBUG: W(",K,",1) =",W(K,1)
		C PRINT *,"DEBUG: W(",KE,",3) =",W(KE,3)
		C PRINT *,"DEBUG: W(",KE,",1) =",W(KE,1)
		S11=S11-W(K,3)W(K,1)RM-W(KE,3)W(KE,1)RE
		C PRINT *,"DEBUG: S11 =",S11
		C PRINT *,"DEBUG: W(",K,",4) =",W(K,4)
		C PRINT *,"DEBUG: W(",KE,",4) =",W(KE,4)
		S21=S21-W(K,4)W(K,1)RM-W(KE,4)W(KE,1)RE
		C PRINT *,"DEBUG: W(",K,",2) =",W(K,2)
		C PRINT *,"DEBUG: W(",KE,",2) =",W(KE,2)
		S12=S12-W(K,3)W(K,2)RM-W(KE,3)W(KE,2)RE
		10 S22=S22-W(K,4)W(K,2)RM-W(KE,4)W(KE,2)RE
		C PRINT *,"DEBUG: CSAM =",CSAM
		C PRINT *,"DEBUG: S11 =",S11
		C PRINT *,"DEBUG: S21 =",S21
		C PRINT *,"DEBUG: S12 =",S12
		C PRINT *,"DEBUG: S22 =",S22
		SAS(I,1,1)=S11*CSAM
		SAS(I,2,1)=S21*CSAM
		SAS(I,1,2)=S12*CSAM
		@@ -756,6 +781,10 @@ CCC CG1(LMPML,MU,L,M)=CLGO(1,LMP,L;MU,M-MU,M)
		SINT=DSIN(TH)
		COSP=DCOS(PH)
		SINP=DSIN(PH)
		C PRINT *,"DEBUG: cost =",COST
		C PRINT *,"DEBUG: sint =",SINT
		C PRINT *,"DEBUG: cosp =",COSP
		C PRINT *,"DEBUG: sinp =",SINP
		U(1)=COSP*SINT
		U(2)=SINP*SINT
		U(3)=COST
		@@ -799,6 +828,8 @@ CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
		DO 60 N=1,NSPH
		60 ARG(N)=-ARG(N)
		65 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
		C PRINT *,"DEBUG: in WMAMP and calling PWMA with lm =",LM,
		C 1"and iis =",IIS
		CALL PWMA(UP,UN,YLM,INPOL,LM,IIS)
		RETURN
		END
		@@ -897,9 +928,11 @@ CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
		55 ARGS(N)=-COSTS*RZZ(N)
		60 CONTINUE
		75 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
		C PRINT *,"DEBUG: in WMASP and calling PWMA"
		CALL PWMA(UP,UN,YLM,INPOL,LM,ISQ)
		IF(IBF.LT.0)RETURN
		CALL SPHAR(COSTS,SINTS,COSPS,SINPS,LM,YLM)
		C PRINT *,"DEBUG: in WMASP and calling PWMA with IBF =",IBF
		CALL PWMA(UPS,UNS,YLM,INPOL,LM,2)
		RETURN
		END
		@@ -930,6 +963,7 @@ CCC DIMENSION YLM(NLWM+2)
		CM2=.5D0*DCMPLX(UN(1),UN(2))
		CP2=.5D0*DCMPLX(UN(1),-UN(2))
		CZ2=UN(3)
		C PRINT *,"DEBUG: in PWMA YLM(2) =",YLM(2)
		DO 20 L=1,LW
		LF=L+1
		LFTL=LF*L
		@@ -946,26 +980,35 @@ CCC DIMENSION YLM(NLWM+2)
		CM=DSQRT(X)
		CZ=M
		W(K,ISPO)=DCONJG(CP1CPYLM(K+2)+CM1CMYLM(K)+CZ1CZYLM(K+1))*CL
		C IF(ISPO.EQ.1)PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO)
		20 W(K,ISPT)=DCONJG(CP2CPYLM(K+2)+CM2CMYLM(K)+CZ2CZYLM(K+1))*CL
		C 20 PRINT *,"DEBUG: W(",K,",",ISPT,") =",W(K,ISPT)
		DO 30 K=1,NLWM
		I=K+NLWM
		W(I,ISPO)=UIM*W(K,ISPT)
		C PRINT *,"DEBUG: W(",I,",",ISPO,") =",W(I,ISPO)
		30 W(I,ISPT)=-UIM*W(K,ISPO)
		C 30 PRINT *,"DEBUG: W(",I,",",ISPT,") =",W(I,ISPT)
		IF(INPOL.EQ.0)GO TO 42
		DO 40 K=1,NLWM
		I=K+NLWM
		C1=(W(K,ISPO)+UIMW(K,ISPT))SQRTWI
		C2=(W(K,ISPO)-UIMW(K,ISPT))SQRTWI
		W(K,ISPO)=C2
		C PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO)
		W(I,ISPO)=-C2
		C PRINT *,"DEBUG: W(",I,",",ISPO,") =",W(I,ISPO)
		W(K,ISPT)=C1
		C PRINT *,"DEBUG: W(",K,",",ISPT,") =",W(K,ISPT)
		40 W(I,ISPT)=C1
		C 40 PRINT *,"DEBUG: W(",I,",",ISPT,") =",W(I,ISPT)
		42 IF(ISQ.EQ.0)RETURN
		DO 50 I=1,2
		IPT=I+2
		IPIS=I+IS
		DO 50 K=1,NLWMT
		50 W(K,IPT)=DCONJG(W(K,IPIS))
		C 50 PRINT *,"DEBUG: W(",K,",",IPT,") =",W(K,IPT)
		RETURN
		END
		SUBROUTINE ORUNVE(U1,U2,U3,IORTH,TORTH)
		@@ -1455,13 +1498,18 @@ CCC LL=LM
		PI4=DACOS(0.0D0)*8.0D0
		PI4IRS=1.0D0/DSQRT(PI4)
		X=COSRTH
		C PRINT *,"DEBUG: X =",X
		Y=DABS(SINRTH)
		C PRINT *,"DEBUG: Y =",Y
		CLLMO=3.0D0
		CLL=1.5D0
		YTOL=Y
		PLEGN(1)=1.0D0
		C PRINT *,"DEBUG: PLEGN( 1 ) =",PLEGN(1)
		PLEGN(2)=X*DSQRT(CLLMO)
		C PRINT *,"DEBUG: PLEGN( 2 ) =",PLEGN(2)
		PLEGN(3)=YTOL*DSQRT(CLL)
		C PRINT *,"DEBUG: PLEGN( 3 ) =",PLEGN(3)
		SINRMP(1)=SINRPH
		COSRMP(1)=COSRPH
		IF(LL.LT.2)GO TO 30
		@@ -1481,14 +1529,17 @@ CCC LL=LM
		CDM=(LTMO-2)*LS
		10 PLEGN(MPOPK)=PLEGN(MPOPK-L)XDSQRT(CN/CD)-
		1PLEGN(MPOPK-LTMO)*DSQRT(CNM/CDM)
		C10 PRINT *,"DEBUG: PLEGN(",MPOPK,") =",PLEGN(MPOPK)
		LPK=L+K
		CLTPO=LTPO
		PLEGN(LPK)=PLEGN(K)XDSQRT(CLTPO)
		C PRINT *,"DEBUG: PLEGN(",LPK,") =",PLEGN(LPK)
		K=LPK+1
		CLT=LTPO-1
		CLL=CLL*(CLTPO/CLT)
		YTOL=YTOL*Y
		PLEGN(K)=YTOL*DSQRT(CLL)
		C PRINT *,"DEBUG: PLEGN(",K,") =",PLEGN(K)
		SINRMP(L)=SINRPHCOSRMP(LMO)+COSRPHSINRMP(LMO)
		20 COSRMP(L)=COSRPHCOSRMP(LMO)-SINRPHSINRMP(LMO)
		30 L=0
		@@ -1497,14 +1548,17 @@ CCC LL=LM
		L0Y=K+1
		L0P=K/2+1
		YLM(L0Y)=PI4IRS*PLEGN(L0P)
		C PRINT *, "DEBUG: YLM(",L0Y,") =",YLM(L0Y)
		GO TO 45
		44 LMP=L0P+M
		SAVE=PI4IRS*PLEGN(LMP)
		LMY=L0Y+M
		YLM(LMY)=SAVE*DCMPLX(COSRMP(M),SINRMP(M))
		IF(MOD(M,2).NE.0)YLM(LMY)=-YLM(LMY)
		C PRINT *, "DEBUG: YLM(",LMY,") =",YLM(LMY)
		LMY=L0Y-M
		YLM(LMY)=SAVE*DCMPLX(COSRMP(M),-SINRMP(M))
		C PRINT *, "DEBUG: YLM(",LMY,") =",YLM(LMY)
		45 IF(M.GE.L)GO TO 47
		M=M+1
		GO TO 44

src/sphere/sphere.cpp

+105 −48

Original line number	Diff line number	Diff line
		@@ -23,6 +23,31 @@ void sphere() {
		ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF");
		GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/sphere/DSPH");
		if (sconf->number_of_spheres == gconf->number_of_spheres) {
		int isq, ibf;
		double cost, sint, cosp, sinp;
		double costs, sints, cosps, sinps;
		double scan;
		double *duk = new double[3];
		double *u = new double[3];
		double *us = new double[3];
		double *un = new double[3];
		double *uns = new double[3];
		double *up = new double[3];
		double *ups = new double[3];
		double *upmp = new double[3];
		double *upsmp = new double[3];
		double *unmp = new double[3];
		double *unsmp = new double[3];
		double *cmul = new double[4];
		double *cmullr = new double[4];
		for (int i = 0; i < 4; i++) {
		cmul[i] = new double[4];
		cmullr[i] = new double[4];
		}
		double frx = 0.0, fry = 0.0, frz = 0.0;
		double cfmp, cfsp, sfmp, sfsp;
		complex<double> *vint = new complex<double>[16];
		int jw;
		int nsph = gconf->number_of_spheres;
		C1 *c1 = new C1;
		for (int i = 0; i < nsph; i++) {
		@@ -114,10 +139,10 @@ void sphere() {
		int nk = nth * nph;
		int nks = nths * nphs;
		int nkks = nk * nks;
		int th1 = gconf->in_theta_start;
		int ph1 = gconf->in_phi_start;
		int ths1 = gconf->sc_theta_start;
		int phs1 = gconf->sc_phi_start;
		double th1 = gconf->in_theta_start;
		double ph1 = gconf->in_phi_start;
		double ths1 = gconf->sc_theta_start;
		double phs1 = gconf->sc_phi_start;
		const double half_pi = acos(0.0);
		const double pi = 2.0 * half_pi;
		double gcs = 0.0;
		@@ -144,26 +169,6 @@ void sphere() {
		}
		}
		thdps(gconf->l_max, zpv);
		//DEBUG CODE
		/*
		for (int zi = 0; zi < gconf->l_max; zi++) {
		for (int zj = 0; zj < 3; zj++) {
		for (int zk = 0; zk < 2; zk++) {
		if (zpv[zi][zj][zk][0] != 0.0)
		printf(
		"DEBUG: zpv[%d][%d][%d][0] = %.3lE\n",
		zi, zj, zk, zpv[zi][zj][zk][0]
		);
		if (zpv[zi][zj][zk][1] != 0.0)
		printf(
		"DEBUG: zpv[%d][%d][%d][1] = %.3lE\n",
		zi, zj, zk, zpv[zi][zj][zk][1]
		);
		}
		}
		}
		*/
		//END OF DEBUG CODE
		double exri = sqrt(sconf->exdc);
		fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
		fstream tppoan;
		@@ -265,7 +270,6 @@ void sphere() {
		fclose(output);
		}
		}
		// We are at line 271 of SPH.f
		double cs0 = 0.25 * vk * vk * vk / half_pi;
		sscr0(tfsas, nsph, gconf->l_max, vk, exri, c1);
		printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag());
		@@ -365,27 +369,9 @@ void sphere() {
		qschu, pschu, s0mag
		);
		}
		th = gconf->in_theta_start;
		int isq, ibf;
		double cost, sint, cosp, sinp;
		double costs, sints, cosps, sinps;
		double scan;
		double *duk = new double[3];
		double *u = new double[3];
		double *us = new double[3];
		double *un = new double[3];
		double *uns = new double[3];
		double *up = new double[3];
		double *ups = new double[3];
		double *upmp = new double[3];
		double *upsmp = new double[3];
		double *unmp = new double[3];
		double *unsmp = new double[3];
		double frx, fry, frz;
		double cfmp, cfsp, sfmp, sfsp;
		int jw;
		th = th1;
		for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP parallelizable section
		ph = gconf->in_phi_start;
		ph = ph1;
		for (int jph484 = 1; jph484 <= nph; jph484++) {
		bool goto182 = (nk == 1) && (jxi488 > 1);
		if (!goto182) {
		@@ -397,7 +383,7 @@ void sphere() {
		double rapr = c1->sexs[i183] - gaps[i183];
		frx = rapr * u[0];
		fry = rapr * u[1];
		frx = rapr * u[2];
		frz = rapr * u[2];
		}
		jw = 1;
		}
		@@ -420,7 +406,7 @@ void sphere() {
		phs = ph + phsph;
		if (phs >= 360.0) phs -= 360.0;
		}
		bool goto190 = ((nks == 1) && (jxi488 > 1)) \|\| (jth486 > 1) \|\| (jph484 > 1);
		bool goto190 = (nks == 1) && ((jxi488 > 1) \|\| (jth486 > 1) \|\| (jph484 > 1));
		if (!goto190) {
		upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
		if (gconf->meridional_type >= 0) {
		@@ -472,7 +458,55 @@ void sphere() {
		} else {
		fprintf(output, " UN=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
		}

		sscr2(nsph, gconf->l_max, vk, exri, c1);
		for (int ns226 = 1; ns226 <= nsph; ns226++) {
		fprintf(output, " SPHERE %2d\n", ns226);
		fprintf(
		output, " SAS(1,1)=%15.7lE,%15.7lE, SAS(2,1)=%15.7lE,%15.7lE\n",
		c1->sas[ns226 - 1][0][0].real(), c1->sas[ns226 - 1][0][0].imag(),
		c1->sas[ns226 - 1][1][0].real(), c1->sas[ns226 - 1][1][0].imag()
		);
		fprintf(
		output, " SAS(1,2)=%15.7lE,%15.7lE, SAS(2,2)=%15.7lE,%15.7lE\n",
		c1->sas[ns226 - 1][0][1].real(), c1->sas[ns226 - 1][0][1].imag(),
		c1->sas[ns226 - 1][1][1].real(), c1->sas[ns226 - 1][1][1].imag()
		);
		if (jths482 == 1 && jphs480 == 1)
		fprintf(
		output, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
		frx, fry, frz
		);
		for (int i225 = 0; i225 < 16; i225++)
		vint[i225] = c1->vints[ns226 - 1][i225];
		mmulc(vint, cmullr, cmul);
		fprintf(output, " MULS\n ");
		for (int imul = 0; imul < 4; imul++) {
		for (int jmul = 0; jmul < 4; jmul++) {
		fprintf(output, "%15.7lE", cmul[imul][jmul]);
		}
		if (imul < 3) fprintf(output, "\n ");
		else fprintf(output, "\n");
		}
		fprintf(output, " MULSLR\n ");
		for (int imul = 0; imul < 4; imul++) {
		for (int jmul = 0; jmul < 4; jmul++) {
		fprintf(output, "%15.7lE", cmullr[imul][jmul]);
		}
		if (imul < 3) fprintf(output, "\n ");
		else fprintf(output, "\n");
		}
		for (int vi = 0; vi < 16; vi++) {
		double value = vint[vi].real();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		value = vint[vi].imag();
		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
		}
		for (int imul = 0; imul < 4; imul++) {
		for (int jmul = 0; jmul < 4; jmul++) {
		tppoan.write(reinterpret_cast<char *>(&(cmul[imul][jmul])), sizeof(double));
		}
		}
		} // ns226 loop
		if (gconf->meridional_type < 1) phs += gconf->sc_phi_step;
		} // jphs480 loop
		if (gconf->meridional_type <= 1) thsl += gconf->sc_theta_step;
		@@ -487,6 +521,29 @@ void sphere() {
		printf("ERROR: failed to open TPPOAN file.\n");
		}
		fclose(output);
		delete c1;
		delete c2;
		delete[] duk;
		delete[] u;
		delete[] us;
		delete[] un;
		delete[] uns;
		delete[] up;
		delete[] ups;
		delete[] upmp;
		delete[] upsmp;
		delete[] unmp;
		delete[] unsmp;
		delete[] vint;
		delete[] argi;
		delete[] args;
		delete[] gaps;
		for (int i = 0; i < 4; i++) {
		delete[] cmul[i];
		delete[] cmullr[i];
		}
		delete[] cmul;
		delete[] cmullr;
		} else { // NSPH mismatch between geometry and scatterer configurations.
		throw UnrecognizedConfigurationException(
		"Inconsistent geometry and scatterer configurations."