diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
new file mode 100644
index 0000000000000000000000000000000000000000..97ea5d0873c1df715c069af076afe7503e1472e1
--- /dev/null
+++ b/src/include/sph_subs.h
@@ -0,0 +1,701 @@
+/*! \file sph_subs.h
+ *
+ * \brief C++ porting of SPH functions, subroutines and data structures.
+ *
+ * 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 SRC_INCLUDE_SPH_SUBS_H_
+#define SRC_INCLUDE_SPH_SUBS_H_
+
+#include <complex>
+
+//! \brief A structure representing the common block C1
+struct C1 {
+ std::complex<double> rmi[40][2];
+ std::complex<double> rei[40][2];
+ std::complex<double> w[3360][4];
+ std::complex<double> fsas[2];
+ std::complex<double> sas[2][2][2];
+ std::complex<double> vints[2][16];
+ double sscs[2];
+ double sexs[2];
+ double sabs[2];
+ double sqscs[2];
+ double sqexs[2];
+ double sqabs[2];
+ double gcsv[2];
+ double rxx[1];
+ double ryy[1];
+ double rzz[1];
+ double ros[2];
+ double rc[2][8];
+ int iog[2];
+ int nshl[2];
+};
+
+//! \brief A structure representing the common block C1
+struct C2 {
+ std::complex<double> ris[999];
+ std::complex<double> dlri[999];
+ std::complex<double> dc0[5];
+ std::complex<double> vkt[2];
+ double vsz[2];
+};
+
+//! \brief Write a message to test FORTRAN -> C++ communication
+extern "C" void testme_();
+
+/*! \brief C++ porting of DIEL
+ *
+ * \param npntmo: `int`
+ * \param ns: `int`
+ * \param i: `int`
+ * \param ic: `int`
+ * \param vk: `double`
+ * \param c1: `C1 *`
+ * \param c2: `C2 *`
+ */
+void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
+ double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1];
+ double half_step = 0.5 * dif / npntmo;
+ double rr = c1->rc[i - 1][ns - 1];
+ std::complex<double> delta = c2->dc0[ic] - c2->dc0[ic - 1];
+ int kpnt = npntmo + npntmo;
+ c2->ris[kpnt] = c2->dc0[ic];
+ c2->dlri[kpnt] = std::complex<double>(0.0, 0.0);
+ for (int np90 = 0; np90 < kpnt; np90++) {
+ double ff = (rr - c1->rc[i - 1][ns - 1]) / dif;
+ c2->ris[np90] = delta * ff * ff * (-2.0 * ff + 3.0) + c2->dc0[ic - 1];
+ c2->dlri[np90] = 3.0 * delta * ff * (1.0 - ff) / (dif * vk * c2->ris[np90]);
+ rr += half_step;
+ }
+}
+
+/*! \brief C++ porting of ENVJ
+ *
+ * \param n: `int`
+ * \param x: `double`
+ * \return result: `double`
+ */
+double envj(int n, double x) {
+ double result = 0.0;
+ double xn;
+ if (n == 0) {
+ xn = 1.0e-100;
+ result = 0.5 * log10(6.28 * xn) - xn * log10(1.36 * x / xn);
+ } else {
+ result = 0.5 * log10(6.28 * n) - n * log10(1.36 * x / n);
+ }
+ return result;
+}
+
+/*! \brief C++ porting of MSTA1
+ *
+ * \param x: `double`
+ * \param mp: `int`
+ * \return result: `int`
+ */
+int msta1(double x, int mp) {
+ int result = 0;
+ double a0 = x;
+ if (a0 < 0.0) a0 *= -1.0;
+ int n0 = (int)(1.1 * a0) + 1;
+ double f0 = envj(n0, a0) - mp;
+ int n1 = n0 + 5;
+ double f1 = envj(n1, a0) - mp;
+ for (int it10 = 0; it10 < 20; it10++) {
+ int nn = n1 - (int)((n1 - n0) / (1.0 - f0 / f1));
+ double f = envj(nn, a0) - mp;
+ int test_n = nn - n1;
+ if (test_n < 0) test_n *= -1;
+ if (test_n < 1) {
+ return nn;
+ }
+ n0 = n1;
+ f0 = f1;
+ n1 = nn;
+ f1 = f;
+ result = nn;
+ }
+ return result;
+}
+
+/*! \brief C++ porting of MSTA2
+ *
+ * \param x: `double`
+ * \param n: `int`
+ * \param mp: `int`
+ * \return result: `int`
+ */
+int msta2(double x, int n, int mp) {
+ int result = 0;
+ double a0 = x;
+ if (a0 < 0) a0 *= -1.0;
+ double half_mp = 0.5 * mp;
+ double ejn = envj(n, a0);
+ double obj;
+ int n0;
+ if (ejn <= half_mp) {
+ obj = 1.0 * mp;
+ n0 = (int)(1.1 * a0) + 1;
+ } else {
+ obj = half_mp + ejn;
+ n0 = n;
+ }
+ double f0 = envj(n0, a0) - obj;
+ int n1 = n0 + 5;
+ double f1 = envj(n1, a0) - obj;
+ for (int it10 = 0; it10 < 20; it10 ++) {
+ int nn = n1 - (int)((n1 - n0) / (1.0 - f0 / f1));
+ double f = envj(nn, a0) - obj;
+ int test_n = nn - n1;
+ if (test_n < 0) test_n *= -1;
+ if (test_n < 1) return (nn + 10);
+ n0 = n1;
+ f0 = f1;
+ n1 = nn;
+ f1 = f;
+ result = nn + 10;
+ }
+ return result;
+}
+
+/*! \brief C++ porting of CBF
+ *
+ * This is the Complex Bessel Function.
+ *
+ * \param n: `int`
+ * \param z: `complex\<double\>`
+ * \param nm: `int &`
+ * \param csj: Vector of complex.
+ */
+void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
+ /*
+ * FROM CSPHJY OF LIBRARY specfun
+ *
+ * ==========================================================
+ * Purpose: Compute spherical Bessel functions j
+ * Input : z --- Complex argument of j
+ * n --- Order of j ( n = 0,1,2,... )
+ * Output: csj(n+1) --- j
+ * nm --- Highest order computed
+ * Routines called:
+ * msta1 and msta2 for computing the starting
+ * point for backward recurrence
+ * ==========================================================
+ */
+ double zz = z.real() * z.real() + z.imag() * z.imag();
+ double a0 = sqrt(zz);
+ nm = n;
+ if (a0 < 1.0e-60) {
+ for (int k = 2; k <= n + 1; k++) {
+ csj[k - 1] = 0.0;
+ }
+ csj[0] = std::complex<double>(1.0, 0.0);
+ return;
+ }
+ csj[0] = std::sin(z) / z;
+ if (n == 0) {
+ return;
+ }
+ csj[1] = (csj[0] -std::cos(z)) / z;
+ if (n == 1) {
+ return;
+ }
+ std::complex<double> csa = csj[0];
+ std::complex<double> csb = csj[1];
+ int m = msta1(a0, 200);
+ if (m < n) nm = m;
+ else m = msta2(a0, n, 15);
+ std::complex<double> cf0 = 0.0;
+ std::complex<double> cf1 = 1.0e-100;
+ std::complex<double> cf, cs;
+ for (int k = m; k >= 0; k--) {
+ cf = (2.0 * k + 3.0) * cf1 / z - cf0;
+ if (k <= nm) csj[k] = cf;
+ cf0 = cf1;
+ cf1 = cf;
+ }
+ double abs_csa = sqrt(csa.real() * csa.real() + csa.imag() * csa.imag());
+ double abs_csb = sqrt(csb.real() * csb.real() + csb.imag() * csb.imag());
+ if (abs_csa > abs_csb) cs = csa / cf;
+ else cs = csb / cf0;
+ for (int k = 0; k <= nm; k++) {
+ csj[k] = cs * csj[k];
+ }
+}
+
+/*! \brief C++ porting of RBF
+ *
+ * This is the Real Bessel Function.
+ *
+ * \param n: `int`
+ * \param x: `double`
+ * \param nm: `int &`
+ * \param sj: `double[]`
+ */
+void rbf(int n, double x, int &nm, double sj[]) {
+ /*
+ * FROM SPHJ OF LIBRARY specfun
+ *
+ * ==========================================================
+ * Purpose: Compute spherical Bessel functions j
+ * Input : x --- Argument of j
+ * n --- Order of j ( n = 0,1,2,... )
+ * Output: sj(n+1) --- j
+ * nm --- Highest order computed
+ * Routines called:
+ * msta1 and msta2 for computing the starting
+ * point for backward recurrence
+ * ==========================================================
+ */
+ double a0 = x;
+ if (a0 < 0.0) a0 *= -1.0;
+ nm = n;
+ if (a0 < 1.0e-60) {
+ for (int k = 2; k <= n + 1; k++)
+ sj[k - 1] = 0.0;
+ sj[0] = 1.0;
+ return;
+ }
+ sj[0] = sin(x) / x;
+ if (n == 0) {
+ return;
+ }
+ sj[1] = (sj[0] - cos(x)) / x;
+ if (n == 1) {
+ return;
+ }
+ double sa = sj[0];
+ double sb = sj[1];
+ int m = msta1(a0, 200);
+ if (m < n) nm = m;
+ else m = msta2(a0, n, 15);
+ double f0 = 0.0;
+ double f1 = 1.0e-100;
+ double f;
+ for (int k = m; k >= 0; k--) {
+ f = (2.0 * k +3.0) * f1 / x - f0;
+ if (k <= nm) sj[k] = f;
+ f0 = f1;
+ f1 = f;
+ }
+ double cs;
+ double abs_sa = sa;
+ if (abs_sa < 0.0) abs_sa *= -1.0;
+ double abs_sb = sb;
+ if (abs_sb < 0.0) abs_sb *= -1.0;
+ if (abs_sa > abs_sb) cs = sa / f;
+ else cs = sb / f0;
+ for (int k = 0; k <= nm; k++) {
+ sj[k] = cs * sj[k];
+ }
+}
+
+/*! \brief C++ porting of RKC
+ *
+ * \param npntmo: `int`
+ * \param step: `double`
+ * \param dcc: `complex\<double\>`
+ * \param x: `double &`
+ * \param lpo: `int`
+ * \param y1: `complex\<double\> &`
+ * \param y2: `complex\<double\> &`
+ * \param dy1: `complex\<double\> &`
+ * \param dy2: `complex\<double\> &`
+ */
+void rkc(
+ int npntmo, double step, std::complex<double> dcc, double &x, int lpo,
+ std::complex<double> &y1, std::complex<double> &y2, std::complex<double> &dy1,
+ std::complex<double> &dy2
+) {
+ std::complex<double> cy1, cdy1, c11, cy23, yc2, c12, c13;
+ std::complex<double> cy4, yy, c14, c21, c22, c23, c24;
+ double half_step = 0.5 * step;
+ double cl = 1.0 * lpo * (lpo - 1);
+ for (int ipnt60 = 1; ipnt60 <= npntmo; ipnt60++) {
+ cy1 = cl / (x * x) - dcc;
+ cdy1 = -2.0 / x;
+ c11 = (cy1 * y1 + cdy1 * dy1) * step;
+ double xh = x + half_step;
+ cy23 = cl / (xh * xh) - dcc;
+ double cdy23 = -2.0 / xh;
+ yc2 = y1 + dy1 * half_step;
+ c12 = (cy23 * yc2 + cdy23 * (dy1 + 0.5 * c11)) * step;
+ c13 = (cy23 * (yc2 + 0.25 * c11 * step) + cdy23 * (dy1 + 0.5 * c12)) * step;
+ double xn = x + step;
+ cy4 = cl / (xn * xn) - dcc;
+ double cdy4 = -2.0 / xn;
+ yy = y1 + dy1 * step;
+ c14 = (cy4 * (yy + 0.5 * c12 * step) + cdy4 * (dy1 + c13)) * step;
+ y1 = yy + (c11 + c12 + c13) * step / 6.0;
+ dy1 += (0.5 * c11 + c12 + c13 + 0.5 * c14) / 3.0;
+ c21 = (cy1 * y2 + cdy1 * dy2) * step;
+ yc2 = y2 + dy2 * half_step;
+ c22 = (cy23 * yc2 + cdy23 * (dy2 + 0.5 * c21)) * step;
+ c23= (cy23 * (yc2 + 0.25 * c21 * step) + cdy23 * (dy2 + 0.5 * c22)) * step;
+ yy = y2 + dy2 * step;
+ c24 = (cy4 * (yc2 + 0.5 * c22 * step) + cdy4 * (dy2 + c23)) * step;
+ y2 = yy + (c21 + c22 + c23) * step / 6.0;
+ dy2 += (0.5 * c21 + c22 + c23 + 0.5 * c24) / 3.0;
+ x = xn;
+ }
+}
+
+/*! \brief C++ porting of RKT
+ *
+ * \param npntmo: `int`
+ * \param step: `double`
+ * \param x: `double &`
+ * \param lpo: `int`
+ * \param y1: `complex\<double\> &`
+ * \param y2: `complex\<double\> &`
+ * \param dy1: `complex\<double\> &`
+ * \param dy2: `complex\<double\> &`
+ * \param c2: `C2 *`
+ */
+void rkt(
+ int npntmo, double step, double &x, int lpo, std::complex<double> &y1,
+ std::complex<double> &y2, std::complex<double> &dy1, std::complex<double> &dy2,
+ C2 *c2
+) {
+ std::complex<double> cy1, cdy1, c11, cy23, cdy23, yc2, c12, c13;
+ std::complex<double> cy4, cdy4, yy, c14, c21, c22, c23, c24;
+ double half_step = 0.5 * step;
+ double cl = 1.0 * lpo * (lpo - 1);
+ for (int ipnt60 = 1; ipnt60 <= npntmo; ipnt60++) {
+ int jpnt = ipnt60 + ipnt60 - 1;
+ cy1 = cl / (x * x) - c2->ris[jpnt - 1];
+ cdy1 = -2.0 / x;
+ c11 = (cy1 * y1 + cdy1 * dy1) * step;
+ double xh = x + half_step;
+ int jpntpo = jpnt + 1;
+ cy23 = cl / (xh * xh) - c2->ris[jpntpo - 1];
+ cdy23 = -2.0 / xh;
+ yc2 = y1 + dy1 * half_step;
+ c12 = (cy23 * yc2 + cdy23 * (dy1 + 0.5 * c11)) * step;
+ c13= (cy23 * (yc2 + 0.25 * c11 *step) + cdy23 * (dy1 + 0.5 * c12)) * step;
+ double xn = x + step;
+ int jpntpt = jpnt + 2;
+ cy4 = cl / (xn * xn) - c2->ris[jpntpt - 1];
+ cdy4 = -2.0 / xn;
+ yy = y1 + dy1 * step;
+ c14 = (cy4 * (yy + 0.5 * c12 * step) + cdy4 * (dy1 + c13)) * step;
+ y1= yy + (c11 + c12 + c13) * step / 6.0;
+ dy1 += (0.5 * c11 + c12 + c13 + 0.5 *c14) /3.0;
+ cy1 -= cdy1 * c2->dlri[jpnt - 1];
+ cdy1 += 2.0 * c2->dlri[jpnt - 1];
+ c21 = (cy1 * y2 + cdy1 * dy2) * step;
+ cy23 -= cdy23 * c2->dlri[jpntpo - 1];
+ cdy23 += 2.0 * c2->dlri[jpntpo - 1];
+ yc2 = y2 + dy2 * half_step;
+ c22 = (cy23 * yc2 + cdy23 * (dy2 + 0.5 * c21)) * step;
+ c23 = (cy23 * (yc2 + 0.25 * c21 * step) + cdy23 * (dy2 + 0.5 * c22)) * step;
+ cy4 -= cdy4 * c2->dlri[jpntpt - 1];
+ cdy4 += 2.0 * c2->dlri[jpntpt - 1];
+ yy = y2 + dy2 * step;
+ c24 = (cy4 * (yc2 + 0.5 * c22 * step) + cdy4 * (dy2 + c23)) * step;
+ y2 = yy + (c21 + c22 + c23) * step / 6.0;
+ dy2 += (0.5 * c21 + c22 + c23 + 0.5 * c24) / 3.0;
+ x = xn;
+ }
+}
+
+/*! \brief C++ porting of RNF
+ *
+ * This is a real spherical Bessel function.
+ *
+ * \param n: `int`
+ * \param x: `double`
+ * \param nm: `int &`
+ * \param sj: `double[]`
+ */
+void rnf(int n, double x, int &nm, double sy[]) {
+ /*
+ * FROM SPHJY OF LIBRARY specfun
+ *
+ * ==========================================================
+ * Purpose: Compute spherical Bessel functions y
+ * Input : x --- Argument of y ( x > 0 )
+ * n --- Order of y ( n = 0,1,2,... )
+ * Output: sy(n+1) --- y
+ * nm --- Highest order computed
+ * ==========================================================
+ */
+ if (x < 1.0e-60) {
+ for (int k = 0; k <= n; k++)
+ sy[k] = -1.0e300;
+ return;
+ }
+ sy[0] = -1.0 * cos(x) / x;
+ if (n == 0) {
+ return;
+ }
+ sy[1] = (sy[0] - sin(x)) / x;
+ if (n == 1) {
+ return;
+ }
+ double f0 = sy[0];
+ double f1 = sy[1];
+ double f;
+ for (int k = 2; k <= n; k++) {
+ f = (2.0 * k - 1.0) * f1 / x - f0;
+ sy[k] = f;
+ double abs_f = f;
+ if (abs_f < 0.0) abs_f *= -1.0;
+ if (abs_f >= 1.0e300) {
+ nm = k;
+ break;
+ }
+ f0 = f1;
+ f1 = f;
+ nm = k;
+ }
+ return;
+}
+
+/*! \brief C++ porting of SSCR0
+ *
+ * \param tfsas: `complex<double> &`
+ * \param nsph: `int`
+ * \param lm: `int`
+ * \param vk: `double`
+ * \param exri: `double`
+ * \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);
+ 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);
+ tfsas = cc0;
+ for (int i12 = 1; i12 <= nsph; i12++) {
+ int iogi = c1->iog[i12 - 1];
+ if (iogi >= i12) {
+ double sums = 0.0;
+ std::complex<double> sum21 = cc0;
+ for (int l10 = 1; l10 <= lm; l10++) {
+ double fl = 1.0 + l10 + l10;
+ rm = 1.0 / c1->rmi[l10 - 1][i12 - 1];
+ re = 1.0 / c1->rei[l10 - 1][i12 - 1];
+ std::complex<double> rm_cnjg = std::complex<double>(rm.real(), -rm.imag());
+ std::complex<double> re_cnjg = std::complex<double>(re.real(), -re.imag());
+ sums += (rm_cnjg * rm + re_cnjg * re).real() * fl;
+ sum21 += (rm + re) * fl;
+ }
+ sum21 *= -1.0;
+ double scasec = cccs * sums;
+ double extsec = -cccs * sum21.real();
+ double abssec = extsec - scasec;
+ c1->sscs[i12 - 1] = scasec;
+ c1->sexs[i12 - 1] = extsec;
+ c1->sabs[i12 - 1] = abssec;
+ double gcss = c1-> gcsv[i12 - 1];
+ c1->sqscs[i12 - 1] = scasec / gcss;
+ c1->sqexs[i12 - 1] = extsec / gcss;
+ c1->sqabs[i12 - 1] = abssec / gcss;
+ c1->fsas[i12 - 1] = sum21 * csam;
+ }
+ tfsas += c1->fsas[iogi - 1];
+ }
+}
+
+/*! \brief C++ porting of THDPS
+ *
+ * \param lm: `int`
+ * \param zpv: `double ****`
+ */
+void thdps(int lm, double ****zpv) {
+ // WARNING: unclear nested loop in THDPS
+ // The optimized interpretation was adopted here.
+ for (int l10 = 1; l10 <= lm; l10++) {
+ for (int ilmp = 1; ilmp <= 3; ilmp++) {
+ zpv[l10 - 1][ilmp - 1][0][0] = 0.0;
+ zpv[l10 - 1][ilmp - 1][0][1] = 0.0;
+ zpv[l10 - 1][ilmp - 1][1][0] = 0.0;
+ zpv[l10 - 1][ilmp - 1][1][1] = 0.0;
+ }
+ }
+ for (int l15 = 1; l15 <= lm; l15++) {
+ double xd = 1.0 * l15 * (l15 + 1);
+ double zp = -1.0 / sqrt(xd);
+ zpv[l15 - 1][1][0][1] = zp;
+ zpv[l15 - 1][1][1][0] = zp;
+ }
+ if (lm != 1) {
+ for (int l20 = 2; l20 <= lm; l20++) {
+ double xn = 1.0 * (l20 - 1) * (l20 + 1);
+ double xd = 1.0 * l20 * (l20 + l20 + 1);
+ double zp = sqrt(xn / xd);
+ zpv[l20 - 1][0][0][0] = zp;
+ zpv[l20 - 1][0][1][1] = zp;
+ }
+ int lmmo = lm - 1;
+ for (int l25 = 1; l25 <= lmmo; l25++) {
+ double xn = 1.0 * l25 * (l25 + 2);
+ double xd = (l25 + 1) * (l25 + l25 + 1);
+ double zp = -1.0 * sqrt(xn / xd);
+ zpv[l25 - 1][2][0][0] = zp;
+ zpv[l25 - 1][2][1][1] = zp;
+ }
+ }
+}
+
+/*! \brief C++ porting of DME
+ *
+ * \param li: `int`
+ * \param i: `int`
+ * \param npnt: `int`
+ * \param npntts: `int`
+ * \param vk: `double`
+ * \param exdc: `double`
+ * \param exri: `double`
+ * \param c1: `C1 *`
+ * \param c2: `C2 *`
+ * \param jer: `int &`
+ * \param lcalc: `int &`
+ * \param arg: `complex<double> &`.
+ */
+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) {
+ //double rfj[42], rfn[42];
+ double *rfj = new double[42];
+ double *rfn = new double[42];
+ std::complex<double> cfj[42], fbi[42], fb[42], fn[42];
+ std::complex<double> rmf[40], drmf[40], ref[40], dref[40];
+ std::complex<double> dfbi, dfb, dfn, ccna, ccnb, ccnc, ccnd;
+ std::complex<double> y1, dy1, y2, dy2, arin, cri, uim;
+ jer = 0;
+ uim = std::complex<double>(0.0, 1.0);
+ int nstp = npnt - 1;
+ int nstpts = npntts - 1;
+ int lipo = li + 1;
+ int lipt = li + 2;
+ double sz = vk * c1->ros[i - 1];
+ c2->vsz[i - 1] = sz;
+ double vkr1 = vk * c1->rc[i - 1][0];
+ int nsh = c1->nshl[i - 1];
+ c2->vkt[i - 1] = std::sqrt(c2->dc0[0]);
+ arg = vkr1 * c2->vkt[i - 1];
+ arin = arg;
+ bool goto32 = false;
+ if (arg.imag() != 0.0) {
+ cbf(lipo, arg, lcalc, cfj);
+ if (lcalc < lipo) {
+ jer = 5;
+ return;
+ }
+ for (int j24 = 1; j24 <= lipt; j24++) fbi[j24 - 1] = cfj[j24 - 1];
+ goto32 = true;
+ }
+ if (!goto32) {
+ rbf(lipo, arg.real(), lcalc, rfj);
+ if (lcalc < lipo) {
+ jer = 5;
+ return;
+ }
+ for (int j30 = 1; j30 <= lipt; j30++) fbi[j30 - 1] = rfj[j30 - 1];
+ }
+ double arex = sz * exri;
+ arg = arex;
+ rbf(lipo, arex, lcalc, rfj);
+ if (lcalc < lipo) {
+ jer = 7;
+ return;
+ }
+ rnf(lipo, arex, lcalc, rfn);
+ if (lcalc < lipo) {
+ jer = 8;
+ return;
+ }
+ for (int j43 = 1; j43 <= lipt; j43++) {
+ fb[j43 - 1] = rfj[j43 - 1];
+ //printf("DEBUG: fb[%d] = (%lE,%lE)\n", j43, fb[j43 - 1].real(), fb[j43 - 1].imag());
+ fn[j43 - 1] = rfn[j43 - 1];
+ //printf("DEBUG: fn[%d] = (%lE,%lE)\n", j43, fn[j43 - 1].real(), fn[j43 - 1].imag());
+ }
+ if (nsh <= 1) {
+ cri = c2->dc0[0] / exdc;
+ for (int l60 = 1; l60 <= li; l60++) {
+ int lpo = l60 + 1;
+ int ltpo = lpo + l60;
+ int lpt = lpo + 1;
+ dfbi = ((1.0 * l60) * fbi[l60 - 1] - (1.0 * lpo) * fbi[lpt - 1]) * arin + fbi[lpo - 1] * (1.0 * ltpo);
+ dfb = ((1.0 * l60) * fb[l60 - 1] - (1.0 * lpo) * fb[lpt - 1]) * arex + fb[lpo - 1] * (1.0 * ltpo);
+ dfn = ((1.0 * l60) * fn[l60 - 1] - (1.0 * lpo) * fn[lpt - 1]) * arex + fn[lpo - 1] * (1.0 * ltpo);
+ ccna = fbi[lpo - 1] * dfn;
+ ccnb = fn[lpo - 1] * dfbi;
+ ccnc = fbi[lpo - 1] * dfb;
+ ccnd = fb[lpo - 1] * dfbi;
+ c1->rmi[l60 - 1][i - 1] = 1.0 + uim * (ccna - ccnb) / (ccnc - ccnd);
+ c1->rei[l60 - 1][i - 1] = 1.0 + uim * (cri * ccna - ccnb) / (cri * ccnc - ccnd);
+ }
+ } else { // nsh > 1
+ int ic = 1;
+ for (int l80 = 1; l80 <= li; l80++) {
+ int lpo = l80 + 1;
+ int ltpo = lpo + l80;
+ int lpt = lpo + 1;
+ int dltpo = ltpo;
+ y1 = fbi[lpo - 1];
+ dy1 = ((1.0 * l80) * fbi[l80 - 1] - (1.0 * lpo) * fbi[lpt - 1]) * c2->vkt[i - 1] / (1.0 * dltpo);
+ y2 = y1;
+ dy2 = dy1;
+ ic = 1;
+ for (int ns76 = 2; ns76 <= nsh; ns76++) {
+ int nsmo = ns76 - 1;
+ double vkr = vk * c1->rc[i - 1][nsmo - 1];
+ if (ns76 % 2 != 0) {
+ ic += 1;
+ double step = 1.0 * nstp;
+ step = vk * (c1->rc[i - 1][ns76 - 1] - c1->rc[i - 1][nsmo - 1]) / step;
+ arg = c2->dc0[ic - 1];
+ rkc(nstp, step, arg, vkr, lpo, y1, y2, dy1, dy2);
+ } else {
+ diel(nstpts, nsmo, i, ic, vk, c1, c2);
+ double stepts = 1.0 * nstpts;
+ stepts = vk * (c1->rc[i - 1][ns76 - 1] - c1->rc[i - 1][nsmo - 1]) / stepts;
+ rkt(nstpts, stepts, vkr, lpo, y1, y2, dy1, dy2, c2);
+ }
+ }
+ rmf[l80 - 1] = y1 * sz;
+ drmf[l80 - 1] = dy1 * sz + y1;
+ ref[l80 - 1] = y2 * sz;
+ dref[l80 - 1] = dy2 * sz + y2;
+ }
+ cri = 1.0 + uim * 0.0;
+ if (nsh % 2 != 0) cri = c2->dc0[ic - 1] / exdc;
+ for (int l90 = 1; l90 <= li; l90++) {
+ int lpo = l90 + 1;
+ int ltpo = lpo + l90;
+ int lpt = lpo + 1;
+ dfb = ((1.0 * l90) * fb[l90 - 1] - (1.0 * lpo) * fb[lpt - 1]) * arex + fb[lpo - 1] * (1.0 * ltpo);
+ dfn = ((1.0 * l90) * fn[l90 - 1] - (1.0 * lpo) * fn[lpt - 1]) * arex + fn[lpo - 1] * (1.0 * ltpo);
+ ccna = rmf[l90 - 1] * dfn;
+ ccnb = drmf[l90 - 1] * fn[lpo - 1] * (1.0 * sz * ltpo);
+ ccnc = rmf[l90 - 1] * dfb;
+ ccnd = drmf[l90 - 1] * fb[lpo -1] * (1.0 * sz * ltpo);
+ c1->rmi[l90 - 1][i - 1] = 1.0 + uim *(ccna - ccnb) / (ccnc - ccnd);
+ //printf("DEBUG: gone 90, rmi[%d][%d] = (%lE,%lE)\n", l90, i, c1->rmi[l90 - 1][i - 1].real(), c1->rmi[l90 - 1][i - 1].imag());
+ ccna = ref[l90 - 1] * dfn;
+ ccnb = dref[l90 - 1] * fn[lpo - 1] * (1.0 * sz * ltpo);
+ ccnc = ref[l90 - 1] * dfb;
+ ccnd = dref[l90 - 1] *fb[lpo - 1] * (1.0 * sz * ltpo);
+ c1->rei[l90 - 1][i - 1] = 1.0 + uim * (cri * ccna - ccnb) / (cri * ccnc - ccnd);
+ //printf("DEBUG: gone 90, rei[%d][%d] = (%lE,%lE)\n", l90, i, c1->rei[l90 - 1][i - 1].real(), c1->rei[l90 - 1][i - 1].imag());
+ }
+ } // nsh <= 1 ?
+ return;
+}
+
+#endif /* SRC_INCLUDE_SPH_SUBS_H_ */