diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 6e6311aa50266c4c62ba56fbec7a140404adaeee..7b7e27c5629f0792363ab7cf3b733b462a4c83d8 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -94,19 +94,49 @@ void cluster() {
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;
+ double **tqse, **tqss, **tqce, **tqcs;
+ complex<double> **tqspe, **tqsps, **tqcpe, **tqcps;
tqse = new double*[2];
tqspe = new complex<double>*[2];
tqss = new double*[2];
tqsps = new complex<double>*[2];
+ tqce = new double*[2];
+ tqcpe = new complex<double>*[2];
+ tqcs = new double*[2];
+ tqcps = 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]();
+ tqce[ti] = new double[3]();
+ tqcpe[ti] = new complex<double>[3]();
+ tqcs[ti] = new double[3]();
+ tqcps[ti] = new complex<double>[3]();
}
- // End of global variables for CLU
+ complex<double> **gapp, **gappm;
+ double **gap, **gapm;
+ gapp = new complex<double>*[3];
+ gappm = new complex<double>*[3];
+ gap = new double*[3];
+ gapm = new double*[3];
+ for (int gi = 0; gi < 3; gi++) {
+ gapp[gi] = new complex<double>[2]();
+ gappm[gi] = new complex<double>[2]();
+ gap[gi] = new double[2]();
+ gapm[gi] = new double[2]();
+ }
+ 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 *unmp = new double[3]();
+ double *unsmp = new double[3]();
+ double *upmp = new double[3]();
+ double *upsmp = new double[3]();
+ // 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",
nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam
@@ -289,7 +319,8 @@ void cluster() {
}
// label 160
double cs0 = 0.25 * vk * vk * vk / acos(0.0);
- double csch;
+ double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
+ std::complex<double> s0(0.0, 0.0);
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;
@@ -316,10 +347,10 @@ void cluster() {
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;
+ s0 = c1->fsas[i] * exri;
+ qschu = s0.imag() * csch;
+ pschu = s0.real() * csch;
+ 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];
@@ -329,6 +360,48 @@ void cluster() {
}
} // i170 loop
fprintf(output, " FSAT=%15.7lE%15.7lE\n", c3->tfsas.real(), c3->tfsas.imag());
+ csch = 2.0 * vk * sqsfi / c3->gcs;
+ s0 = c3->tfsas * exri;
+ qschu = s0.imag() * csch;
+ pschu = s0.real() * csch;
+ 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);
+ tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
+ pcrsm0(vk, exri, inpol, c1, c1ao, c4);
+ apcra(zpv, c4->le, c1ao->am0m, inpol, sqk, gapm, gappm);
+ th = th1;
+ for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP portable?
+ ph = ph1;
+ double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
+ // argi[NSPEF], with NSPEF=1 if IDOT=0, else NSPEF=NSPH
+ double *argi;
+ for (int jph484 = 1; jph484 <= nph; jph484++) {
+ int jw = 0;
+ if (nk != 1 || jxi488 <= 1) {
+ upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
+ if (isam >= 0) {
+ argi = new double[1];
+ wmamp(
+ 0, cost, sint, cosp, sinp, inpol, c4->le, 0,
+ nsph, argi, u, upmp, unmp, c1
+ );
+ // label 182
+ apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
+ raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
+ jw = 1;
+ }
+ } else { // label 180, NK == 1 AND JXI488 == 1
+ if (isam >= 0) {
+ // label 182
+ apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
+ raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
+ jw = 1;
+ }
+ }
+ // label 184
+ double thsl = ths1;
+ } // jph484 loop
+ } // jth486 loop
printf("INFO: done jxi488 iteration.\n");
} // jxi488 loop
tppoan.close();
@@ -360,11 +433,39 @@ void cluster() {
delete[] tqss[ti];
delete[] tqspe[ti];
delete[] tqsps[ti];
+ delete[] tqce[ti];
+ delete[] tqcpe[ti];
+ delete[] tqcs[ti];
+ delete[] tqcps[ti];
}
delete[] tqse;
delete[] tqss;
delete[] tqspe;
delete[] tqsps;
+ delete[] tqce;
+ delete[] tqcpe;
+ delete[] tqcs;
+ delete[] tqcps;
+ for (int gi = 2; gi > -1; gi--) {
+ delete[] gapp[gi];
+ delete[] gappm[gi];
+ delete[] gap[gi];
+ delete[] gapm[gi];
+ }
+ delete[] gapp;
+ delete[] gappm;
+ delete[] gap;
+ delete[] gapm;
+ delete[] u;
+ delete[] us;
+ delete[] un;
+ delete[] uns;
+ delete[] up;
+ delete[] ups;
+ delete[] unmp;
+ delete[] unsmp;
+ delete[] upmp;
+ delete[] upsmp;
} else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException(
"Inconsistent geometry and scatterer configurations."
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 01571dfb8af7e0a5b9c0dec92c86719de498dc8b..81f0cf1b64483ec46eed7a275c7b681fcf0de3d2 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -238,6 +238,10 @@ public:
//! \brief QUESTION: definition?
std::complex<double> *scscp;
//! \brief QUESTION: definition?
+ double *ecscm;
+ //! \brief QUESTION: definition?
+ double *scscm;
+ //! \brief QUESTION: definition?
std::complex<double> *ecscp;
//! \brief QUESTION: definition?
std::complex<double> *scscpm;
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 6e91baa4b13e47e2e77b420e4fd6f3ca6af9af60..f49166b69ab84903f9e30b6cfc073bb413684fef 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -23,6 +23,7 @@
// >>> 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 double cg1(int lmpml, int mu, int l, int m);
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
@@ -35,6 +36,15 @@ 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);
+extern void upvmp(
+ double thd, double phd, int icspnv, double &cost, double &sint,
+ double &cosp, double &sinp, double *u, double *up, double *un
+);
+extern 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
+);
// >>> END OF SPH_SUBS DECLARATION <<<
void logmat(std::complex<double> **mat, int rows, int cols);
@@ -235,7 +245,6 @@ void r3jjr(int j2, int j3, int m2, int m3, C6 *c6) {
}
}
-
/*! \brief C++ porting of GHIT
*
* \param ihi: `int`
@@ -463,6 +472,365 @@ std::complex<double> ghit(
return result;
}
+/*! \brief C++ porting of APC
+ *
+ * \param zpv: `double ****`
+ * \param le: `int`
+ * \param am0m: Matrix of complex.
+ * \param w: Matrix of complex.
+ * \param sqk: `double`
+ * \param gapr: `double **`
+ * \param gapp: Matrix of complex.
+ */
+void apc(
+ double ****zpv, int le, std::complex<double> **am0m, std::complex<double> **w,
+ double sqk, double **gapr, std::complex<double> **gapp
+) {
+ std::complex<double> **ac, **gap;
+ const std::complex<double> cc0(0.0, 0.0);
+ const std::complex<double> uim(0.0, 1.0);
+ std::complex<double> uimmp, summ, sume, suem, suee, summp, sumep;
+ std::complex<double> suemp, sueep;
+ double cof = 1.0 / sqk;
+ double cimu = cof / sqrt(2.0);
+ int nlem = le * (le + 2);
+ const int nlemt = nlem + nlem;
+ ac = new std::complex<double>*[nlemt];
+ gap = new std::complex<double>*[3];
+ for (int ai = 0; ai < nlemt; ai++) ac[ai] = new std::complex<double>[2]();
+ for (int gi = 0; gi < 3; gi++) gap[gi] = new std::complex<double>[2]();
+ for (int j45 = 1; j45 <= nlemt; j45++) {
+ int j = j45 - 1;
+ ac[j][0] = cc0;
+ ac[j][1] = cc0;
+ for (int i45 = 1; i45 <= nlemt; i45++) {
+ int i = i45 - 1;
+ ac[j][0] += (am0m[j][i] * w[i][0]);
+ ac[j][1] += (am0m[j][i] * w[i][1]);
+ } //i45 loop
+ } //j45 loop
+ for (int imu90 = 1; imu90 <=3; imu90++) {
+ int mu = imu90 - 2;
+ gap[imu90 - 1][0] = cc0;
+ gap[imu90 - 1][1] = cc0;
+ gapp[imu90 - 1][0] = cc0;
+ gapp[imu90 - 1][1] = cc0;
+ for (int l80 =1; l80 <= le; l80++) {
+ int lpo = l80 + 1;
+ int ltpo = lpo + l80;
+ int imm = l80 * lpo;
+ for (int ilmp = 1; ilmp <= 3; ilmp++) {
+ if ((l80 == 1 && ilmp == 1) || (l80 == le && ilmp == 3)) continue; // ilmp loop
+ int lmpml = ilmp - 2;
+ int lmp = l80 + lmpml;
+ uimmp = (-1.0 * lmpml) * uim;
+ int impmmmp = lmp * (lmp + 1);
+ for (int im70 = 1; im70 <= ltpo; im70++) {
+ int m = im70 - lpo;
+ int mmp = m - mu;
+ int abs_mmp = (mmp > 0) ? mmp : -mmp;
+ if (abs_mmp <= lmp) {
+ int i = imm + m;
+ int ie = i + nlem;
+ int imp = impmmmp + mmp;
+ int impe = imp + nlem;
+ double cgc = cg1(lmpml, mu, l80, m);
+ int jpo = 2;
+ for (int ipo = 1; ipo <= 2; ipo++) {
+ if (ipo == 2) jpo = 1;
+ //printf("DEBUG: i=%d, ipo=%d, imp=%d\n", i, ipo, imp);
+ //fflush(stdout);
+ summ = dconjg(ac[i - 1][ipo - 1]) * ac[imp - 1][ipo - 1];
+ sume = dconjg(ac[i - 1][ipo - 1]) * ac[impe - 1][ipo - 1];
+ suem = dconjg(ac[ie - 1][ipo - 1]) * ac[imp - 1][ipo - 1];
+ suee = dconjg(ac[ie - 1][ipo - 1]) * ac[impe - 1][ipo - 1];
+ summp = dconjg(ac[i - 1][jpo - 1]) * ac[imp - 1][ipo - 1];
+ sumep = dconjg(ac[i - 1][jpo - 1]) * ac[impe - 1][ipo - 1];
+ suemp = dconjg(ac[ie - 1][jpo - 1]) * ac[imp - 1][ipo - 1];
+ sueep = dconjg(ac[ie - 1][jpo - 1]) * ac[impe - 1][ipo - 1];
+ if (lmpml != 0) {
+ summ *= uimmp;
+ sume *= uimmp;
+ suem *= uimmp;
+ suee *= uimmp;
+ summp *= uimmp;
+ sumep *= uimmp;
+ suemp *= uimmp;
+ sueep *= uimmp;
+ }
+ // label 55
+ gap[imu90 - 1][ipo - 1] += (
+ (
+ summ * zpv[l80 - 1][ilmp - 1][0][0]
+ + sume * zpv[l80 - 1][ilmp - 1][0][1]
+ + suem * zpv[l80 - 1][ilmp - 1][1][0]
+ + suee * zpv[l80 - 1][ilmp - 1][1][1]
+ ) * cgc
+ );
+ gapp[imu90 - 1][ipo - 1] += (
+ (
+ summp * zpv[l80 - 1][ilmp - 1][0][0]
+ + sumep * zpv[l80 - 1][ilmp - 1][0][1]
+ + suemp * zpv[l80 - 1][ilmp - 1][1][0]
+ + sueep * zpv[l80 - 1][ilmp - 1][1][1]
+ ) * cgc
+ );
+ } // ipo loop
+ } // ends im70 loop
+ } // im70 loop
+ } // ilmp loop
+ } // l80 loop
+ } // imu90 loop
+ for (int ipo95 = 1; ipo95 <= 2; ipo95++) {
+ sume = gap[0][ipo95 - 1] * cimu;
+ suee = gap[1][ipo95 - 1] * cof;
+ suem = gap[2][ipo95 - 1] * cimu;
+ gapr[0][ipo95 - 1] = (sume - suem).real();
+ gapr[1][ipo95 - 1] = ((sume + suem) * uim).real();
+ gapr[2][ipo95 - 1] = suee.real();
+ sumep = gapp[0][ipo95 - 1] * cimu;
+ sueep = gapp[1][ipo95 - 1] * cof;
+ suemp = gapp[2][ipo95 - 1] * cimu;
+ gapp[0][ipo95 - 1] = sumep - suemp;
+ gapp[1][ipo95 - 1] = (sumep + suemp) * uim;
+ gapp[2][ipo95 - 1] = sueep;
+ } // ipo95 loop
+ // Clean memory
+ for (int ai = nlemt - 1; ai > -1; ai--) delete[] ac[ai];
+ for (int gi = 2; gi > -1; gi--) delete[] gap[gi];
+ delete[] ac;
+ delete[] gap;
+}
+
+/*! \brief C++ porting of APCRA
+ *
+ * \param zpv: `double ****`
+ * \param le: `int`
+ * \param am0m: Matrix of complex.
+ * \param inpol: `int` Polarization type.
+ * \param sqk: `double`
+ * \param gaprm: `double **`
+ * \param gappm: Matrix of complex.
+ */
+void apcra(
+ double ****zpv, const int le, std::complex<double> **am0m, int inpol, double sqk,
+ double **gaprm, std::complex<double> **gappm
+) {
+ const std::complex<double> cc0(0.0, 0.0);
+ const std::complex<double> uim(0.0, 1.0);
+ std::complex<double> uimtl, uimtls, ca11, ca12, ca21, ca22;
+ std::complex<double> a11, a12, a21, a22, sum1, sum2, fc;
+ double ****svw = new double***[le];
+ std::complex<double> ****svs = new std::complex<double>***[le];
+ for (int i = 0; i < le; i++) {
+ svw[i] = new double**[3];
+ svs[i] = new std::complex<double>**[3];
+ for (int j = 0; j < 3; j++) {
+ svw[i][j] = new double*[2];
+ svs[i][j] = new std::complex<double>*[2];
+ for (int k = 0; k < 2; k++) {
+ svw[i][j][k] = new double[2]();
+ svs[i][j][k] = new std::complex<double>[2]();
+ }
+ }
+ }
+ int nlem = le * (le + 2);
+ for (int l28 = 1; l28 <= le; l28++) {
+ int lpo = l28 + 1;
+ int ltpo = lpo + l28;
+ double fl = sqrt(1.0 * ltpo);
+ for (int ilmp = 1; ilmp <= 3; ilmp++) {
+ if ((l28 == 1 && ilmp == 1) || (l28 == le && ilmp == 3)) continue; // ilmp loop
+ int lmpml = ilmp - 2;
+ int lmp = l28 + lmpml;
+ double flmp = sqrt(1.0 * (lmp + lmp + 1));
+ double fllmp = flmp / fl;
+ double cgmmo = fllmp * cg1(lmpml, 0, l28, 1);
+ double cgmpo = fllmp * cg1(lmpml, 0, l28, -1);
+ if (inpol == 0) {
+ double cgs = cgmpo + cgmmo;
+ double cgd = cgmpo - cgmmo;
+ svw[l28 - 1][ilmp - 1][0][0] = cgs;
+ svw[l28 - 1][ilmp - 1][0][1] = cgd;
+ svw[l28 - 1][ilmp - 1][1][0] = cgd;
+ svw[l28 - 1][ilmp - 1][1][1] = cgs;
+ } else { // label 22
+ svw[l28 - 1][ilmp - 1][0][0] = cgmpo;
+ svw[l28 - 1][ilmp - 1][1][0] = cgmpo;
+ svw[l28 - 1][ilmp - 1][0][1] = -cgmmo;
+ svw[l28 - 1][ilmp - 1][1][1] = cgmmo;
+ }
+ // label 26
+ } // ilmp loop
+ } // l28 loop
+ for (int l30 = 1; l30 <= le; l30++) { // This can be omitted, since svs was initialized at 0
+ for (int ilmp = 1; ilmp <= 3; ilmp++) {
+ for (int ipa = 1; ipa <= 2; ipa++) {
+ for (int ipamp = 1; ipamp <= 2; ipamp++) {
+ svs[l30 - 1][ilmp - 1][ipa - 1][ipamp - 1] = cc0;
+ }
+ } // ipa loop
+ } // ilmp loop
+ } // l30 loop
+ for (int l58 = 1; l58 <= le; l58 ++) {
+ int lpo = l58 + 1;
+ int ltpo = l58 + lpo;
+ int imm = l58 * lpo;
+ for (int ilmp = 1; ilmp <= 3; ilmp++) {
+ if ((l58 == 1 && ilmp == 1) || (l58 == le && ilmp == 3)) continue; // ilmp loop
+ int lmpml = ilmp - 2;
+ int lmp = l58 + lmpml;
+ int impmm = lmp * (lmp + 1);
+ uimtl = uim * (1.0 * lmpml);
+ if (lmpml == 0) uimtl = std::complex<double>(1.0, 0.0);
+ for (int im54 = 1; im54 <= ltpo; im54++) {
+ int m = im54 - lpo;
+ int i = imm + m;
+ int ie = i + nlem;
+ for (int imu52 = 1; imu52 <= 3; imu52++) {
+ int mu = imu52 - 2;
+ int mmp = m - mu;
+ int abs_mmp = (mmp > 0) ? mmp : -mmp;
+ if (abs_mmp <= lmp) {
+ int imp = impmm + mmp;
+ int impe = imp + nlem;
+ double cgc = cg1(lmpml, -mu, l58, -m);
+ for (int ls = 1; ls <= le; ls++) {
+ int lspo = ls + 1;
+ int lstpo = ls + lspo;
+ int ismm = ls * lspo;
+ for (int ilsmp = 1; ilsmp <= 3; ilsmp++) {
+ if ((ls == 1 && ilsmp == 1) || (ls == le && ilsmp == 3)) continue; // ilsmp loop
+ int lsmpml = ilsmp - 2;
+ int lsmp = ls + lsmpml;
+ int ismpmm = lsmp * (lsmp + 1);
+ uimtls = -uim * (1.0 * lsmpml);
+ if (lsmpml == 0) uimtls = std::complex<double>(1.0, 0.0);
+ for (int ims = 1; ims <= lstpo; ims++) {
+ int ms = ims - lspo;
+ int msmp = ms - mu;
+ int abs_msmp = (msmp > 0) ? msmp : -msmp;
+ if (abs_msmp <= lsmp) {
+ int is = ismm + ms;
+ int ise = is + nlem;
+ int ismp = ismpmm + msmp;
+ int ismpe = ismp + nlem;
+ double cgcs = cg1(lsmpml, mu, ls, ms);
+ fc = (uimtl * uimtls) * (cgc * cgcs);
+ ca11 = dconjg(am0m[is - 1][i - 1]);
+ ca12 = dconjg(am0m[is - 1][ie - 1]);
+ ca21 = dconjg(am0m[ise - 1][i - 1]);
+ ca22 = dconjg(am0m[ise - 1][ie - 1]);
+ a11 = am0m[ismp - 1][imp - 1];
+ a12 = am0m[ismp - 1][impe - 1];
+ a21 = am0m[ismpe - 1][imp - 1];
+ a22 = am0m[ismpe - 1][impe - 1];
+ double z11 = zpv[ls - 1][ilsmp - 1][0][0];
+ double z12 = zpv[ls - 1][ilsmp - 1][0][1];
+ double z21 = zpv[ls - 1][ilsmp - 1][1][0];
+ double z22 = zpv[ls - 1][ilsmp - 1][1][1];
+ svs[l58 - 1][ilmp - 1][0][0] += ((ca11 * a11 * z11
+ + ca11 * a21 * z12
+ + ca21 * a11 * z21
+ + ca21 * a21 * z22) * fc);
+ svs[l58 - 1][ilmp - 1][0][1] += ((ca11 * a12 * z11
+ + ca11 * a22 * z12
+ + ca21 * a12 * z21
+ + ca21 * a22 * z22) * fc);
+ svs[l58 - 1][ilmp - 1][1][0] += ((ca12 * a11 * z11
+ + ca12 * a21 * z12
+ + ca22 * a11 * z21
+ + ca22 * a21 * z22) * fc);
+ svs[l58 - 1][ilmp - 1][1][1] += ((ca12 * a12 * z11
+ + ca12 * a22 * z12
+ + ca22 * a12 * z21
+ + ca22 * a22 * z22) * fc);
+ } // ends ims loop
+ } // ims loop
+ } // ilsmp loop
+ } // ls loop
+ } // ends imu52 loop
+ } // imu52 loop
+ } // im54 loop
+ } // ilmp loop
+ } // l58 loop
+ sum1 = cc0;
+ sum2 = cc0;
+ for (int l68 = 1; l68 <= le; l68++) {
+ int lpo = l68 + 1;
+ int ltpo = l68 + lpo;
+ int imm = l68 * lpo;
+ for (int ilmp = 1; ilmp <= 3; ilmp++) {
+ if ((l68 == 1 && ilmp == 1) || (l68 == le && ilmp == 3)) continue; // ilmp loop
+ if (inpol == 0) {
+ sum1 += (
+ svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][0][0]
+ + svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][0][1]
+ + svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][1][0]
+ + svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][1][1]
+ );
+ sum2 += (
+ svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][0][0]
+ + svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][0][1]
+ + svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][1][0]
+ + svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][1][1]
+ );
+ } else { // label 62
+ sum1 += (
+ svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][0][0]
+ + svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][0][1]
+ + svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][1][0]
+ + svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][1][1]
+ );
+ sum2 += (
+ svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][0][0]
+ + svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][0][1]
+ + svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][1][0]
+ + svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][1][1]
+ );
+ } // label 66, ends ilmp loop
+ } // ilmp loop
+ } // l68 loop
+ const double half_pi = acos(0.0);
+ double cofs = half_pi * 2.0 / sqk;
+ gaprm[0][0] = 0.0;
+ gaprm[0][1] = 0.0;
+ gaprm[1][0] = 0.0;
+ gaprm[1][1] = 0.0;
+ gappm[0][0] = cc0;
+ gappm[0][1] = cc0;
+ gappm[1][0] = cc0;
+ gappm[1][1] = cc0;
+ if (inpol == 0) {
+ sum1 *= cofs;
+ sum2 *= cofs;
+ gaprm[2][0] = sum1.real();
+ gaprm[2][1] = sum1.real();
+ gappm[2][0] = sum2 * uim;
+ gappm[2][1] = -gappm[2][0];
+ } else { // label 72
+ gaprm[2][0] = sum1.real() * cofs;
+ gaprm[2][1] = sum2.real() * cofs;
+ gappm[2][0] = cc0;
+ gappm[2][1] = cc0;
+ }
+ // Clean memory
+ for (int i = le - 1; i > -1; i--) {
+ for (int j = 2; j > -1; j--) {
+ for (int k = 1; k > -1; k--) {
+ delete[] svw[i][j][k];
+ delete[] svs[i][j][k];
+ }
+ delete[] svw[i][j];
+ delete[] svs[i][j];
+ }
+ delete[] svw[i];
+ delete[] svs[i];
+ }
+ delete[] svw;
+ delete[] svs;
+}
+
/*! \brief C++ porting of CMS
*
* \param am: Matrix of complex.
@@ -746,6 +1114,81 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
delete[] v;
}
+/*! \brief C++ porting of PCRSM0
+ *
+ * \param vk: `double`
+ * \param exri: `double`
+ * \param inpol: `int`
+ * \param c1: `C1 *`
+ * \param c1ao: `C1_AddOns *`
+ * \param c4: `C4 *`
+ */
+void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
+ const std::complex<double> cc0(0.0, 0.0);
+ const std::complex<double> uim(0.0, 1.0);
+ std::complex<double> sum1, sum2, sum3, sum4, sumpd;
+ std::complex<double> sums1, sums2, sums3, sums4, csam;
+ double exdc = exri * exri;
+ double ccs = 4.0 * acos(0.0) / (vk * vk);
+ double cccs = ccs / exdc;
+ csam = -(ccs / (exri * vk)) * std::complex<double>(0.0, 0.5);
+ sum2 = cc0;
+ sum3 = cc0;
+ for (int i14 = 1; i14 <= c4->nlem; i14++) { // GPU portable?
+ int ie = i14 + c4->nlem;
+ sum2 += (c1ao->am0m[i14 - 1][i14 - 1] + c1ao->am0m[ie - 1][ie - 1]);
+ sum3 += (c1ao->am0m[i14 - 1][ie - 1] + c1ao->am0m[ie - 1][i14 - 1]);
+ } // i14 loop
+ double sumpi = 0.0;
+ sumpd = cc0;
+ int nlemt = c4->nlem + c4->nlem;
+ for (int i16 = 1; i16 <= nlemt; i16++) {
+ for (int j16 = 1; j16 <= c4->nlem; j16++) {
+ int je = j16 + c4->nlem;
+ double rvalue = (
+ dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][j16 - 1]
+ + dconjg(c1ao->am0m[i16 - 1][je - 1]) * c1ao->am0m[i16 - 1][je - 1]
+ ).real();
+ sumpi += rvalue;
+ sumpd += (
+ dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][je - 1]
+ + dconjg(c1ao->am0m[i16 - 1][je - 1]) * c1ao->am0m[i16 - 1][j16 - 1]
+ );
+ } // j16 loop
+ } // i16 loop
+ if (inpol == 0) {
+ sum1 = sum2;
+ sum4 = sum3 * uim;
+ sum3 = -sum4;
+ sums1 = sumpi;
+ sums2 = sumpi;
+ sums3 = sumpd * uim;
+ sums4 = -sums3;
+ } else { // label 18
+ sum1 = sum2 + sum3;
+ sum2 = sum2 - sum3;
+ sum3 = cc0;
+ sum4 = cc0;
+ sums1 = sumpi - sumpd;
+ sums2 = sumpi + sumpd;
+ sums3 = cc0;
+ sums4 = cc0;
+ }
+ // label 20
+ c1ao->ecscm[0] = -cccs * sum2.real();
+ c1ao->ecscm[1] = -cccs * sum1.real();
+ c1ao->ecscpm[0] = -cccs * sum4;
+ c1ao->ecscpm[1] = -cccs * sum3;
+ c1ao->fsacm[0][0] = csam * sum2;
+ c1ao->fsacm[1][0] = csam * sum4;
+ c1ao->fsacm[1][1] = csam * sum1;
+ c1ao->fsacm[0][1] = csam * sum3;
+ c1ao->scscm[0] = cccs * sums1.real();
+ c1ao->scscm[1] = cccs * sums2.real();
+ c1ao->scscpm[0] = cccs * sums3;
+ c1ao->scscpm[1] = cccs * sums4;
+}
+
/*! \brief C++ porting of POLAR
*
* \param x: `double`
@@ -915,6 +1358,147 @@ void r3j000(int j2, int j3, C6 *c6) {
}
}
+/*! \brief C++ porting of RABA
+ *
+ * \param le: `int`
+ * \param am0m: Matrix of complex.
+ * \param w: Matrix of complex.
+ * \param tqce: `double **`
+ * \param tqcpe: Matrix of complex.
+ * \param tqcs: `double **`
+ * \param tqcps: Matrix of complex.
+ */
+void raba(
+ int le, std::complex<double> **am0m, std::complex<double> **w, double **tqce,
+ std::complex<double> **tqcpe, double **tqcs, std::complex<double> **tqcps
+) {
+ std::complex<double> **a, **ctqce, **ctqcs;
+ std::complex<double> acw, acwp, aca, acap, c1, c2, c3;
+ const std::complex<double> cc0(0.0, 0.0);
+ const std::complex<double> uim(0.0, 1.0);
+ const double sq2i = 1.0 / sqrt(2.0);
+ int nlem = le * (le + 2);
+ const int nlemt = nlem + nlem;
+ a = new std::complex<double>*[nlemt];
+ ctqce = new std::complex<double>*[2];
+ ctqcs = new std::complex<double>*[2];
+ for (int ai = 0; ai < nlemt; ai++) a[ai] = new std::complex<double>[2]();
+ for (int ci = 0; ci < 2; ci++) {
+ ctqce[ci] = new std::complex<double>[3]();
+ ctqcs[ci] = new std::complex<double>[3]();
+ }
+ for (int i20 = 1; i20 <= nlemt; i20++) {
+ int i = i20 - 1;
+ c1 = cc0;
+ c2 = cc0;
+ for (int j10 = 1; j10 <= nlemt; j10++) {
+ int j = j10 - 1;
+ c1 += (am0m[i][j] * w[j][0]);
+ c2 += (am0m[i][j] * w[j][1]);
+ } // j10 loop
+ a[i][0] = c1;
+ a[i][1] = c2;
+ } //i20 loop
+ int jpo = 2;
+ for (int ipo70 = 1; ipo70 <= 2; ipo70++) {
+ if (ipo70 == 2) jpo = 1;
+ int ipo = ipo70 - 1;
+ ctqce[ipo][0] = cc0;
+ ctqce[ipo][1] = cc0;
+ ctqce[ipo][2] = cc0;
+ tqcpe[ipo][0] = cc0;
+ tqcpe[ipo][1] = cc0;
+ tqcpe[ipo][2] = cc0;
+ ctqcs[ipo][0] = cc0;
+ ctqcs[ipo][1] = cc0;
+ ctqcs[ipo][2] = cc0;
+ tqcps[ipo][0] = cc0;
+ tqcps[ipo][1] = cc0;
+ tqcps[ipo][2] = cc0;
+ for (int l60 = 1; l60 <= le; l60 ++) {
+ int lpo = l60 + 1;
+ int il = l60 * lpo;
+ int ltpo = l60 + lpo;
+ for (int im60 = 1; im60 <= ltpo; im60++) {
+ int m = im60 - lpo;
+ int i = m + il;
+ int ie = i + nlem;
+ int mmmu = m + 1;
+ int mmmmu = (mmmu > 0) ? mmmu : -mmmu;
+ double rmu = 0.0;
+ if (mmmmu <= l60) {
+ int immu = mmmu + il;
+ int immue = immu + nlem;
+ rmu = -sqrt(1.0 * (l60 + mmmu) * (l60 - m)) * sq2i;
+ acw = dconjg(a[i - 1][ipo]) * w[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * w[immue - 1][ipo];
+ acwp = dconjg(a[i - 1][ipo]) * w[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * w[immue - 1][jpo - 1];
+ aca = dconjg(a[i - 1][ipo]) * a[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * a[immue - 1][ipo];
+ acap = dconjg(a[i - 1][ipo]) * a[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * a[immue - 1][jpo - 1];
+ ctqce[ipo][0] += (acw * rmu);
+ tqcpe[ipo][0] += (acwp * rmu);
+ ctqcs[ipo][0] += (aca * rmu);
+ tqcps[ipo][0] += (acap * rmu);
+ }
+ // label 30
+ rmu = -1.0 * m;
+ acw = dconjg(a[i - 1][ipo]) * w[i - 1][ipo] + dconjg(a[ie - 1][ipo]) * w[ie - 1][ipo];
+ acwp = dconjg(a[i - 1][ipo]) * w[i - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * w[ie - 1][jpo - 1];
+ aca = dconjg(a[i - 1][ipo]) * a[i - 1][ipo] + dconjg(a[ie - 1][ipo]) * a[ie - 1][ipo];
+ acap = dconjg(a[i - 1][ipo]) * a[i - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * a[ie - 1][jpo - 1];
+ ctqce[ipo][1] += (acw * rmu);
+ tqcpe[ipo][1] += (acwp * rmu);
+ ctqcs[ipo][1] += (aca * rmu);
+ tqcps[ipo][1] += (acap * rmu);
+ mmmu = m - 1;
+ mmmmu = (mmmu > 0) ? mmmu : -mmmu;
+ if (mmmmu <= l60) {
+ int immu = mmmu + il;
+ int immue = immu + nlem;
+ rmu = sqrt(1.0 * (l60 - mmmu) * (l60 + m)) * sq2i;
+ acw = dconjg(a[i - 1][ipo]) * w[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * w[immue - 1][ipo];
+ acwp = dconjg(a[i - 1][ipo]) * w[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * w[immue - 1][jpo - 1];
+ aca = dconjg(a[i - 1][ipo]) * a[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * a[immue - 1][ipo];
+ acap = dconjg(a[i - 1][ipo]) * a[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * a[immue - 1][jpo - 1];
+ ctqce[ipo][2] += (acw * rmu);
+ tqcpe[ipo][2] += (acwp * rmu);
+ ctqcs[ipo][2] += (aca * rmu);
+ tqcps[ipo][2] += (acap * rmu);
+ } // ends im60 loop
+ } // im60 loop
+ } // l60 loop
+ } // ipo70 loop
+ for (int ipo78 = 1; ipo78 <= 2; ipo78++) {
+ int ipo = ipo78 - 1;
+ tqce[ipo][0] = (ctqce[ipo][0] - ctqce[ipo][2]).real() * sq2i;
+ tqce[ipo][1] = ((ctqce[ipo][0] + ctqce[ipo][2]) * uim).real() * sq2i;
+ tqce[ipo][2] = ctqce[ipo][1].real();
+ c1 = tqcpe[ipo][0];
+ c2 = tqcpe[ipo][1];
+ c3 = tqcpe[ipo][2];
+ tqcpe[ipo][0] = (c1 - c3) * sq2i;
+ tqcpe[ipo][1] = (c1 + c3) * (uim * sq2i);
+ tqcpe[ipo][2] = c2;
+ tqcs[ipo][0] = -sq2i * (ctqcs[ipo][0] - ctqcs[ipo][2]).real();
+ tqcs[ipo][1] = -sq2i * ((ctqcs[ipo][0] + ctqcs[ipo][2]) * uim).real();
+ tqcs[ipo][2] = -1.0 * ctqcs[ipo][1].real();
+ c1 = tqcps[ipo][0];
+ c2 = tqcps[ipo][1];
+ c3 = tqcps[ipo][2];
+ tqcps[ipo][0] = -(c1 - c3) * sq2i;
+ tqcps[ipo][1] = -(c1 + c3) * (uim * sq2i);
+ tqcps[ipo][2] = -c2;
+ } // ipo78 loop
+ // Clean memory
+ for (int ai = 0; ai < nlemt; ai++) delete[] a[ai];
+ for (int ci = 0; ci < 2; ci++) {
+ delete[] ctqce[ci];
+ delete[] ctqcs[ci];
+ }
+ delete[] a;
+ delete[] ctqce;
+ delete[] ctqcs;
+}
+
/*! \brief C++ porting of SCR0
*
* \param vk: `double` QUESTION: definition?
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index af2b222ea0a3948bbfc412247c6f30a9966a353a..621918756faec5e05838cd11b191e63abd85a179 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -123,6 +123,8 @@ C1_AddOns::C1_AddOns(C4 *c4) {
ecscpm = new complex<double>[2]();
allocate_vectors(c4);
sscs = new double[nsph]();
+ ecscm = new double[2]();
+ scscm = new double[2]();
}
C1_AddOns::~C1_AddOns() {
@@ -153,6 +155,8 @@ C1_AddOns::~C1_AddOns() {
delete[] ecscp;
delete[] scscpm;
delete[] ecscpm;
+ delete[] ecscm;
+ delete[] scscm;
}
void C1_AddOns::allocate_vectors(C4 *c4) {
diff --git a/src/np_cluster.cpp b/src/np_cluster.cpp
index 16abf18d024e5c208e7f251cf1cc183e2548545b..e151aed745f5e7633bd23531593cc11f13441dcf 100644
--- a/src/np_cluster.cpp
+++ b/src/np_cluster.cpp
@@ -10,7 +10,6 @@
using namespace std;
extern void cluster();
-extern void sphere();
/*! \brief Main program entry point.
*
@@ -21,8 +20,6 @@ extern void sphere();
* the configuration and runs the main program.
*/
int main(int argc, char **argv) {
- bool is_sphere = false;
- if (is_sphere) sphere();
- else cluster();
+ cluster();
return 0;
}
diff --git a/src/np_sphere.cpp b/src/np_sphere.cpp
index b6ba389eb191b88c11c44887522cc0f114226913..2e7d05f251b44f4d66550efce2523893e197971e 100644
--- a/src/np_sphere.cpp
+++ b/src/np_sphere.cpp
@@ -9,7 +9,6 @@
using namespace std;
-extern void cluster();
extern void sphere();
/*! \brief Main program entry point.
@@ -21,8 +20,6 @@ extern void sphere();
* the configuration and runs the main program.
*/
int main(int argc, char **argv) {
- bool is_sphere = true;
- if (is_sphere) sphere();
- else cluster();
+ sphere();
return 0;
}
diff --git a/src/sphere/Makefile b/src/sphere/Makefile
index 76e16f0b52ee8cd5e4c91254ef2bd20cd0e0e48e..c42d9b04f701cfb823a0c443581c96c23207d920 100644
--- a/src/sphere/Makefile
+++ b/src/sphere/Makefile
@@ -14,8 +14,8 @@ edfb: edfb.o
sph: sph.o
$(FC) $(FCFLAGS) -o $(BUILDDIR)/sph $(BUILDDIR)/sph.o $(LFLAGS)
-np_sphere: $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o $(BUILDDIR)/cluster.o
- $(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_sphere $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o $(BUILDDIR)/cluster.o
+np_sphere: $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o
+ $(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_sphere $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o
$(BUILDDIR)/np_sphere.o:
$(CXX) $(CXXFLAGS) -c ../np_sphere.cpp -o $(BUILDDIR)/np_sphere.o
@@ -29,9 +29,6 @@ $(BUILDDIR)/Configuration.o:
$(BUILDDIR)/Parsers.o:
$(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o
-$(BUILDDIR)/cluster.o:
- $(CXX) $(CXXFLAGS) -c ../cluster/cluster.cpp -o $(BUILDDIR)/cluster.o
-
$(BUILDDIR)/sphere.o:
$(CXX) $(CXXFLAGS) -c sphere.cpp -o $(BUILDDIR)/sphere.o