From 432f58092ba9c36c5e7fdd19c646983253a1c2f2 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Tue, 31 Oct 2023 12:41:10 +0100
Subject: [PATCH] Complete the migration of sph.f to C++
---
src/include/sph_subs.h | 253 +++++++++++++++++++++++++++++++++--------
src/sphere/sph.f | 68 +++++++++--
src/sphere/sphere.cpp | 153 +++++++++++++++++--------
3 files changed, 372 insertions(+), 102 deletions(-)
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index 49bd8dcc..3cdfb083 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -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++) {
@@ -430,7 +490,7 @@ void pwma(
int k = lftl + m;
x = 1.0 * (lftl - m * (m + 1));
double cp = sqrt(x);
- x = 1.0 * (lftl - m * (m -1));
+ x = 1.0 * (lftl - m * (m - 1));
double cm = sqrt(x);
double cz = 1.0 * m;
c1->w[k - 1][ispo - 1] = dconjg(
@@ -438,37 +498,52 @@ 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];
- c1->w[i - 1][ispt - 1] = -uim * c1->w[k30 - 1][ispo - 1];
+ }
+ 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++) {
+ 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;
+ //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());
}
- 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;
- }
+ } 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++)
- c1->w[k50 - 1][ipt - 1] = dconjg(c1->w[k50 - 1][ipis - 1]);
+ }
+ 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]);
+ //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;
- 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());
- }
- }
+ // label 30
+ int l = 0;
+ int m, k, l0y, l0p, lmy, lmp;
+ double save;
+label40:
+ m = 0;
+ k = l * (l + 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());
+ 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;
diff --git a/src/sphere/sph.f b/src/sphere/sph.f
index 77a6df97..bebacd74 100644
--- a/src/sphere/sph.f
+++ b/src/sphere/sph.f
@@ -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+1*DACOS(C0)**2
CFSQ=4.0D0/(PIGFSQ*CCS*CCS)
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
@@ -798,7 +827,9 @@ CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
55 IF(IIS.NE.2)GO TO 65
DO 60 N=1,NSPH
60 ARG(N)=-ARG(N)
- 65 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
+ 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
@@ -896,10 +927,12 @@ CCC NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
GO TO 60
55 ARGS(N)=-COSTS*RZZ(N)
60 CONTINUE
- 75 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
+ 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(CP1*CP*YLM(K+2)+CM1*CM*YLM(K)+CZ1*CZ*YLM(K+1))*CL
- 20 W(K,ISPT)=DCONJG(CP2*CP*YLM(K+2)+CM2*CM*YLM(K)+CZ2*CZ*YLM(K+1))*CL
+C IF(ISPO.EQ.1)PRINT *,"DEBUG: W(",K,",",ISPO,") =",W(K,ISPO)
+ 20 W(K,ISPT)=DCONJG(CP2*CP*YLM(K+2)+CM2*CM*YLM(K)+CZ2*CZ*YLM(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)
- 30 W(I,ISPT)=-UIM*W(K,ISPO)
+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)+UIM*W(K,ISPT))*SQRTWI
C2=(W(K,ISPO)-UIM*W(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
- 40 W(I,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))
+ 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
@@ -1479,16 +1527,19 @@ CCC LL=LM
CD=LS
CNM=LTPO*(LMO+M)*(L-MPO)
CDM=(LTMO-2)*LS
- 10 PLEGN(MPOPK)=PLEGN(MPOPK-L)*X*DSQRT(CN/CD)-
+ 10 PLEGN(MPOPK)=PLEGN(MPOPK-L)*X*DSQRT(CN/CD)-
1PLEGN(MPOPK-LTMO)*DSQRT(CNM/CDM)
+C10 PRINT *,"DEBUG: PLEGN(",MPOPK,") =",PLEGN(MPOPK)
LPK=L+K
CLTPO=LTPO
PLEGN(LPK)=PLEGN(K)*X*DSQRT(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)=SINRPH*COSRMP(LMO)+COSRPH*SINRMP(LMO)
20 COSRMP(L)=COSRPH*COSRMP(LMO)-SINRPH*SINRMP(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
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index d1c9ffb1..fbf8bc77 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -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."
--
GitLab