From 685d3e1ba4d3b80e30389434f6790af65a15a647 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Tue, 28 Nov 2023 12:48:21 +0100
Subject: [PATCH] Extend conversion of CLU to C++
---
src/cluster/cluster.cpp | 81 ++++++++++++++++++++++++++++++++++++++---
src/include/clu_subs.h | 59 +++++++++++++++++++++++++++++-
2 files changed, 133 insertions(+), 7 deletions(-)
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index ad983afc..6e6311aa 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -19,13 +19,13 @@ using namespace std;
//! \brief C++ implementation of CLU
void cluster() {
printf("INFO: making legacy configuration ...\n");
- ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("DEDFB_clu");
+ ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/cluster/DEDFB");
conf->write_formatted("c_OEDFB_clu");
conf->write_binary("c_TEDF_clu");
delete conf;
printf("INFO: reading binary configuration ...\n");
ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF_clu");
- GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("DCLU");
+ GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/cluster/DCLU");
if (sconf->number_of_spheres == gconf->number_of_spheres) {
// Shortcuts to variables stored in configuration objects
int nsph = gconf->number_of_spheres;
@@ -93,6 +93,19 @@ void cluster() {
for (int ai = 0; ai < mxndm; ai++) am[ai] = new complex<double>[mxndm];
const int ndi = c4->nsph * c4->nlim;
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
+ double *gaps = new double[nsph]();
+ double **tqse, **tqss;
+ complex<double> **tqspe, **tqsps;
+ tqse = new double*[2];
+ tqspe = new complex<double>*[2];
+ tqss = new double*[2];
+ tqsps = new complex<double>*[2];
+ for (int ti = 0; ti < 2; ti++) {
+ tqse[ti] = new double[nsph]();
+ tqspe[ti] = new complex<double>[nsph]();
+ tqss[ti] = new double[nsph]();
+ tqsps[ti] = new complex<double>[nsph]();
+ }
// 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",
@@ -155,9 +168,7 @@ void cluster() {
for (int zj = 0; zj < 3; zj++) {
zpv[zi][zj] = new double*[2];
for (int zk = 0; zk < 2; zk++) {
- zpv[zi][zj][zk] = new double[2];
- zpv[zi][zj][zk][0] = 0.0;
- zpv[zi][zj][zk][1] = 0.0;
+ zpv[zi][zj][zk] = new double[2]();
}
}
}
@@ -271,7 +282,54 @@ void cluster() {
}
}
// label 156: continue from here
-
+ if (inpol == 0) {
+ fprintf(output, " LIN\n");
+ } else { // label 158
+ fprintf(output, " CIRC\n");
+ }
+ // label 160
+ double cs0 = 0.25 * vk * vk * vk / acos(0.0);
+ double csch;
+ scr0(vk, exri, c1, c1ao, c3, c4);
+ printf("DEBUG: after SCR0 TFSAS = (%lE, %lE)\n", c3->tfsas.real(), c3->tfsas.imag());
+ double sqk = vk * vk * sconf->exdc;
+ aps(zpv, c4->li, nsph, c1, sqk, gaps);
+ rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps);
+ if (c4->li != c4->le) fprintf(output, " SPHERES; LMX=LI\n");
+ for (int i170 = 1; i170 <= nsph; i170++) {
+ if (c1->iog[i170 - 1] >= i170) {
+ int i = i170 - 1;
+ double albeds = c1->sscs[i] / c1->sexs[i];
+ c1->sqscs[i] *= sqsfi;
+ c1->sqabs[i] *= sqsfi;
+ c1->sqexs[i] *= sqsfi;
+ fprintf(output, " SPHERE %2d\n", i170);
+ 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());
+ }
+ // label 164
+ fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+ fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", c1->sscs[i], c1->sabs[i], c1->sexs[i], albeds);
+ fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
+ fprintf(output, " %14.7lE%15.7lE%15.7lE\n", c1->sqscs[i], c1->sqabs[i], c1->sqexs[i]);
+ fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i].real(), c1->fsas[i].imag());
+ csch = 2.0 * vk * sqsfi / c1->gcsv[i];
+ std::complex<double> s0 = c1->fsas[i] * exri;
+ double qschu = s0.imag() * csch;
+ double pschu = s0.real() * csch;
+ double s0mag = sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag())) * cs0;
+ fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
+ double rapr = c1->sexs[i] - gaps[i];
+ double cosav = gaps[i] / c1->sscs[i];
+ fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
+ fprintf(output, " IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[0][i], tqss[0][i]);
+ fprintf(output, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[1][i], tqss[1][i]);
+ }
+ } // i170 loop
+ fprintf(output, " FSAT=%15.7lE%15.7lE\n", c3->tfsas.real(), c3->tfsas.imag());
+ printf("INFO: done jxi488 iteration.\n");
} // jxi488 loop
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
@@ -296,6 +354,17 @@ void cluster() {
delete[] zpv;
for (int ai = mxndm - 1; ai > -1; ai--) delete[] am[ai];
delete[] am;
+ delete[] gaps;
+ for (int ti = 1; ti > -1; ti--) {
+ delete[] tqse[ti];
+ delete[] tqss[ti];
+ delete[] tqspe[ti];
+ delete[] tqsps[ti];
+ }
+ delete[] tqse;
+ delete[] tqss;
+ delete[] tqspe;
+ delete[] tqsps;
} else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException(
"Inconsistent geometry and scatterer configurations."
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index c10aecb1..6e91baa4 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -21,14 +21,19 @@
#include <complex>
// >>> DECLARATION OF SPH_SUBS <<<
+extern void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps);
extern std::complex<double> dconjg(std::complex<double> value);
extern void dme(
int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg
);
-extern void sphar(double cth, double sth, double cph, double sph, int lm, std::complex<double> *ylm);
+extern void rabas(
+ int inpol, int li, int nsph, C1 *c1, double **tqse, std::complex<double> **tqspe,
+ double **tqss, std::complex<double> **tqsps
+);
extern void rbf(int n, double x, int &nm, double sj[]);
extern void rnf(int n, double x, int &nm, double sy[]);
+extern void sphar(double cth, double sth, double cph, double sph, int lm, std::complex<double> *ylm);
extern void thdps(int lm, double ****zpv);
// >>> END OF SPH_SUBS DECLARATION <<<
void logmat(std::complex<double> **mat, int rows, int cols);
@@ -910,6 +915,58 @@ void r3j000(int j2, int j3, C6 *c6) {
}
}
+/*! \brief C++ porting of SCR0
+ *
+ * \param vk: `double` QUESTION: definition?
+ * \param exri: `double` External medium refractive index. QUESTION: correct?
+ * \param c1: `C1 *` Pointer to a C1 instance.
+ * \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.
+ * \param c3: `C3 *` Pointer to a C3 instance.
+ * \param c4: `C4 *` Pointer to a C4 structure.
+ */
+void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
+ const std::complex<double> cc0(0.0, 0.0);
+ double exdc = exri * exri;
+ double ccs = 4.0 * acos(0.0) / (vk * vk);
+ 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;
+ c3->tfsas = cc0;
+ for (int i14 = 1; i14 <= c4->nsph; i14++) {
+ int iogi = c1->iog[i14 - 1];
+ if (iogi >= i14) {
+ double sums = 0.0;
+ sum21 = cc0;
+ for (int l10 = 1; l10 <= c4->li; l10++) {
+ double fl = 1.0 * (l10 + l10 + 1);
+ rm = 1.0 / c1->rmi[l10 - 1][i14 - 1];
+ re = 1.0 / c1->rei[l10 - 1][i14 - 1];
+ double rvalue = (dconjg(rm) * rm + dconjg(re) * re).real() * fl;
+ sums += rvalue;
+ sum21 += ((rm + re) * fl);
+ } // l10 loop
+ sum21 *= -1.0;
+ double scasec = cccs * sums;
+ double extsec = -cccs * sum21.real();
+ double abssec = extsec - scasec;
+ c1->sscs[i14 - 1] = scasec;
+ c1->sexs[i14 - 1] = extsec;
+ c1->sabs[i14 - 1] = abssec;
+ double gcss = c1->gcsv[i14 - 1];
+ c1->sqscs[i14 - 1] = scasec / gcss;
+ c1->sqexs[i14 - 1] = extsec / gcss;
+ c1->sqabs[i14 - 1] = abssec / gcss;
+ c1->fsas[i14 - 1] = sum21 * csam;
+ }
+ // label 12
+ scs += c1->sscs[iogi - 1];
+ ecs += c1->sexs[iogi - 1];
+ acs += c1->sabs[iogi - 1];
+ c3->tfsas += c1->fsas[iogi - 1];
+ } // i14 loop
+}
+
/*! \brief C++ porting of STR
*
* \param rcf: `double **` Matrix of double.
--
GitLab