diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e57added3ad108c00998cf37fa57770cba1a9d91
--- /dev/null
+++ b/src/cluster/cluster.cpp
@@ -0,0 +1,279 @@
+#include <cstdio>
+#include <fstream>
+#include <string>
+#include <complex>
+#ifndef INCLUDE_CONFIGURATION_H_
+#include "../include/Configuration.h"
+#endif
+#ifndef INCLUDE_CLU_SUBS_H_
+#include "../include/clu_subs.h"
+#endif
+
+using namespace std;
+
+/*
+ * >>> WARNING: works only for IDFC >= 0, as the original code <<<
+ *
+ */
+
+//! \brief C++ implementation of CLU
+void cluster() {
+ printf("INFO: making legacy configuration ...\n");
+ ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("DEDFB_clu");
+ 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");
+ if (sconf->number_of_spheres == gconf->number_of_spheres) {
+ // Shortcuts to variables stored in configuration objects
+ int nsph = gconf->number_of_spheres;
+ int mxndm = gconf->mxndm;
+ int inpol = gconf->in_pol;
+ int npnt = gconf->npnt;
+ int npntts = gconf->npntts;
+ int iavm = gconf->iavm;
+ int isam = gconf->meridional_type;
+ int nxi = sconf->number_of_scales;
+ double th = gconf->in_theta_start;
+ double thstp = gconf->in_theta_step;
+ double thlst = gconf->in_theta_end;
+ double ths = gconf->sc_theta_start;
+ double thsstp = gconf->sc_theta_step;
+ double thslst = gconf->sc_theta_end;
+ double ph = gconf->in_phi_start;
+ double phstp = gconf->in_phi_step;
+ double phlst = gconf->in_phi_end;
+ double phs = gconf->sc_phi_start;
+ double phsstp = gconf->sc_phi_step;
+ double phslst = gconf->sc_phi_end;
+ double wp = sconf->wp;
+ // Global variables for CLU
+ double pi = 2.0 * acos(0.0);
+ int lm = gconf->l_max;
+ if (gconf->li > lm) lm = gconf->li;
+ if (gconf->le > lm) lm = gconf->le;
+ C1 *c1 = new C1(nsph, lm, sconf->nshl_vec, sconf->iog_vec);
+ C3 *c3 = new C3();
+ for (int c1i = 0; c1i < nsph; c1i++) {
+ c1->rxx[c1i] = gconf->sph_x[c1i];
+ c1->ryy[c1i] = gconf->sph_y[c1i];
+ c1->rzz[c1i] = gconf->sph_z[c1i];
+ c1->ros[c1i] = sconf->radii_of_spheres[c1i];
+ int iogi = c1->iog[c1i];
+ if (iogi >= c1i + 1) {
+ double gcss = pi * c1->ros[c1i] * c1->ros[c1i];
+ c1->gcsv[c1i] = gcss;
+ int nsh = c1->nshl[c1i];
+ for (int j16 = 1; j16 <= nsh; j16++) {
+ c1->rc[c1i][j16- 1] = sconf->rcf[c1i][j16] * c1->ros[c1i];
+ }
+ c3->gcs += c1->gcsv[c1i - 1];
+ }
+ }
+ C4 *c4 = new C4;
+ c4->li = gconf->li;
+ c4->le = gconf->le;
+ c4->lm = lm;
+ c4->lmpo = c4->lm + 1;
+ c4->litpo = 2 * gconf->li + 1;
+ c4->litpos = c4->litpo * c4->litpo;
+ c4->lmtpo = gconf->li + gconf->le + 1;
+ c4->lmtpos = c4->lmtpo * c4->lmtpo;
+ c4->nlim = c4->li * (c4->li + 2);
+ c4->nlem = c4->le * (c4->le + 2);
+ c4->nsph = nsph;
+ C6 *c6 = new C6(c4->lmtpo);
+ C1_AddOns *c1ao = new C1_AddOns(c4);
+ FILE *output = fopen("c_OCLU", "w");
+ double ****zpv = new double***[c4->lm];
+ for (int zi = 0; zi < c4->lm; zi++) {
+ zpv[zi] = new double**[3];
+ for (int zj = 0; zj < 3; zj++) {
+ zpv[zi][zj] = new double*[2];
+ zpv[zi][zj][0] = new double[2];
+ zpv[zi][zj][1] = new double[2];
+ }
+ }
+ int jer = 0, lcalc = 0;
+ complex<double> arg = complex<double>(0.0, 0.0);
+ int max_ici = 0;
+ for (int insh = 0; insh < nsph; insh++) {
+ int nsh = sconf->nshl_vec[insh];
+ int ici = (nsh + 1) / 2;
+ if (ici > max_ici) max_ici = ici;
+ }
+ C2 *c2 = new C2(nsph, max_ici, npnt, npntts);
+ complex<double> **am = new complex<double>*[mxndm];
+ for (int ai = 0; ai < mxndm; ai++) am[ai] = new complex<double>[mxndm];
+ // 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
+ );
+ fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
+ for (int ri = 0; ri < nsph; ri++) fprintf(output, "%17.8lE%17.8lE%17.8lE\n",
+ gconf->sph_x[ri], gconf->sph_y[ri], gconf->sph_z[ri]
+ );
+ fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
+ fprintf(
+ output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
+ th, thstp, thlst, ths, thsstp, thslst
+ );
+ fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
+ fprintf(
+ output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
+ ph, phstp, phlst, phs, phsstp, phslst
+ );
+ fprintf(output, " READ(IR,*)JWTM\n");
+ fprintf(output, " %5d\n", gconf->jwtm);
+ fprintf(output, " READ(ITIN)NSPHT\n");
+ fprintf(output, " READ(ITIN)(IOG(I),I=1,NSPH)\n");
+ fprintf(output, " READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
+ fprintf(output, " READ(ITIN)(XIV(I),I=1,NXI)\n");
+ fprintf(output, " READ(ITIN)NSHL(I),ROS(I)\n");
+ fprintf(output, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");
+ fprintf(output, " \n");
+ double small = 1.0e-3;
+ int nth = 0, nph = 0;
+ if (thstp != 0.0) nth = int((thlst - th) / thstp + small);
+ nth++;
+ if (phstp != 0.0) nph = int((phlst - ph) / phstp + small);
+ nph++;
+ int nths = 0, nphs = 0;
+ double thsca = 0.0;
+ if (isam > 1) {
+ nths = 1;
+ thsca = ths - th;
+ } else { // ISAM <= 1
+ if (thsstp == 0.0) nths = 0;
+ else nths = int ((thslst - ths) / thsstp + small);
+ nths++;
+ }
+ if (isam >= 1) {
+ nphs = 1;
+ } else {
+ if (phsstp == 0.0) nphs = 0;
+ else nphs = int((phslst - phs) / phsstp + small);
+ nphs++;
+ }
+ int nk = nth * nph;
+ int nks = nths * nphs;
+ int nkks = nk * nks;
+ double th1 = th, ph1 = ph, ths1 = ths, phs1 = phs;
+ str(c1, c1ao, c3, c4, c6);
+ thdps(c4->lm, zpv);
+ double exri = sqrt(sconf->exdc);
+ double vk = 0.0; // NOTE: Not really sure it should be initialized at 0
+ fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
+ fstream tppoan;
+ tppoan.open("c_TPPOAN", ios::out | ios::binary);
+ if (tppoan.is_open()) {
+ tppoan.write(reinterpret_cast<char *>(&iavm), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&inpol), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&nxi), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&nth), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&nph), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
+ int nlemt = c4->nlem + c4->nlem;
+ double wn = wp / 3.0e8;
+ double sqsfi = 1.0;
+ if (sconf->idfc < 0) {
+ vk = sconf->xip * wn;
+ fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
+ fprintf(output, " \n");
+ }
+ for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {
+ int jaw = 1;
+ fprintf(output, "========== JXI =%3d ====================\n", jxi488);
+ double xi = sconf->scale_vec[jxi488 - 1];
+ double vkarg = 0.0;
+ if (sconf->idfc >= 0) {
+ vk = xi * wn;
+ vkarg = vk;
+ fprintf(output, " VK=%15.7lE, XI=%15.7lE\n", vk, xi);
+ } else {
+ vkarg = xi * vk;
+ sqsfi = 1.0 / (xi * xi);
+ fprintf(output, " XI=%15.7lE\n", xi);
+ }
+ // Would call HJV(EXRI, VKARG, JER, LCALC, ARG)
+ hjv(exri, vkarg, jer, lcalc, arg, c1, c1ao, c4);
+ //printf("INFO: calculation went up to %d and jer = %d\n", lcalc, jer);
+ //printf("INFO: arg = (%lE,%lE)\n", arg.real(), arg.imag());
+ if (jer != 0) {
+ fprintf(output, " STOP IN HJV\n");
+ break; // jxi488 loop: goes to memory cleaning and return
+ }
+ for (int i132 = 1; i132 <= nsph; i132++) {
+ int iogi = c1->iog[i132 - 1];
+ if (iogi != i132) {
+ for (int l123 = 1; l123 <= gconf->li; l123++) {
+ c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];
+ c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];
+ } // l123 loop
+ } else {
+ int nsh = sconf->nshl_vec[i132 - 1];
+ int ici = (nsh + 1) / 2;
+ int size_dc0 = (nsh % 2 == 0) ? ici + 1 : ici;
+ if (sconf->idfc == 0) {
+ for (int ic = 0; ic < ici; ic++)
+ c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][jxi488 - 1];
+ } else {
+ if (jxi488 == 1) {
+ for (int ic = 0; ic < ici; ic++)
+ c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][0];
+ }
+ }
+ if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
+ dme(
+ c4->li, i132, npnt, npntts, vkarg, sconf->exdc, exri,
+ c1, c2, jer, lcalc, arg
+ );
+ //printf("INFO: DME returned jer = %d , lcalc = %d and arg = (%lE, %lE)\n",
+ // jer, lcalc, arg.real(), arg.imag());
+ if (jer != 0) {
+ fprintf(output, " STOP IN DME\n");
+ break;
+ }
+ }
+ if (jer != 0) break;
+ } // i132 loop
+ // Would call CMS(AM)
+
+ if (jer != 0) break;
+ } // jxi488 loop
+ tppoan.close();
+ } else { // In case TPPOAN could not be opened. Should never happen.
+ printf("ERROR: failed to open TPPOAN file.\n");
+ }
+ fclose(output);
+ // Clean memory
+ delete c1;
+ delete c1ao;
+ delete c3;
+ delete c4;
+ delete c6;
+ for (int zi = c4->lm - 1; zi > -1; zi--) {
+ for (int zj = 2; zj > -1; zj--) {
+ delete[] zpv[zi][zj][1];
+ delete[] zpv[zi][zj][0];
+ delete[] zpv[zi][zj];
+ }
+ delete[] zpv[zi];
+ }
+ delete[] zpv;
+ for (int ai = mxndm - 1; ai > -1; ai--) delete[] am[ai];
+ delete[] am;
+ } else { // NSPH mismatch between geometry and scatterer configurations.
+ throw UnrecognizedConfigurationException(
+ "Inconsistent geometry and scatterer configurations."
+ );
+ }
+ delete sconf;
+ delete gconf;
+ printf("Done.\n");
+}
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
new file mode 100644
index 0000000000000000000000000000000000000000..52568582d0184916ddc56603b05b538fff51f3b6
--- /dev/null
+++ b/src/include/clu_subs.h
@@ -0,0 +1,740 @@
+/*! \file clu_subs.h
+ *
+ * \brief C++ porting of CLU functions and subroutines.
+ *
+ * Remember that FORTRAN passes arguments by reference, so, every time we use
+ * a subroutine call, we need to add a referencing layer to the C++ variable.
+ * All the functions defined below need to be properly documented and ported
+ * to C++.
+ *
+ * Currently, only basic documenting information about functions and parameter
+ * types are given, to avoid doxygen warning messages.
+ */
+
+#ifndef INCLUDE_COMMONS_H_
+#include "Commons.h"
+#endif
+
+#ifndef INCLUDE_CLU_SUBS_H_
+#define INCLUDE_CLU_SUBS_H_
+
+#include <complex>
+
+// >>> DECLARATION OF SPH_SUBS <<<
+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 rbf(int n, double x, int &nm, double sj[]);
+extern void rnf(int n, double x, int &nm, double sy[]);
+extern void thdps(int lm, double ****zpv);
+// >>> END OF SPH_SUBS DECLARATION <<<
+
+/*! \brief C++ porting of CGEV
+ *
+ * \param ipamo: `int`
+ * \param mu: `int`
+ * \param l: `int`
+ * \param m: `int`
+ */
+double cgev(int ipamo, int mu, int l, int m) {
+ double result = 0.0;
+ double xd = 0.0, xn = 0.0;
+ if (ipamo == 0) {
+ if (m != 0 || mu != 0) { // label 10
+ if (mu != 0) {
+ xd = 2.0 * l * (l + 1);
+ if (mu <= 0) {
+ xn = 1.0 * (l + m) * (l - m + 1);
+ result = sqrt(xn / xd);
+ } else { // label 15
+ xn = 1.0 * (l - m) * (l + m + 1);
+ result = -sqrt(xn / xd);
+ }
+ } else { // label 20
+ xd = 1.0 * (l + 1) * l;
+ xn = -1.0 * m;
+ result = xn / sqrt(xd);
+ }
+ }
+ } else { // label 30
+ xd = 2.0 * l * (l * 2 - 1);
+ if (mu < 0) { // label 35. CHECK: is clu.f code line 2466 a switch?
+ xn = 1.0 * (l - 1 + m) * (l + m);
+ } else if (mu == 0) { // label 40
+ xn = 2.0 * (l - m) * (l + m);
+ } else { // mu > 0, label 45
+ xn = 1.0 * (l - 1 - m) * (l - m);
+ }
+ result = sqrt(xn / xd);
+ }
+ return result;
+}
+
+/*! \brief C++ porting of R3JJR
+ *
+ * \param j2: `int`
+ * \param j3: `int`
+ * \param m2: `int`
+ * \param m3: `int`
+ * \param c6: `C6 *`
+ */
+void r3jjr(int j2, int j3, int m2, int m3, C6 *c6) {
+
+}
+
+
+/*! \brief C++ porting of GHIT
+ *
+ * \param ihi: `int`
+ * \param ipamo: `int`
+ * \param nbl: `int`
+ * \param l1: `int`
+ * \param m1: `int`
+ * \param l2: `int`
+ * \param m2: `int`
+ * \param c1: `C1 *`
+ * \param c1ao: `C1_AddOns *`
+ * \param c4: `C4 *`
+ * \param c6: `C6 *`
+ */
+std::complex<double> ghit(
+ int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
+ C1_AddOns *c1ao, C4 *c4, C6 *c6
+) {
+ /* NBL identifies transfer vector going from N2 to N1;
+ * IHI=0 for Hankel, IHI=1 for Bessel, IHI=2 for Bessel from origin;
+ * depending on IHI, IPAM=0 gives H or I, IPAM= 1 gives K or L. */
+ const std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
+ const std::complex<double> uim = std::complex<double>(0.0, 1.0);
+ std::complex<double> result = cc0;
+ if (ihi == 2) {
+ if (!(c1->rxx[nbl - 1] != 0.0 || c1->ryy[nbl - 1] != 0.0 || c1->rzz[nbl - 1] != 0.0)) {
+ if (ipamo != 0) return result;
+ if (l1 == l2 && m1 == m2) result = std::complex<double>(1.0, 0.0);
+ return result;
+ }
+ }
+ int l1mp = l1 - ipamo;
+ int l1po = l1 + 1;
+ int m1mm2 = m1 - m2;
+ int m1mm2m = (m1mm2 > 0) ? m1mm2 + 1 : -m1mm2 + 1;
+ int lminpo = (l2 - l1mp > 0) ? l2 - l1mp + 1 : l1mp - l2 + 1;
+ int lmaxpo = l2 + l1mp + 1;
+ int i3j0in = c1ao->ind3j[l1mp][l2 - 1];
+ int ilin = -1;
+ if (m1mm2m > lminpo && (m1mm2m - lminpo) % 2 != 0) ilin = 0;
+ int isn = 1;
+ if (m1 % 2 != 0) isn *= -1;
+ if (lminpo % 2 == 0) {
+ isn *= -1;
+ if (l2 > l1mp) isn *= -1;
+ }
+ int nblmo = nbl - 1;
+ if (ihi != 2) {
+ int nbhj = nblmo * c4->litpo;
+ int nby = nblmo * c4->litpos;
+ if (ihi != 1) {
+ for (int jm24 = 1; jm24 <= 3; jm24++) {
+ std::complex<double> csum = cc0;
+ int mu = jm24 - 2;
+ int mupm1 = mu + m1;
+ int mupm2 = mu + m2;
+ if (mupm1 < -l1mp || mupm1 > l1mp || mupm2 < -l2 || mupm2 > l2)
+ continue; //jm24 loop
+ int jsn = -isn;
+ if (mu == 0) jsn = isn;
+ double cr = cgev(ipamo, mu, l1, m1) * cgev(0, mu, l2, m2);
+ int i3j0 = i3j0in;
+ if (mupm1 == 0 && mupm2 == 0) {
+ int lt14 = lminpo;
+ while (lt14 <= lmaxpo) {
+ i3j0 += 1;
+ int l3 = lt14 - 1;
+ int ny = l3 * l3 + lt14;
+ double aors = 1.0 * (l3 + lt14);
+ double f3j = (
+ c1ao->v3j0[i3j0 - 1] * c1ao->v3j0[i3j0 - 1] *
+ sqrt(aors)
+ ) * jsn;
+ std::complex<double> cfun = (
+ c1ao->vh[nbhj + lt14 - 1] * c1ao->vyhj[nby + ny - 1]
+ ) * f3j;
+ csum += cfun;
+ jsn *= - 1;
+ lt14 += 2;
+ }
+ } else { // label 16
+ // Would call R3JJR(L1MP, L2, -MUPM1, MUPM2)
+ int il = ilin;
+ int lt20 = lminpo;
+ while (lt20 <= lmaxpo) {
+ i3j0 += 1;
+ if (m1mm2m <= lt20) {
+ il = il + 2;
+ int l3 = lt20 - 1;
+ int ny = l3 * l3 + lt20 + m1mm2;
+ double aors = 1.0 * (l3 + lt20);
+ double f3j = (
+ c6->rac3j[il - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)
+ ) * jsn;
+ std::complex<double> cfun = (
+ c1ao->vh[nbhj + lt20 - 1] * c1ao->vyhj[nby + ny - 1]
+ ) * f3j;
+ csum += cfun;
+ }
+ jsn *= -1;
+ lt20 += 2;
+ }
+ }
+ csum *= cr;
+ result += csum;
+ } // jm24 loop
+ } else { // label 30, IHI = 1
+ for (int jm44 = 1; jm44 <= 3; jm44 ++) {
+ std::complex<double> csum = cc0;
+ int mu = jm44 - 2;
+ int mupm1 = mu + m1;
+ int mupm2 = mu + m2;
+ if (mupm1 >= -l1mp && mupm1 <= l1mp && mupm2 >= -l2 && mupm2 >= l2) {
+ int jsn = -isn;
+ if (mu == 0) jsn = isn;
+ double cr = cgev(ipamo, mu, l1, m1) * cgev(0, mu, l2, m2);
+ int i3j0 = i3j0in;
+ if (mupm1 == 0 && mupm2 == 0) {
+ int lt34 = lminpo;
+ while (lt34 <= lmaxpo) {
+ i3j0++;
+ int l3 = lt34 - 1;
+ int ny = l3 * l3 + lt34;
+ double aors = 1.0 * (l3 + lt34);
+ std::complex<double> f3j = (
+ c1ao->v3j0[i3j0 - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)
+ ) * jsn;
+ std::complex<double> cfun = (
+ c1ao->vj[nbhj + lt34 - 1] * c1ao->vyhj[nby + ny - 1]
+ ) * f3j;
+ csum += cfun;
+ jsn *= -1;
+ lt34 += 2;
+ }
+ } else { // label 36
+ // Would call R3JJR
+ int il = ilin;
+ int lt40 = lminpo;
+ while (lt40 <= lmaxpo) {
+ i3j0++;
+ if (m1mm2m <= lt40) {
+ il += 2;
+ int l3 = lt40 - 1;
+ int ny = l3 * l3 + lt40 + m1mm2;
+ double aors = 1.0 * (l3 + lt40);
+ std::complex<double> f3j = (
+ c6->rac3j[il - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)
+ ) * jsn;
+ std::complex<double> cfun = (
+ c1ao->vj[nbhj + lt40 - 1] * c1ao->vyhj[nby + ny - 1]
+ ) * f3j;
+ csum += cfun;
+ }
+ jsn *= -1;
+ lt40 += 2;
+ }
+ }
+ // label 42
+ csum *= cr;
+ }
+ result += csum;
+ } // jm44 loop
+ }
+ } else { // label 50, IHI = 2
+ int nbhj = nblmo * c4->lmtpo;
+ int nby = nblmo * c4->lmtpos;
+ for (int jm64 =1; jm64 <= 3; jm64++) {
+ std::complex<double> csum = cc0;
+ int mu = jm64 - 2;
+ int mupm1 = mu + m1;
+ int mupm2 = mu + m2;
+ if (mupm1 >= -l1mp && mupm1 < l1mp && mupm2 >= -l2 && mupm2 <= l2) {
+ int jsn = - isn;
+ if (mu == 0) jsn = isn;
+ double cr = cgev(ipamo, mu, l1, m1) * cgev(0, mu, l2, m2);
+ int i3j0 = i3j0in;
+ if (mupm1 == 0 && mupm2 == 0) {
+ int lt54 = lminpo;
+ while (lt54 <= lmaxpo) {
+ i3j0++;
+ int l3 = lt54 - 1;
+ int ny = l3 * l3 + lt54;
+ double aors = 1.0 * (l3 + lt54);
+ std::complex<double> f3j = (
+ c1ao->v3j0[i3j0 - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)
+ ) * jsn;
+ std::complex<double> cfun = (
+ c1ao->vj0[nbhj + lt54 - 1] * c1ao->vyj0[nby + ny - 1]
+ ) * f3j;
+ csum += cfun;
+ jsn *= -1;
+ lt54 += 2;
+ }
+ } else { // label 56
+ // Would call R3JJR(L1MP, L2, -MUPM1, MUPM2)
+ int il = ilin;
+ int lt60 = lminpo;
+ while (lt60 <= lmaxpo) {
+ i3j0++;
+ if (m1mm2m <= lt60) {
+ il += 2;
+ int l3 = lt60 - 1;
+ int ny = l3 * l3 + lt60 + m1mm2;
+ double aors = 1.0 * (l3 + lt60);
+ std::complex<double> f3j = (
+ c6->rac3j[il - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)
+ ) * jsn;
+ std::complex<double> cfun = (
+ c1ao->vj0[nbhj + lt60 - 1] * c1ao->vyj0[nby + ny - 1]
+ ) * f3j;
+ csum += cfun;
+ }
+ jsn *= -1;
+ lt60 += 2;
+ }
+ }
+ // label 62
+ csum *= cr;
+ }
+ result += csum;
+ } // jm64 loop
+ }
+ // label 70
+ const double four_pi = acos(0.0) * 8.0;
+ if (ipamo != 1) {
+ double cr = sqrt(four_pi * (l1 + l1po) * (l2 + l2 + 1));
+ result *= cr;
+ } else {
+ double cr = sqrt(four_pi * (l1 + l1mp) * (l1 + l1po) * (l2 + l2 + 1) / l1po);
+ result *= (cr * uim);
+ }
+ return result;
+}
+
+/*! \brief C++ porting of CMS
+ *
+ * \param am: Matrix of complex.
+ * \param c1: `C1 *`
+ * \param c1ao: `C1_AddOns *`
+ * \param c4: `C4 *`
+ * \param c6: `C6 *`
+ */
+void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
+ std::complex<double> dm, de, cgh, cgk;
+ const std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
+ int ndi = c4->nsph * c4->nlim;
+ int nbl = 0;
+ int nsphmo = c4->nsph - 1;
+ for (int n1 = 1; n1 <= nsphmo; n1++) { // GPU portable?
+ int in1 = (n1 - 1) * c4->nlim;
+ int n1po = n1 + 1;
+ for (int n2 = n1po; n2 <= c4->nsph; n2++) {
+ int in2 = (n2 - 1) * c4->nlim;
+ nbl++;
+ for (int l1 = 1; l1 <= c4->li; l1++) {
+ int l1po = l1 + 1;
+ int il1 = l1po * l1;
+ int l1tpo = l1po + l1;
+ for (int im1 =1; im1 <= l1tpo; il1++) {
+ int m1 = im1 - l1po;
+ int ilm1 = il1 + m1;
+ int ilm1e = ilm1 + ndi;
+ int i1 = in1 + ilm1;
+ int i1e = in1 + ilm1e;
+ int j1 = in2 + ilm1;
+ int j1e = in2 + ilm1e;
+ for (int l2 =1; l2 <= c4->li; l2++) {
+ int l2po = l2 + 1;
+ int il2 = l2po * l2;
+ int l2tpo = l2po + l2;
+ int ish = ((l2 + l1) % 2 == 0) ? 1 : -1;
+ int isk = -ish;
+ for (int im2 = 1; im2 <= l2tpo; im2++) {
+ int m2 = im2 - l2po;
+ int ilm2 = il2 + m2;
+ int ilm2e = ilm2 + ndi;
+ int i2 = in2 + ilm2;
+ int i2e = in2 + ilm2e;
+ int j2 = in1 + ilm2;
+ int j2e = in1 + ilm2e;
+ cgh = ghit(0, 0, nbl, l1, m1, l2, m2, c1, c1ao, c4, c6);
+ cgk = ghit(0, 1, nbl, l1, m1, l2, m2, c1, c1ao, c4, c6);
+ am[i1 - 1][i2 - 1] = cgh;
+ am[i1 - 1][i2e - 1] = cgk;
+ am[i1e - 1][i2 - 1] = cgk;
+ am[i1e - 1][i2e - 1] = cgh;
+ am[j1 - 1][j2 - 1] = cgh * (1.0 * ish);
+ am[j1 - 1][j2e - 1] = cgk * (1.0 * isk);
+ am[j1e - 1][j2 - 1] = cgk * (1.0 * isk);
+ am[j1e - 1][j2e - 1] = cgh * (1.0 * ish);
+ }
+ }
+ } // im1 loop
+ } // l1 loop
+ } // n2 loop
+ } // n1 loop
+ for (int n1 = 1; n1 <= c4->nsph; n1++) { // GPU portable?
+ int in1 = (n1 - 1) * c4->nlim;
+ for (int l1 = 1; l1 <= c4->li; l1++) {
+ dm = c1->rmi[l1 - 1][n1 - 1];
+ de = c1->rei[l1 - 1][n1 - 1];
+ int l1po = l1 + 1;
+ int il1 = l1po * l1;
+ int l1tpo = l1po + l1;
+ for (int im1 = 1; im1 <= l1tpo; im1++) {
+ int m1 = im1 - l1po;
+ int ilm1 = il1 + m1;
+ int i1 = in1 + ilm1;
+ int i1e = i1 + ndi;
+ for (int ilm2 = 1; ilm2 <= c4->nlim; ilm2++) {
+ int i2 = in1 + ilm2;
+ int i2e = i2 + ndi;
+ am[i1 - 1][i2 - 1] = cc0;
+ am[i1 - 1][i2e - 1] = cc0;
+ am[i1e - 1][i2 - 1] = cc0;
+ am[i1e - 1][i2e - 1] = cc0;
+ }
+ am[i1 - 1][i1 - 1] = dm;
+ am[i1e - 1][i1e - 1] = de;
+ } // im1 loop
+ } // l1 loop
+ } // n1 loop
+}
+
+/*! \brief C++ porting of HJV
+ *
+ * \param exri: `double`
+ * \param vk: `double`
+ * \param jer: `int &`
+ * \param lcalc: `int &`
+ * \param arg: `complex\<double\> &`
+ * \param c1: `C1 *`
+ * \param c1ao: `C1_AddOns *`
+ * \param c4: `C4 *`
+ */
+void hjv(
+ double exri, double vk, int &jer, int &lcalc, std::complex<double> &arg,
+ C1 *c1, C1_AddOns *c1ao, C4 *c4
+) {
+ int nsphmo = c4->nsph - 1;
+ int lit = c4->li + c4->li;
+ int lmt = c4->lm + c4->lm;
+ const int rfj_size = (lit > lmt) ? lit : lmt;
+ const int rfn_size = c4->litpo;
+ double *rfj, *rfn;
+ rfj = new double[rfj_size];
+ rfn = new double[rfn_size];
+ jer = 0;
+ int ivhb = 0;
+ for (int nf40 = 1; nf40 <= nsphmo; nf40++) { // GPU portable?
+ int nfpo = nf40 + 1;
+ for (int ns40 = nfpo; ns40 <= c4->nsph; ns40++) {
+ double rx = c1->rxx[nf40 - 1] - c1->rxx[ns40 - 1];
+ double ry = c1->ryy[nf40 - 1] - c1->ryy[ns40 - 1];
+ double rz = c1->rzz[nf40 - 1] - c1->rzz[ns40 - 1];
+ double rr = sqrt(rx * rx + ry * ry + rz * rz);
+ double rarg = rr * vk * exri;
+ arg = std::complex<double>(rarg, 0.0);
+ rbf(lit, rarg, lcalc, rfj);
+ if (lcalc < lit) {
+ jer = 1;
+ delete[] rfj;
+ delete[] rfn;
+ return;
+ }
+ rnf(lit, rarg, lcalc, rfn);
+ if (lcalc < lit) {
+ jer = 2;
+ delete[] rfj;
+ delete[] rfn;
+ return;
+ }
+ for (int lpo38 = 1; lpo38 <= c4->litpo; lpo38++) {
+ c1ao->vh[lpo38 + ivhb - 1] = std::complex<double>(rfj[lpo38 - 1], rfn[lpo38 - 1]);
+ }
+ ivhb += c4->litpo;
+ } // ns40 loop
+ } // nf40 loop
+ ivhb = 0;
+ for (int nf50 = 1; nf50 <= c4->nsph; nf50++) {
+ double rx = c1->rxx[nf50 - 1];
+ double ry = c1->ryy[nf50 - 1];
+ double rz = c1->rzz[nf50 - 1];
+ if (!(rx == 0.0 && ry == 0.0 && rz == 0.0)) {
+ double rr = sqrt(rx * rx + ry * ry + rz * rz);
+ double rarg = rr * vk * exri;
+ rbf(lmt, rarg, lcalc, rfj);
+ if (lcalc < lmt) {
+ jer = 3;
+ delete[] rfj;
+ delete[] rfn;
+ return;
+ }
+ for (int lpo47 = 1; lpo47 <= c4->lmtpo; lpo47++) {
+ c1ao->vj0[lpo47 + ivhb - 1] = rfj[lpo47 - 1];
+ }
+ }
+ ivhb += c4->lmtpo;
+ } // nf50 loop
+ delete[] rfj;
+ delete[] rfn;
+}
+
+/*! \brief C++ porting of POLAR
+ *
+ * \param x: `double`
+ * \param y: `double`
+ * \param z: `double`
+ * \param r: `double &`
+ * \param cth: `double &`
+ * \param sth: `double &`
+ * \param cph: `double &`
+ * \param sph: `double &`
+ */
+void polar(
+ double x, double y, double z, double &r, double &cth, double &sth,
+ double &cph, double &sph
+) {
+ bool onx = (y == 0.0);
+ bool ony = (x == 0.0);
+ bool onz = (onx && ony);
+ double rho;
+ if (!(onz || onx || ony)) {
+ rho = sqrt(x * x + y * y);
+ cph = x / rho;
+ sph = y / rho;
+ } else {
+ if (onz) {
+ cph = 1.0;
+ sph = 0.0;
+ } else if (onx) {
+ rho = x;
+ cph = 1.0;
+ if (x < 0.0) {
+ rho *= -1.0;
+ cph *= -1.0;
+ }
+ sph = 0.0;
+ } else if (ony) {
+ rho = y;
+ sph = 1.0;
+ if (y < 0.0) {
+ rho *= -1.0;
+ sph *= -1.0;
+ }
+ cph = 0.0;
+ }
+ }
+ if (z == 0.0) {
+ if (!onz) {
+ r = rho;
+ cth = 0.0;
+ sth = 1.0;
+ return;
+ } else {
+ r = 0.0;
+ cth = 1.0;
+ sth = 0.0;
+ return;
+ }
+ } else {
+ if (!onz) {
+ r = sqrt(rho + z * z);
+ cth = z / r;
+ sth = rho / r;
+ } else {
+ r = z;
+ cth = 1.0;
+ sth = 0.0;
+ if (z < 0.0) {
+ r *= -1.0;
+ cth *= -1.0;
+ }
+ }
+ }
+}
+
+/*! \brief C++ porting of R3J000
+ *
+ * \param j2: `int`
+ * \param j3: `int`
+ * \param c6: `C6 *` Pointer to a C6 instance.
+ */
+void r3j000(int j2, int j3, C6 *c6) {
+ int jmx = j3 + j2;
+ if (jmx <= 0) {
+ c6->rac3j[0] = 1.0;
+ //printf("DEBUG: RAC3J(1) = %lE\n", c6->rac3j[0]);
+ return;
+ }
+ int jmn = j3 - j2;
+ if (jmn < 0) jmn *= -1;
+ int njmo = (jmx - jmn) / 2;
+ int jf = jmx + jmx + 1;
+ int isn = 1;
+ if (jmn % 2 != 0) isn = -1;
+ if (njmo <= 0) {
+ double sj = 1.0 * jf;
+ double cnr = (1 / sqrt(sj)) * isn;
+ c6->rac3j[0] = cnr;
+ //printf("DEBUG: RAC3J(1) = %lE\n", c6->rac3j[0]);
+ return;
+ }
+ double sjr = jf;
+ int jmxpos = (jmx + 1) * (jmx + 1);
+ int jmns = jmn * jmn;
+ int j1mo = jmx - 1;
+ int j1s = jmx * jmx; // a slight optimization was applied here
+ double cj = sqrt(1.0 * (jmxpos - j1s) * (j1s - jmns));
+ int j1mos = j1mo * j1mo;
+ double cjmo = sqrt(1.0 * (jmxpos - j1mos) * (j1mos - jmns));
+ if (njmo <= 1) {
+ c6->rac3j[0] = -cj / cjmo;
+ double sj = sjr + (c6->rac3j[0] * c6->rac3j[0]) * (jf - 4);
+ double cnr = (1.0 / sqrt(sj)) * isn;
+ c6->rac3j[1] = cnr;
+ c6->rac3j[0] *= cnr;
+ //printf("DEBUG: RAC3J(1) = %lE\n", c6->rac3j[0]);
+ //printf("DEBUG: RAC3J(2) = %lE\n", c6->rac3j[1]);
+ return;
+ }
+ int nj = njmo + 1;
+ int nmat = (nj + 1) / 2;
+ c6->rac3j[nj - 1] = 1.0;
+ c6->rac3j[njmo - 1] = -cj / cjmo;
+ if (nmat != njmo) {
+ int nbr = njmo - nmat;
+ for (int ibr45 = 1; ibr45 <= nbr; ibr45++) {
+ int irr = nj - ibr45;
+ jf = jf -4;
+ j1mo = j1mo - 2;
+ j1s = (j1mo + 1) * (j1mo + 1);
+ cj = sqrt(1.0 * (jmxpos - j1s) * (j1s - jmns));
+ j1mos = j1mo * j1mo;
+ cjmo = sqrt(1.0 * (jmxpos - j1mos) * (j1s - jmns));
+ c6->rac3j[irr - 2] = c6->rac3j[irr - 1] * (-cj / cjmo);
+ sjr = sjr + (c6->rac3j[irr - 1] * c6->rac3j[irr - 1]) * jf;
+ }
+ }
+ double racmat = c6->rac3j[nmat - 1];
+ sjr = sjr + (racmat * racmat) * (jf - 4);
+ c6->rac3j[0] = 1.0;
+ jf = jmn + jmn + 1;
+ double sjl = 1.0 * jf;
+ int j1pt = jmn + 2;
+ int j1pos = (j1pt - 1) * (j1pt - 1);
+ double cjpo = sqrt(1.0 * (jmxpos - j1pos) * (j1pos - jmns));
+ int j1pts = j1pt * j1pt;
+ double cjpt = sqrt(1.0 * (jmxpos - j1pts) * (j1pts - jmns));
+ c6->rac3j[1] = -cjpo / cjpt;
+ int nmatmo = nmat - 1;
+ if (nmatmo >= 2) {
+ for (int irl70 = 2; irl70 <= nmatmo; irl70++) {
+ jf = jf + 4;
+ j1pt = j1pt + 2;
+ j1pos = (j1pt - 1) * (j1pt - 1);
+ cjpo = sqrt(1.0 * (jmxpos - j1pos) * (j1pos - jmns));
+ j1pts = j1pt * j1pt;
+ cjpt = sqrt(1.0 * (jmxpos - j1pts) * (j1pts - jmns));
+ c6->rac3j[irl70] = c6->rac3j[irl70 - 1] * (-cjpo / cjpt);
+ sjl = sjl + (c6->rac3j[irl70 - 1] * c6->rac3j[irl70 - 1]) * jf;
+ }
+ }
+ double ratrac = racmat / c6->rac3j[nmat - 1];
+ double rats = ratrac * ratrac;
+ double sj = sjr + sjl * rats;
+ c6->rac3j[nmat - 1] = racmat;
+ double cnr = (1.0 / sqrt(sj)) * isn;
+ for (int irr80 = nmat; irr80 <= nj; irr80++) {
+ c6->rac3j[irr80 - 1] *= cnr;
+ }
+ double cnl = cnr * ratrac;
+ for (int irl85 = 1; irl85 <= nmatmo; irl85++) {
+ c6->rac3j[irl85 - 1] *= cnl;
+ }
+}
+
+/*! \brief C++ porting of STR
+ *
+ * \param c1: `C1 *` Pointer to a C1 instance.
+ * \param c1ao: `C1_AddOns *` Pointer to a C1_AddOns instance.
+ * \param c3: `C3 *` Pointer to a C3 instance.
+ * \param c4: `C4 *` Pointer to a C4 structure.
+ * \param c6: `C6 *` Pointer to a C6 instance.
+ */
+void str(C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
+ std::complex<double> *ylm;
+ c3 = new C3(); // this makes up for GCS = 0.0D0
+ int i = 0;
+ for (int l1po = 1; l1po <= c4->lmpo; l1po++) {
+ int l1 = l1po - 1;
+ for (int l2 = 1; l2 <= c4->lm; l2++) {
+ r3j000(l1, l2, c6);
+ c1ao->ind3j[l1po - 1][l2 - 1] = i;
+ int iabs = l2 - l1;
+ if (iabs < 0) iabs *= -1;
+ int lmnpo = iabs + 1;
+ int lmxpo = l2 + l1 + 1;
+ int il = 0;
+ int lpo = lmnpo;
+ while (lpo <= lmxpo) {
+ i++;
+ il++;
+ c1ao->v3j0[i - 1] = c6->rac3j[il - 1];
+ lpo += 2;
+ }
+ }
+ } // l1po loop
+ int nsphmo = c4->nsph - 1;
+ int lit = c4->li + c4-> li;
+ int lmt = c4->li + c4->le;
+ int size40 = nsphmo * c4->nsph * c4->litpos;
+ int size50 = c4->nsph * c4->lmtpos;
+ int ylm_size = (size40 > size50) ? size40 : size50;
+ ylm = new std::complex<double>[ylm_size];
+ int ivy = 0;
+ for (int nf40 = 1; nf40 <= nsphmo; nf40++) {
+ int nfpo = nf40 + 1;
+ for (int ns40 = nfpo; ns40 <= c4->nsph; ns40++) {
+ double rx = c1->rxx[nf40 - 1] - c1->rxx[ns40 - 1];
+ double ry = c1->ryy[nf40 - 1] - c1->ryy[ns40 - 1];
+ double rz = c1->rzz[nf40 - 1] - c1->rzz[ns40 - 1];
+ double rr, crth, srth, crph, srph;
+ polar(rx, ry, rz, rr, crth, srth, crph, srph);
+ //printf("DEBUG: crth = %lE\n", crth);
+ //printf("DEBUG: srth = %lE\n", srth);
+ //printf("DEBUG: crph = %lE\n", crph);
+ //printf("DEBUG: srph = %lE\n", srph);
+ sphar(crth, srth, crph, srph, lit, ylm);
+ for (int iv38 = 1; iv38 <= c4->litpos; iv38++) {
+ c1ao->vyhj[iv38 + ivy - 1] = dconjg(ylm[iv38 - 1]);
+ } // iv38 loop
+ ivy += c4->litpos;
+ } // ns40 loop
+ } // nf40 loop
+ ivy = 0;
+ for (int nf50 = 1; nf50 <= c4->nsph; nf50++) {
+ double rx = c1->rxx[nf50 - 1];
+ double ry = c1->ryy[nf50 - 1];
+ double rz = c1->rzz[nf50 - 1];
+ double rr, crth, srth, crph, srph;
+ if (rx == 0.0 && ry == 0.0 && rz == 0.0) continue;
+ polar(rx, ry, rz, rr, crth, srth, crph, srph);
+ sphar(crth, srth, crph, srph, lmt, ylm);
+ for (int iv48 = 1; iv48 <= c4->lmtpos; iv48++) {
+ c1ao->vyj0[iv48 + ivy - 1] = dconjg(ylm[iv48 - 1]);
+ } // iv48 loop
+ ivy += c4->lmtpos;
+ } // nf50 loop
+}
+
+#endif