From 94f581284a7876d3215b02aa63004564600e4f8c Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Thu, 26 Oct 2023 17:51:36 +0200
Subject: [PATCH] Introduce FORTRAN/C++ correspondence logs
---
src/include/sph_subs.h | 133 +++++++++++++++++++++++++++++++++++++++++
src/sphere/sph.f | 1 +
src/sphere/sphere.cpp | 7 ++-
3 files changed, 138 insertions(+), 3 deletions(-)
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index 411a8ec0..081b13ac 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -49,6 +49,139 @@ struct C2 {
double vsz[2];
};
+/*! \brief Conjugate of a double precision complex number
+ *
+ * \param z: `std::complex\<double\>` The input complex number.
+ * \return result: `std::complex\<double\>` The conjugate of the input number.
+ */
+std::complex<double> dconjg(std::complex<double> z) {
+ double zreal = z.real();
+ double zimag = z.imag();
+ return std::complex<double>(zreal, -zimag);
+}
+
+/*! \brief C++ porting of CG1
+ *
+ * \param li: `int`
+ * \param mu: `int`
+ * \param l: `int`
+ * \param m: `int`
+ * \return result: `double`
+ */
+double cg1(int lmpml, int mu, int l, int m) {
+ double result = 0.0;
+ double xd, xn;
+ if (lmpml == -1) { // Interpreted as GOTO 30
+ xd = 2.0 * l * (2 * l - 1);
+ if (mu == -1) {
+ xn = 1.0 * (l - 1 - m) * (l - m);
+ } else if (mu == 0) {
+ xn = 2.0 * (l - m) * (l + m);
+ } else if (mu == 1) {
+ xn = 1.0 * (l - 1 + m) * (l + m);
+ } else {
+ throw 111; // Need an exception for unpredicted indices.
+ }
+ result = sqrt(xn / xd);
+ } else if (lmpml == 0) { // Interpreted as GOTO 5
+ bool goto10 = (m != 0) || (mu != 0);
+ if (!goto10) {
+ result = 0.0;
+ return result;
+ }
+ if (mu != 0) {
+ xd = 2.0 * l * (l + 1);
+ if (mu == -1) {
+ xn = 1.0 * (l - m) * (l + m + 1);
+ result = -sqrt(xn / xd);
+ } else if (mu == 1) { // mu > 0
+ xn = 1.0 * (l + m) * (l - m + 1);
+ result = sqrt(xn / xd);
+ } else {
+ throw 111; // Need an exception for unpredicted indices.
+ }
+ } else { // mu = 0
+ xd = 1.0 * l * (l + 1);
+ xn = -1.0 * m;
+ result = xn / sqrt(xd);
+ }
+ } else if (lmpml == 1) { // Interpreted as GOTO 60
+ xd = 2.0 * (l * 2 + 3) * (l + 1);
+ if (mu == -1) {
+ xn = 1.0 * (l + 1 + m) * (l + 2 + m);
+ result = sqrt(xn / xd);
+ } else if (mu == 0) {
+ xn = 2.0 * (l + 1 - m) * (l + 1 + m);
+ result = -sqrt(xn / xd);
+ } else if (mu == 1) {
+ xn = 1.0 * (l + 1 - m) * (l + 2 - m);
+ result = sqrt(xn / xd);
+ } else { // mu was not recognized.
+ throw 111; // Need an exception for unpredicted indices.
+ }
+ } else { // lmpml was not recognized
+ throw 111; // Need an exception for unpredicted indices.
+ }
+ return result;
+}
+
+/*! \brief C++ porting of APS
+ *
+ * \param zpv: `double ****`
+ * \param li: `int`
+ * \param nsph: `int`
+ * \param c1: `C1 *`
+ * \param sqk: `double`
+ * \param gaps: `double *`
+ */
+void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
+ std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
+ std::complex<double> summ, sume, suem, suee, sum;
+ double half_pi = acos(0.0);
+ double cofs = half_pi * 2.0 / sqk;
+ for (int i40 = 1; i40 <= nsph; i40++) {
+ int iogi = c1->iog[i40 - 1];
+ if (iogi >= i40) {
+ sum = cc0;
+ for (int l30 = 1; l30 <= li; l30++) {
+ int ltpo = l30 + l30 + 1;
+ for (int ilmp = 1; ilmp <= 3; ilmp++) {
+ bool goto30 = (l30 == 1 && ilmp == 1) || (l30 == li && ilmp == 3);
+ if (!goto30) {
+ int lmpml = ilmp - 2;
+ int lmp = l30 + lmpml;
+ double cofl = sqrt(1.0 * (ltpo * (lmp + lmp + 1)));
+ summ = zpv[l30 - 1][ilmp - 1][0][0] /
+ (
+ dconjg(c1->rmi[l30 - 1][i40 - 1]) *
+ c1->rmi[lmp - 1][i40 - 1]
+ );
+ sume = zpv[l30 - 1][ilmp - 1][0][1] /
+ (
+ dconjg(c1->rmi[l30 - 1][i40 - 1]) *
+ c1->rei[lmp - 1][i40 - 1]
+ );
+ suem = zpv[l30 - 1][ilmp - 1][1][0] /
+ (
+ dconjg(c1->rei[l30 - 1][i40 - 1]) *
+ c1->rmi[lmp - 1][i40 - 1]
+ );
+ suee = zpv[l30 - 1][ilmp - 1][1][1] /
+ (
+ dconjg(c1->rei[l30 - 1][i40 - 1]) *
+ c1->rei[lmp - 1][i40 - 1]
+ );
+ sum += (cg1(lmpml, 0, l30, -1) * (summ - sume - suem + suee) +
+ cg1(lmpml, 0, l30, 1) * (summ + sume + suem + suee)) * cofl;
+ }
+ }
+ }
+ }
+ gaps[i40 - 1] = sum.real() * cofs;
+ printf("DEBUG: gaps[%d] = %lE\n", i40, gaps[i40 - 1]);
+ }
+}
+
/*! \brief C++ porting of DIEL
*
* \param npntmo: `int`
diff --git a/src/sphere/sph.f b/src/sphere/sph.f
index 69d8e3b3..77a6df97 100644
--- a/src/sphere/sph.f
+++ b/src/sphere/sph.f
@@ -645,6 +645,7 @@ CCC 1RMI(LI,NSPH),REI(LI,NSPH),GAPS(NSPH)
1CG1(LMPML,0,L,1)*(SUMM+SUME+SUEM+SUEE))*COFL
30 CONTINUE
GAPS(I)=DREAL(SUM)*COFS
+ PRINT *,"DEBUG: GAPS(",I,") =",GAPS(I)
40 CONTINUE
RETURN
END
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 3f87c90d..e422603d 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -9,9 +9,8 @@ using namespace std;
//! \brief C++ implementation of SPH
void sphere() {
- //complex<double> dc0[5];
- //complex<double> dc0m[2][4];
complex<double> arg, s0, tfsas;
+ double gaps[2];
printf("INFO: making legacy configuration ...\n");
ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/sphere/DEDFB");
conf->write_formatted("c_OEDFB");
@@ -204,7 +203,6 @@ void sphere() {
int iogi = sconf->iog_vec[i132 - 1];
if (iogi != i132) {
for (int l123 = 1; l123 <= gconf->l_max; l123++) {
- // QUESTION: aren't RMI and REI still empty?
c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];
c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];
}
@@ -265,6 +263,9 @@ void sphere() {
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());
+ double sqk = vk * vk * sconf->exdc;
+ aps(zpv, gconf->l_max, nsph, c1, sqk, gaps);
+ // Would call RABAS
} //jxi488
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
--
GitLab