Create sphere subroutines module

/*! \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);