From 2f144ba92e54cc97e1952b8f6e2696e12dec2b0c Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Tue, 9 Jan 2024 10:36:41 +0100
Subject: [PATCH] Move SPH subroutine implementation to sph_subs.cpp
---
src/cluster/Makefile | 6 +-
src/include/sph_subs.h | 1283 +++-----------------------------------
src/libnptm/sph_subs.cpp | 1167 ++++++++++++++++++++++++++++++++++
src/sphere/Makefile | 7 +-
4 files changed, 1275 insertions(+), 1188 deletions(-)
create mode 100644 src/libnptm/sph_subs.cpp
diff --git a/src/cluster/Makefile b/src/cluster/Makefile
index 694bbc0d..71a106ed 100644
--- a/src/cluster/Makefile
+++ b/src/cluster/Makefile
@@ -15,7 +15,7 @@ clu: clu.o
edfb: edfb.o
$(FC) $(FCFLAGS) -o $(BUILDDIR)/edfb $(BUILDDIR)/edfb.o $(LFLAGS)
-np_cluster: $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o $(BUILDDIR)/cluster.o
+np_cluster: $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/cluster.o
$(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_cluster $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o $(BUILDDIR)/cluster.o
$(BUILDDIR)/np_cluster.o:
@@ -33,8 +33,8 @@ $(BUILDDIR)/Parsers.o:
$(BUILDDIR)/cluster.o:
$(CXX) $(CXXFLAGS) -c cluster.cpp -o $(BUILDDIR)/cluster.o
-$(BUILDDIR)/sphere.o:
- $(CXX) $(CXXFLAGS) -c ../sphere/sphere.cpp -o $(BUILDDIR)/sphere.o
+$(BUILDDIR)/sph_subs.o:
+ $(CXX) $(CXXFLAGS) -c ../libnptm/sph_subs.cpp -o $(BUILDDIR)/sphere.o
clean:
rm -f $(BUILDDIR)/*.o
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index accc3ea0..781f6795 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -20,16 +20,32 @@
#include <complex>
-/*! \brief Conjugate of a double precision complex number
+/*! \brief Compute the asymmetry-corrected scattering cross-section.
*
- * \param z: `complex<double>` The input complex number.
- * \return result: `complex<double>` The conjugate of the input number.
+ * This function computes the product between the geometrical asymmetry parameter and
+ * the scattering cross-section. See Sec. 3.2.1 of Borghese, Denti & Saija (2007).
+ *
+ * \param zpv: `double ****` Geometrical asymmetry parameter coefficients matrix.
+ * \param li: `int` Maximum field expansion order.
+ * \param nsph: `int` Number of spheres.
+ * \param c1: `C1 *` Pointer to a `C1` data structure.
+ * \param sqk: `double`
+ * \param gaps: `double *` Geometrical asymmetry parameter-corrected cross-section.
+ */
+void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps);
+
+/*! \brief Complex Bessel Function.
+ *
+ * This function computes the complex spherical Bessel funtions \f$j\f$. It uses the
+ * auxiliary functions `msta1()` and `msta2()` to determine the starting point for
+ * backward recurrence. This is the `CSPHJ` implementation of the `specfun` library.
+ *
+ * \param n: `int` Order of the function.
+ * \param z: `complex<double>` Argumento of the function.
+ * \param nm: `int &` Highest computed order.
+ * \param csj: Vector of complex. The desired function \f$j\f$.
*/
-std::complex<double> dconjg(std::complex<double> z) {
- double zreal = z.real();
- double zimag = z.imag();
- return std::complex<double>(zreal, -zimag);
-}
+void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]);
/*! \brief Clebsch-Gordan coefficients.
*
@@ -42,124 +58,14 @@ std::complex<double> dconjg(std::complex<double> z) {
* \param m: `int`
* \return result: `double` Clebsh-Gordan coefficient.
*/
-double cg1(int lmpml, int mu, int l, int m) {
- double result = 0.0;
- double xd, xn;
- if (lmpml == -1) { // Interpreted as GOTO 30
- xd = 2.0 * l * (2 * l - 1);
- if (mu == -1) {
- xn = 1.0 * (l - 1 - m) * (l - m);
- } else if (mu == 0) {
- xn = 2.0 * (l - m) * (l + m);
- } else if (mu == 1) {
- xn = 1.0 * (l - 1 + m) * (l + m);
- } else {
- throw 111; // Need an exception for unpredicted indices.
- }
- result = sqrt(xn / xd);
- } else if (lmpml == 0) { // Interpreted as GOTO 5
- bool goto10 = (m != 0) || (mu != 0);
- if (!goto10) {
- result = 0.0;
- return result;
- }
- if (mu != 0) {
- xd = 2.0 * l * (l + 1);
- if (mu == -1) {
- xn = 1.0 * (l - m) * (l + m + 1);
- result = -sqrt(xn / xd);
- } else if (mu == 1) { // mu > 0
- xn = 1.0 * (l + m) * (l - m + 1);
- result = sqrt(xn / xd);
- } else {
- throw 111; // Need an exception for unpredicted indices.
- }
- } else { // mu = 0
- xd = 1.0 * l * (l + 1);
- xn = -1.0 * m;
- result = xn / sqrt(xd);
- }
- } else if (lmpml == 1) { // Interpreted as GOTO 60
- xd = 2.0 * (l * 2 + 3) * (l + 1);
- if (mu == -1) {
- xn = 1.0 * (l + 1 + m) * (l + 2 + m);
- result = sqrt(xn / xd);
- } else if (mu == 0) {
- xn = 2.0 * (l + 1 - m) * (l + 1 + m);
- result = -sqrt(xn / xd);
- } else if (mu == 1) {
- xn = 1.0 * (l + 1 - m) * (l + 2 - m);
- result = sqrt(xn / xd);
- } else { // mu was not recognized.
- throw 111; // Need an exception for unpredicted indices.
- }
- } else { // lmpml was not recognized
- throw 111; // Need an exception for unpredicted indices.
- }
- return result;
-}
+double cg1(int lmpml, int mu, int l, int m);
-/*! \brief Compute the asymmetry-corrected scattering cross-section.
- *
- * This function computes the product between the geometrical asymmetry parameter and
- * the scattering cross-section. See Sec. 3.2.1 of Borghese, Denti & Saija (2007).
+/*! \brief Conjugate of a double precision complex number
*
- * \param zpv: `double ****` Geometrical asymmetry parameter coefficients matrix.
- * \param li: `int` Maximum field expansion order.
- * \param nsph: `int` Number of spheres.
- * \param c1: `C1 *` Pointer to a `C1` data structure.
- * \param sqk: `double`
- * \param gaps: `double *` Geometrical asymmetry parameter-corrected cross-section.
+ * \param z: `complex<double>` The input complex number.
+ * \return result: `complex<double>` The conjugate of the input number.
*/
-void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
- std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
- std::complex<double> summ, sume, suem, suee, sum;
- double half_pi = acos(0.0);
- double cofs = half_pi * 2.0 / sqk;
- for (int i40 = 0; i40 < nsph; i40++) {
- int i = i40 + 1;
- int iogi = c1->iog[i40];
- if (iogi >= i) {
- sum = cc0;
- for (int l30 = 0; l30 < li; l30++) {
- int l = l30 + 1;
- int ltpo = l + l + 1;
- for (int ilmp = 1; ilmp < 4; ilmp++) {
- int ilmp30 = ilmp - 1;
- bool goto30 = (l == 1 && ilmp == 1) || (l == li && ilmp == 3);
- if (!goto30) {
- int lmpml = ilmp - 2;
- int lmp = l + lmpml;
- double cofl = sqrt(1.0 * (ltpo * (lmp + lmp + 1)));
- summ = zpv[l30][ilmp30][0][0] /
- (
- dconjg(c1->rmi[l30][i40]) *
- c1->rmi[lmp - 1][i40]
- );
- sume = zpv[l30][ilmp30][0][1] /
- (
- dconjg(c1->rmi[l30][i40]) *
- c1->rei[lmp - 1][i40]
- );
- suem = zpv[l30][ilmp30][1][0] /
- (
- dconjg(c1->rei[l30][i40]) *
- c1->rmi[lmp - 1][i40]
- );
- suee = zpv[l30][ilmp30][1][1] /
- (
- dconjg(c1->rei[l30][i40]) *
- c1->rei[lmp - 1][i40]
- );
- sum += (cg1(lmpml, 0, l, -1) * (summ - sume - suem + suee) +
- cg1(lmpml, 0, l, 1) * (summ + sume + suem + suee)) * cofl;
- }
- }
- }
- }
- gaps[i40] = sum.real() * cofs;
- }
-}
+std::complex<double> dconjg(std::complex<double> z);
/*! \brief Comute the continuous variation of the refractive index and of its radial derivative.
*
@@ -175,24 +81,32 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
* \param c1: `C1 *` Pointer to `C1` data structure.
* \param c2: `C2 *` Pointer to `C2` data structure.
*/
-void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
- const double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1];
- const double half_step = 0.5 * dif / npntmo;
- double rr = c1->rc[i - 1][ns - 1];
- const std::complex<double> delta = c2->dc0[ic] - c2->dc0[ic - 1];
- const int kpnt = npntmo + npntmo;
- c2->ris[kpnt] = c2->dc0[ic];
- c2->dlri[kpnt] = std::complex<double>(0.0, 0.0);
- const int i90 = i - 1;
- const int ns90 = ns - 1;
- const int ic90 = ic - 1;
- for (int np90 = 0; np90 < kpnt; np90++) {
- double ff = (rr - c1->rc[i90][ns90]) / dif;
- c2->ris[np90] = delta * ff * ff * (-2.0 * ff + 3.0) + c2->dc0[ic90];
- c2->dlri[np90] = 3.0 * delta * ff * (1.0 - ff) / (dif * vk * c2->ris[np90]);
- rr += half_step;
- }
-}
+void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2);
+
+/*! \brief Compute Mie scattering coefficients.
+ *
+ * This function determines the L-dependent Mie scattering coefficients \f$a_l\f$ and \f$b_l\f$
+ * for the cases of homogeneous spheres, radially non-homogeneous spheres and, in case of sphere
+ * with dielectric function, sustaining longitudinal waves. See Sec. 5.1 in Borghese, Denti
+ * & Saija (2007).
+ *
+ * \param li: `int` Maximum field expansion order.
+ * \param i: `int`
+ * \param npnt: `int`
+ * \param npntts: `int`
+ * \param vk: `double` Wave number in scale units.
+ * \param exdc: `double` External medium dielectric constant.
+ * \param exri: `double` External medium refractive index.
+ * \param c1: `C1 *` Pointer to a `C1` data structure.
+ * \param c2: `C2 *` Pointer to a `C2` data structure.
+ * \param jer: `int &` Reference to integer error code variable.
+ * \param lcalc: `int &` Reference to integer variable recording the maximum expansion order accounted for.
+ * \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
+);
/*! \brief Bessel function calculation control parameters.
*
@@ -203,17 +117,18 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
* \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;
-}
+double envj(int n, double x);
+
+/*! \brief Compute the Mueller Transformation Matrix.
+ *
+ * This function computes the Mueller Transformation Matrix, or Phase Matrix. See
+ * Sec. 2.8.1 of Borghese, Denti & Saija (2007).
+ *
+ * \param vint: Vector of complex.
+ * \param cmullr: `double **`
+ * \param cmul: `double **`
+ */
+void mmulc(std::complex<double> *vint, double **cmullr, double **cmul);
/*! \brief Starting point for Bessel function magnitude.
*
@@ -224,30 +139,7 @@ double envj(int n, double x) {
* \param mp: `int` Requested order of magnitude.
* \return result: `int` The necessary starting point.
*/
-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;
-}
+int msta1(double x, int mp);
/*! \brief Starting point for Bessel function precision.
*
@@ -259,165 +151,7 @@ int msta1(double x, int mp) {
* \param mp: `int` Requested number of significant digits.
* \return result: `int` The necessary starting point.
*/
-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 Complex Bessel Function.
- *
- * This function computes the complex spherical Bessel funtions \f$j\f$. It uses the
- * auxiliary functions `msta1()` and `msta2()` to determine the starting point for
- * backward recurrence. This is the `CSPHJ` implementation of the `specfun` library.
- *
- * \param n: `int` Order of the function.
- * \param z: `complex<double>` Argumento of the function.
- * \param nm: `int &` Highest computed order.
- * \param csj: Vector of complex. The desired function \f$j\f$.
- */
-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 = abs(csa);
- double abs_csb = abs(csb);
- 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 Compute the Mueller Transformation Matrix.
- *
- * This function computes the Mueller Transformation Matrix, or Phase Matrix. See
- * Sec. 2.8.1 of Borghese, Denti & Saija (2007).
- *
- * \param vint: Vector of complex.
- * \param cmullr: `double **`
- * \param cmul: `double **`
- */
-void mmulc(std::complex<double> *vint, double **cmullr, double **cmul) {
- double sm2 = vint[0].real();
- double s24 = vint[1].real();
- double d24 = vint[1].imag();
- double sm4 = vint[5].real();
- double s23 = vint[8].real();
- double d32 = vint[8].imag();
- double s34 = vint[9].real();
- double d34 = vint[9].imag();
- double sm3 = vint[10].real();
- double s31 = vint[11].real();
- double d31 = vint[11].imag();
- double s21 = vint[12].real();
- double d12 = vint[12].imag();
- double s41 = vint[13].real();
- double d14 = vint[13].imag();
- double sm1 = vint[15].real();
- cmullr[0][0] = sm2;
- cmullr[0][1] = sm3;
- cmullr[0][2] = -s23;
- cmullr[0][3] = -d32;
- cmullr[1][0] = sm4;
- cmullr[1][1] = sm1;
- cmullr[1][2] = -s41;
- cmullr[1][3] = -d14;
- cmullr[2][0] = -s24 * 2.0;
- cmullr[2][1] = -s31 * 2.0;
- cmullr[2][2] = s21 + s34;
- cmullr[2][3] = d34 + d12;
- cmullr[3][0] = -d24 * 2.0;
- cmullr[3][1] = -d31 * 2.0;
- cmullr[3][2] = d34 - d12;
- cmullr[3][3] = s21 - s34;
- cmul[0][0] = (sm2 + sm3 + sm4 + sm1) * 0.5;
- cmul[0][1] = (sm2 - sm3 + sm4 - sm1) * 0.5;
- cmul[0][2] = -s23 - s41;
- cmul[0][3] = -d32 - d14;
- cmul[1][0] = (sm2 + sm3 - sm4 - sm1) * 0.5;
- cmul[1][1] = (sm2 - sm3 - sm4 + sm1) * 0.5;
- cmul[1][2] = -s23 + s41;
- cmul[1][3] = -d32 + d14;
- cmul[2][0] = -s24 - s31;
- cmul[2][1] = -s24 + s31;
- cmul[2][2] = s21 + s34;
- cmul[2][3] = d34 + d12;
- cmul[3][0] = -d24 - d31;
- cmul[3][1] = -d24 + d31;
- cmul[3][2] = d34 - d12;
- cmul[3][3] = s21 - s34;
-}
+int msta2(double x, int n, int mp);
/*! \brief Compute the amplitude of the orthogonal unit vector.
*
@@ -431,23 +165,7 @@ void mmulc(std::complex<double> *vint, double **cmullr, double **cmul) {
* \param iorth: `int`
* \param torth: `double`
*/
-void orunve( double *u1, double *u2, double *u3, int iorth, double torth) {
- if (iorth <= 0) {
- double cp = u1[0] * u2[0] + u1[1] * u2[1] + u1[2] * u2[2];
- double abs_cp = cp;
- if (abs_cp < 0.0) abs_cp *= -1.0;
- if (iorth == 0 || abs_cp >= torth) {
- double fn = 1.0 / sqrt(1.0 - cp * cp);
- u3[0] = (u1[1] * u2[2] - u1[2] * u2[1]) * fn;
- u3[1] = (u1[2] * u2[0] - u1[0] * u2[2]) * fn;
- u3[2] = (u1[0] * u2[1] - u1[1] * u2[0]) * fn;
- return;
- }
- }
- u3[0] = u1[1] * u2[2] - u1[2] * u2[1];
- u3[1] = u1[2] * u2[0] - u1[0] * u2[2];
- u3[2] = u1[0] * u2[1] - u1[1] * u2[0];
-}
+void orunve(double *u1, double *u2, double *u3, int iorth, double torth);
/*! \brief Compute incident and scattered field amplitudes.
*
@@ -465,81 +183,7 @@ void orunve( double *u1, double *u2, double *u3, int iorth, double torth) {
void pwma(
double *up, double *un, std::complex<double> *ylm, int inpol, int lw,
int isq, C1 *c1
- ) {
- const double four_pi = 8.0 * acos(0.0);
- int is = isq;
- if (isq == -1) is = 0;
- int ispo = is + 1;
- int ispt = is + 2;
- int nlwm = lw * (lw + 2);
- int nlwmt = nlwm + nlwm;
- const double sqrtwi = 1.0 / sqrt(2.0);
- const std::complex<double> uim(0.0, 1.0);
- std::complex<double> cm1 = 0.5 * std::complex<double>(up[0], up[1]);
- std::complex<double> cp1 = 0.5 * std::complex<double>(up[0], -up[1]);
- double cz1 = up[2];
- std::complex<double> cm2 = 0.5 * std::complex<double>(un[0], un[1]);
- std::complex<double> cp2 = 0.5 * std::complex<double>(un[0], -un[1]);
- double cz2 = un[2];
- for (int l20 = 0; l20 < lw; l20++) {
- int l = l20 + 1;
- int lf = l + 1;
- int lftl = lf * l;
- double x = 1.0 * lftl;
- std::complex<double> cl = std::complex<double>(four_pi / sqrt(x), 0.0);
- for (int powi = 1; powi <= l; powi++) cl *= uim;
- int mv = l + lf;
- int m = -lf;
- for (int mf20 = 0; mf20 < mv; mf20++) {
- m += 1;
- int k = lftl + m;
- x = 1.0 * (lftl - m * (m + 1));
- double cp = sqrt(x);
- x = 1.0 * (lftl - m * (m - 1));
- double cm = sqrt(x);
- double cz = 1.0 * m;
- c1->w[k - 1][ispo - 1] = dconjg(
- cp1 * cp * ylm[k + 1] +
- cm1 * cm * ylm[k - 1] +
- cz1 * cz * ylm[k]
- ) * cl;
- c1->w[k - 1][ispt - 1] = dconjg(
- cp2 * cp * ylm[k + 1] +
- cm2 * cm * ylm[k - 1] +
- cz2 * cz * ylm[k]
- ) * cl;
- }
- }
- for (int k30 = 0; k30 < nlwm; k30++) {
- int i = k30 + nlwm;
- c1->w[i][ispo - 1] = uim * c1->w[k30][ispt - 1];
- c1->w[i][ispt - 1] = -uim * c1->w[k30][ispo - 1];
- }
- if (inpol != 0) {
- for (int k40 = 0; k40 < nlwm; k40++) {
- int i = k40 + nlwm;
- std::complex<double> cc1 = sqrtwi * (c1->w[k40][ispo - 1] + uim * c1->w[k40][ispt - 1]);
- std::complex<double> cc2 = sqrtwi * (c1->w[k40][ispo - 1] - uim * c1->w[k40][ispt - 1]);
- c1->w[k40][ispo - 1] = cc2;
- c1->w[i][ispo - 1] = -cc2;
- c1->w[k40][ispt - 1] = cc1;
- c1->w[i][ispt - 1] = cc1;
- }
- } else {
- if (isq == 0) {
- return;
- }
- }
- if (isq != 0) {
- for (int i50 = 0; i50 < 2; i50++) {
- int ipt = i50 + 2;
- int ipis = i50 + is;
- for (int k50 = 0; k50 < nlwmt; k50++) {
- c1->w[k50][ipt] = dconjg(c1->w[k50][ipis]);
- }
- }
- }
-}
+);
/*! \brief Compute radiation torques on particles.
*
@@ -559,58 +203,7 @@ void pwma(
void rabas(
int inpol, int li, int nsph, C1 *c1, double **tqse, std::complex<double> **tqspe,
double **tqss, std::complex<double> **tqsps
- ) {
- std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
- std::complex<double> uim = std::complex<double>(0.0, 1.0);
- double two_pi = 4.0 * acos(0.0);
- for (int i80 = 0; i80 < nsph; i80++) {
- int i = i80 + 1;
- if(c1->iog[i80] >= i) {
- tqse[0][i80] = 0.0;
- tqse[1][i80] = 0.0;
- tqspe[0][i80] = cc0;
- tqspe[1][i80] = cc0;
- tqss[0][i80] = 0.0;
- tqss[1][i80] = 0.0;
- tqsps[0][i80] = cc0;
- tqsps[1][i80] = cc0;
- for (int l70 = 0; l70 < li; l70++) {
- int l = l70 + 1;
- double fl = 1.0 * (l + l + 1);
- std::complex<double> rm = 1.0 / c1->rmi[l70][i80];
- double rmm = (rm * dconjg(rm)).real();
- std::complex<double> re = 1.0 / c1->rei[l70][i80];
- double rem = (re * dconjg(re)).real();
- if (inpol == 0) {
- std::complex<double> pce = ((rm + re) * uim) * fl;
- std::complex<double> pcs = ((rmm + rem) * fl) * uim;
- tqspe[0][i80] -= pce;
- tqspe[1][i80] += pce;
- tqsps[0][i80] -= pcs;
- tqsps[1][i80] += pcs;
- } else {
- double ce = (rm + re).real() * fl;
- double cs = (rmm + rem) * fl;
- tqse[0][i80] -= ce;
- tqse[1][i80] += ce;
- tqss[0][i80] -= cs;
- tqss[1][i80] += cs;
- }
- }
- if (inpol == 0) {
- tqspe[0][i80] *= two_pi;
- tqspe[1][i80] *= two_pi;
- tqsps[0][i80] *= two_pi;
- tqsps[1][i80] *= two_pi;
- } else {
- tqse[0][i80] *= two_pi;
- tqse[1][i80] *= two_pi;
- tqss[0][i80] *= two_pi;
- tqss[1][i80] *= two_pi;
- }
- }
- }
-}
+);
/*! \brief Real Bessel Function.
*
@@ -623,61 +216,7 @@ void rabas(
* \param nm: `int &` Highest computed order.
* \param sj: `double[]` The desired function \f$j\f$.
*/
-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 = 1; k <= n; k++)
- sj[k] = 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 < 0.0) ? -sa : sa;
- double abs_sb = (sb < 0.0) ? -sb : sb;
- if (abs_sa > abs_sb) cs = sa / f;
- else cs = sb / f0;
- for (int k = 0; k <= nm; k++) {
- sj[k] = cs * sj[k];
- }
-}
+void rbf(int n, double x, int &nm, double sj[]);
/*! \brief Soft layer radial function and derivative.
*
@@ -698,39 +237,7 @@ 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 = 0; 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 Transition layer radial function and derivative.
*
@@ -751,50 +258,7 @@ 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 = 0; ipnt60 < npntmo; ipnt60++) {
- int ipnt = ipnt60 + 1;
- int jpnt = ipnt + ipnt - 1;
- int jpnt60 = jpnt - 1;
- cy1 = cl / (x * x) - c2->ris[jpnt60];
- 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[jpnt];
- 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[jpntpo];
- 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[jpnt60];
- cdy1 += 2.0 * c2->dlri[jpnt60];
- c21 = (cy1 * y2 + cdy1 * dy2) * step;
- cy23 -= cdy23 * c2->dlri[jpnt];
- cdy23 += 2.0 * c2->dlri[jpnt];
- 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[jpntpo];
- cdy4 += 2.0 * c2->dlri[jpntpo];
- 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 Spherical Bessel functions.
*
@@ -806,49 +270,24 @@ void rkt(
* \param nm: `int &` Highest computed order.
* \param sy: `double[]` The desired function \f$y\f$.
*/
-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;
-}
+void rnf(int n, double x, int &nm, double sy[]);
+
+/*! \brief Spherical harmonics for given direction.
+ *
+ * This function computes the field spherical harmonics for a given direction. See Sec.
+ * 1.5.2 in Borghese, Denti & Saija (2007).
+ *
+ * \param cosrth: `double` Cosine of direction's elevation.
+ * \param sinrth: `double` Sine of direction's elevation.
+ * \param cosrph: `double` Cosine of direction's azimuth.
+ * \param sinrph: `double` Sine of direction's azimuth.
+ * \param ll: `int` L value expansion order.
+ * \param ylm: Vector of complex. The requested spherical harmonics.
+ */
+void sphar(
+ double cosrth, double sinrth, double cosrph, double sinrph,
+ int ll, std::complex<double> *ylm
+);
/*! \brief Compute scattering, absorption and extinction cross-sections.
*
@@ -862,46 +301,7 @@ void rnf(int n, double x, int &nm, double sy[]) {
* \param exri: `double` External medium refractive index.
* \param c1: `C1 *` Pointer to a `C1` data structure.
*/
-void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
- std::complex<double> sum21, rm, re, csam;
- const std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
- const 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 = 0; i12 < nsph; i12++) {
- int i = i12 + 1;
- int iogi = c1->iog[i12];
- if (iogi >= i) {
- double sums = 0.0;
- std::complex<double> sum21 = cc0;
- for (int l10 = 0; l10 < lm; l10++) {
- int l = l10 + 1;
- double fl = 1.0 + l + l;
- rm = 1.0 / c1->rmi[l10][i12];
- re = 1.0 / c1->rei[l10][i12];
- std::complex<double> rm_cnjg = dconjg(rm);
- std::complex<double> re_cnjg = dconjg(re);
- 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] = scasec;
- c1->sexs[i12] = extsec;
- c1->sabs[i12] = abssec;
- double gcss = c1->gcsv[i12];
- c1->sqscs[i12] = scasec / gcss;
- c1->sqexs[i12] = extsec / gcss;
- c1->sqabs[i12] = abssec / gcss;
- c1->fsas[i12] = sum21 * csam;
- }
- tfsas += c1->fsas[iogi - 1];
- }
-}
+void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1);
/*! \brief C++ Compute the scattering amplitude and the scattered field intensity.
*
@@ -914,152 +314,7 @@ void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri
* \param exri: `double` External medium refractive index.
* \param c1: `C1 *` Pointer to a `C1` data structure.
*/
-void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
- std::complex<double> s11, s21, s12, s22, rm, re, csam, cc;
- const std::complex<double> cc0(0.0, 0.0);
- double ccs = 1.0 / (vk * vk);
- csam = -(ccs / (exri * vk)) * std::complex<double>(0.0, 0.5);
- const double pigfsq = 64.0 * acos(0.0) * acos(0.0);
- double cfsq = 4.0 / (pigfsq * ccs * ccs);
- int nlmm = lm * (lm + 2);
- for (int i14 = 0; i14 < nsph; i14++) {
- int i = i14 + 1;
- int iogi = c1->iog[i14];
- if (iogi >= i) {
- int k = 0;
- s11 = cc0;
- s21 = cc0;
- s12 = cc0;
- s22 = cc0;
- for (int l10 = 0; l10 < lm; l10++) {
- int l = l10 + 1;
- rm = 1.0 / c1->rmi[l10][i14];
- re = 1.0 / c1->rei[l10][i14];
- int ltpo = l + l + 1;
- for (int im10 = 0; im10 < ltpo; im10++) {
- k += 1;
- int ke = k + nlmm;
- s11 = s11 - c1->w[k - 1][2] * c1->w[k - 1][0] * rm - c1->w[ke - 1][2] * c1->w[ke - 1][0] * re;
- s21 = s21 - c1->w[k - 1][3] * c1->w[k - 1][0] * rm - c1->w[ke - 1][3] * c1->w[ke - 1][0] * re;
- s12 = s12 - c1->w[k - 1][2] * c1->w[k - 1][1] * rm - c1->w[ke - 1][2] * c1->w[ke - 1][1] * re;
- s22 = s22 - c1->w[k - 1][3] * c1->w[k - 1][1] * rm - c1->w[ke - 1][3] * c1->w[ke - 1][1] * re;
- }
- }
- c1->sas[i14][0][0] = s11 * csam;
- c1->sas[i14][1][0] = s21 * csam;
- c1->sas[i14][0][1] = s12 * csam;
- c1->sas[i14][1][1] = s22 * csam;
- }
- } // loop i14
- for (int i24 = 0; i24 < nsph; i24++) {
- int i = i24 + 1;
- int iogi = c1->iog[i24];
- if (iogi >= i) {
- int j = 0;
- for (int ipo1 = 0; ipo1 < 2; ipo1++) {
- for (int jpo1 = 0; jpo1 < 2; jpo1++) {
- std::complex<double> cc = dconjg(c1->sas[i24][jpo1][ipo1]);
- for (int ipo2 = 0; ipo2 < 2; ipo2++) {
- for (int jpo2 = 0; jpo2 < 2; jpo2++) {
- c1->vints[i24][j++] = c1->sas[i24][jpo2][ipo2] * cc * cfsq;
- }
- }
- }
- }
- }
- }
-}
-
-/*! \brief Spherical harmonics for given direction.
- *
- * This function computes the field spherical harmonics for a given direction. See Sec.
- * 1.5.2 in Borghese, Denti & Saija (2007).
- *
- * \param cosrth: `double` Cosine of direction's elevation.
- * \param sinrth: `double` Sine of direction's elevation.
- * \param cosrph: `double` Cosine of direction's azimuth.
- * \param sinrph: `double` Sine of direction's azimuth.
- * \param ll: `int` L value expansion order.
- * \param ylm: Vector of complex. The requested spherical harmonics.
- */
-void sphar(
- double cosrth, double sinrth, double cosrph, double sinrph,
- int ll, std::complex<double> *ylm
- ) {
- const int rmp_size = ll;
- const int plegn_size = (ll + 1) * ll / 2 + ll + 1;
- double sinrmp[rmp_size], cosrmp[rmp_size], plegn[plegn_size];
- double four_pi = 8.0 * acos(0.0);
- double pi4irs = 1.0 / sqrt(four_pi);
- double x = cosrth;
- double y = sinrth;
- if (y < 0.0) y *= -1.0;
- double cllmo = 3.0;
- double cll = 1.5;
- double ytol = y;
- plegn[0] = 1.0;
- plegn[1] = x * sqrt(cllmo);
- plegn[2] = ytol * sqrt(cll);
- sinrmp[0] = sinrph;
- cosrmp[0] = cosrph;
- if (ll >= 2) {
- int k = 3;
- for (int l20 = 2; l20 <= ll; l20++) {
- int lmo = l20 - 1;
- int ltpo = l20 + l20 + 1;
- int ltmo = ltpo - 2;
- int lts = ltpo * ltmo;
- double cn = 1.0 * lts;
- for (int mpo10 = 1; mpo10 <= lmo; mpo10++) {
- int m = mpo10 - 1;
- int mpopk = mpo10 + k;
- int ls = (l20 + m) * (l20 - m);
- double cd = 1.0 * ls;
- double cnm = 1.0 * ltpo * (lmo + m) * (l20 - mpo10);
- double cdm = 1.0 * ls * (ltmo - 2);
- plegn[mpopk - 1] = plegn[mpopk - l20 - 1] * x * sqrt(cn / cd) -
- plegn[mpopk - ltmo - 1] * sqrt(cnm / cdm);
- }
- int lpk = l20 + k;
- double cltpo = 1.0 * ltpo;
- plegn[lpk - 1] = plegn[k - 1] * x * sqrt(cltpo);
- k = lpk + 1;
- double clt = 1.0 * (ltpo - 1);
- cll *= (cltpo / clt);
- ytol *= y;
- plegn[k - 1] = ytol * sqrt(cll);
- sinrmp[l20 - 1] = sinrph * cosrmp[lmo - 1] + cosrph * sinrmp[lmo - 1];
- cosrmp[l20 - 1] = cosrph * cosrmp[lmo - 1] - sinrph * sinrmp[lmo - 1];
- } // end l20 loop
- }
- // label 30
- int l = 0;
- int m, k, l0y, l0p, lmy, lmp;
- double save;
- label40:
- m = 0;
- k = l * (l + 1);
- l0y = k + 1;
- l0p = k / 2 + 1;
- ylm[l0y - 1] = pi4irs * plegn[l0p - 1];
- goto label45;
- label44:
- lmp = l0p + m;
- save = pi4irs * plegn[lmp - 1];
- lmy = l0y + m;
- ylm[lmy - 1] = save * std::complex<double>(cosrmp[m - 1], sinrmp[m - 1]);
- if (m % 2 != 0) ylm[lmy - 1] *= -1.0;
- lmy = l0y - m;
- ylm[lmy - 1] = save * std::complex<double>(cosrmp[m - 1], -sinrmp[m - 1]);
- label45:
- if (m >= l) goto label47;
- m += 1;
- goto label44;
- label47:
- if (l >= ll) return;
- l += 1;
- goto label40;
-}
+void sscr2(int nsph, int lm, double vk, double exri, C1 *c1);
/*! \brief Determine the geometrical asymmetry parameter coefficients.
*
@@ -1070,42 +325,7 @@ void sphar(
* \param lm: `int` Maximum field expansion order.
* \param zpv: `double ****` Matrix of geometrical asymmetry parameter coefficients.
*/
-void thdps(int lm, double ****zpv) {
- //for (int l10 = 0; l10 < lm; l10++) { // 0-init, can be omitted
- // for (int ilmp = 0; ilmp < 3; ilmp++) {
- // zpv[l10][ilmp][0][0] = 0.0;
- // zpv[l10][ilmp][0][1] = 0.0;
- // zpv[l10][ilmp][1][0] = 0.0;
- // zpv[l10][ilmp][1][1] = 0.0;
- // }
- //}
- for (int l15 = 0; l15 < lm; l15++) {
- int l = l15 + 1;
- double xd = 1.0 * l * (l + 1);
- double zp = -1.0 / sqrt(xd);
- zpv[l15][1][0][1] = zp;
- zpv[l15][1][1][0] = zp;
- }
- if (lm != 1) {
- for (int l20 = 1; l20 < lm; l20++) {
- int l = l20 + 1;
- double xn = 1.0 * (l - 1) * (l + 1);
- double xd = 1.0 * l * (l + l + 1);
- double zp = sqrt(xn / xd);
- zpv[l20][0][0][0] = zp;
- zpv[l20][0][1][1] = zp;
- }
- int lmmo = lm - 1;
- for (int l25 = 0; l25 < lmmo; l25++) {
- int l = l25 + 1;
- double xn = 1.0 * l * (l + 2);
- double xd = (l + 1) * (l + l + 1);
- double zp = -1.0 * sqrt(xn / xd);
- zpv[l25][2][0][0] = zp;
- zpv[l25][2][1][1] = zp;
- }
- }
-}
+void thdps(int lm, double ****zpv);
/*! \brief C++ porting of UPVMP
*
@@ -1123,32 +343,7 @@ void thdps(int lm, double ****zpv) {
void upvmp(
double thd, double phd, int icspnv, double &cost, double &sint,
double &cosp, double &sinp, double *u, double *up, double *un
- ) {
- double half_pi = acos(0.0);
- double rdr = half_pi / 90.0;
- double th = thd * rdr;
- double ph = phd * rdr;
- cost = cos(th);
- sint = sin(th);
- cosp = cos(ph);
- sinp = sin(ph);
- u[0] = cosp * sint;
- u[1] = sinp * sint;
- u[2] = cost;
- up[0] = cosp * cost;
- up[1] = sinp * cost;
- up[2] = -sint;
- un[0] = -sinp;
- un[1] = cosp;
- un[2] = 0.0;
- if (icspnv != 0) {
- up[0] *= -1.0;
- up[1] *= -1.0;
- up[2] *= -1.0;
- un[0] *= -1.0;
- un[1] *= -1.0;
- }
-}
+);
/*! \brief Compute the unitary vector perpendicular to incident and scattering plane.
*
@@ -1179,69 +374,7 @@ void upvsp(
double *u, double *upmp, double *unmp, double *us, double *upsmp, double *unsmp,
double *up, double *un, double *ups, double *uns, double *duk, int &isq,
int &ibf, double &scand, double &cfmp, double &sfmp, double &cfsp, double &sfsp
- ) {
- double rdr = acos(0.0) / 90.0;
- double small = 1.0e-6;
- isq = 0;
- scand = u[0] * us[0] + u[1] * us[1] + u[2] * us[2];
- double abs_scand = (scand >= 1.0) ? scand - 1.0 : 1.0 - scand;
- if (abs_scand >= small) {
- abs_scand = scand + 1.0;
- if (abs_scand < 0.0) abs_scand *= -1.0;
- if (abs_scand >= small) {
- scand = acos(scand) / rdr;
- duk[0] = u[0] - us[0];
- duk[1] = u[1] - us[1];
- duk[2] = u[2] - us[2];
- ibf = 0;
- } else { // label 15
- scand = 180.0;
- duk[0] = 2.0 * u[0];
- duk[1] = 2.0 * u[1];
- duk[2] = 2.0 * u[2];
- ibf = 1;
- ups[0] = -upsmp[0];
- ups[1] = -upsmp[1];
- ups[2] = -upsmp[2];
- uns[0] = -unsmp[0];
- uns[1] = -unsmp[1];
- uns[2] = -unsmp[2];
- }
- } else { // label 10
- scand = 0.0;
- duk[0] = 0.0;
- duk[1] = 0.0;
- duk[2] = 0.0;
- ibf = -1;
- isq = -1;
- ups[0] = upsmp[0];
- ups[1] = upsmp[1];
- ups[2] = upsmp[2];
- uns[0] = unsmp[0];
- uns[1] = unsmp[1];
- uns[2] = unsmp[2];
- }
- if (ibf == -1 || ibf == 1) { // label 20
- up[0] = upmp[0];
- up[1] = upmp[1];
- up[2] = upmp[2];
- un[0] = unmp[0];
- un[1] = unmp[1];
- un[2] = unmp[2];
- } else { // label 25
- orunve(u, us, un, -1, small);
- uns[0] = un[0];
- uns[1] = un[1];
- uns[2] = un[2];
- orunve(un, u, up, 1, small);
- orunve(uns, us, ups, 1, small);
- }
- // label 85
- cfmp = upmp[0] * up[0] + upmp[1] * up[1] + upmp[2] * up[2];
- sfmp = unmp[0] * up[0] + unmp[1] * up[1] + unmp[2] * up[2];
- cfsp = ups[0] * upsmp[0] + ups[1] * upsmp[1] + ups[2] * upsmp[2];
- sfsp = uns[0] * upsmp[0] + uns[1] * upsmp[1] + uns[2] * upsmp[2];
-}
+);
/*! \brief Compute meridional plane-referred geometrical asymmetry parameter coefficients.
*
@@ -1268,29 +401,7 @@ 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
- ) {
- const int ylm_size = (lm + 1) * (lm + 1) + 1;
- std::complex<double> *ylm = new std::complex<double>[ylm_size];
- const int nlmp = lm * (lm + 2) + 2;
- ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
- if (idot != 0) {
- if (idot != 1) {
- for (int n40 = 0; n40 < nsph; n40++) {
- arg[n40] = u[0] * c1->rxx[n40] + u[1] * c1->ryy[n40] + u[2] * c1->rzz[n40];
- }
- } else {
- for (int n50 = 0; n50 < nsph; n50++) {
- arg[n50] = c1->rzz[n50];
- }
- }
- if (iis == 2) {
- for (int n60 = 0; n60 < nsph; n60++) arg[n60] *= -1;
- }
- }
- sphar(cost, sint, cosp, sinp, lm, ylm);
- pwma(up, un, ylm, inpol, lm, iis, c1);
- delete[] ylm;
-}
+);
/*! \brief Compute the scattering plane-referred geometrical asymmetry parameter coefficients.
*
@@ -1327,200 +438,6 @@ void wmasp(
double cosps, double sinps, double *u, double *up, double *un, double *us,
double *ups, double *uns, int isq, int ibf, int inpol, int lm, int idot,
int nsph, double *argi, double *args, C1 *c1
- ) {
- const int ylm_size = (lm + 1) * (lm + 1) + 1;
- std::complex<double> *ylm = new std::complex<double>[ylm_size];
- const int nlmp = lm * (lm + 2) + 2;
- ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
- if (idot != 0) {
- if (idot != 1) {
- for (int n40 = 0; n40 < nsph; n40++) {
- argi[n40] = u[0] * c1->rxx[n40] + u[1] * c1->ryy[n40] + u[2] * c1->rzz[n40];
- if (ibf != 0) {
- args[n40] = argi[n40] * ibf;
- } else {
- args[n40] = -1.0 * (us[0] * c1->rxx[n40] + us[1] * c1->ryy[n40] + us[2] * c1->rzz[n40]);
- }
- }
- } else { // label 50
- for (int n60 = 0; n60 < nsph; n60++) {
- argi[n60] = cost * c1->rzz[n60];
- if (ibf != 0) {
- args[n60] = argi[n60] * ibf;
- } else {
- args[n60] = -costs * c1->rzz[n60];
- }
- }
- }
- }
- sphar(cost, sint, cosp, sinp, lm, ylm);
- pwma(up, un, ylm, inpol, lm, isq, c1);
- if (ibf >= 0) {
- sphar(costs, sints, cosps, sinps, lm, ylm);
- pwma(ups, uns, ylm, inpol, lm, 2, c1);
- }
- delete[] ylm;
-}
-
-/*! \brief Compute Mie scattering coefficients.
- *
- * This function determines the L-dependent Mie scattering coefficients \f$a_l\f$ and \f$b_l\f$
- * for the cases of homogeneous spheres, radially non-homogeneous spheres and, in case of sphere
- * with dielectric function, sustaining longitudinal waves. See Sec. 5.1 in Borghese, Denti
- * & Saija (2007).
- *
- * \param li: `int` Maximum field expansion order.
- * \param i: `int`
- * \param npnt: `int`
- * \param npntts: `int`
- * \param vk: `double` Wave number in scale units.
- * \param exdc: `double` External medium dielectric constant.
- * \param exri: `double` External medium refractive index.
- * \param c1: `C1 *` Pointer to a `C1` data structure.
- * \param c2: `C2 *` Pointer to a `C2` data structure.
- * \param jer: `int &` Reference to integer error code variable.
- * \param lcalc: `int &` Reference to integer variable recording the maximum expansion order accounted for.
- * \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
-) {
- const int lipo = li + 1;
- const int lipt = li + 2;
- double *rfj = new double[lipt];
- double *rfn = new double[lipt];
- std::complex<double> cfj[lipt], fbi[lipt], fb[lipt], fn[lipt];
- std::complex<double> rmf[li], drmf[li], ref[li], dref[li];
- 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;
- 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;
- delete[] rfj;
- delete[] rfn;
- 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;
- delete[] rfj;
- delete[] rfn;
- 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;
- delete[] rfj;
- delete[] rfn;
- return;
- }
- rnf(lipo, arex, lcalc, rfn);
- if (lcalc < lipo) {
- jer = 8;
- delete[] rfj;
- delete[] rfn;
- return;
- }
- for (int j43 = 1; j43 <= lipt; j43++) {
- fb[j43 - 1] = rfj[j43 - 1];
- fn[j43 - 1] = rfn[j43 - 1];
- }
- 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 ?
- delete[] rfj;
- delete[] rfn;
- return;
-}
+);
#endif /* SRC_INCLUDE_SPH_SUBS_H_ */
diff --git a/src/libnptm/sph_subs.cpp b/src/libnptm/sph_subs.cpp
new file mode 100644
index 00000000..02c2d5a9
--- /dev/null
+++ b/src/libnptm/sph_subs.cpp
@@ -0,0 +1,1167 @@
+/*! \file sph_subs.cpp
+ *
+ * \brief C++ implementation of SPHERE subroutines.
+ */
+
+#include "../include/sph_subs.h"
+
+using namespace std;
+
+void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
+ complex<double> cc0 = complex<double>(0.0, 0.0);
+ complex<double> summ, sume, suem, suee, sum;
+ double half_pi = acos(0.0);
+ double cofs = half_pi * 2.0 / sqk;
+ for (int i40 = 0; i40 < nsph; i40++) {
+ int i = i40 + 1;
+ int iogi = c1->iog[i40];
+ if (iogi >= i) {
+ sum = cc0;
+ for (int l30 = 0; l30 < li; l30++) {
+ int l = l30 + 1;
+ int ltpo = l + l + 1;
+ for (int ilmp = 1; ilmp < 4; ilmp++) {
+ int ilmp30 = ilmp - 1;
+ bool goto30 = (l == 1 && ilmp == 1) || (l == li && ilmp == 3);
+ if (!goto30) {
+ int lmpml = ilmp - 2;
+ int lmp = l + lmpml;
+ double cofl = sqrt(1.0 * (ltpo * (lmp + lmp + 1)));
+ summ = zpv[l30][ilmp30][0][0] /
+ (
+ dconjg(c1->rmi[l30][i40]) *
+ c1->rmi[lmp - 1][i40]
+ );
+ sume = zpv[l30][ilmp30][0][1] /
+ (
+ dconjg(c1->rmi[l30][i40]) *
+ c1->rei[lmp - 1][i40]
+ );
+ suem = zpv[l30][ilmp30][1][0] /
+ (
+ dconjg(c1->rei[l30][i40]) *
+ c1->rmi[lmp - 1][i40]
+ );
+ suee = zpv[l30][ilmp30][1][1] /
+ (
+ dconjg(c1->rei[l30][i40]) *
+ c1->rei[lmp - 1][i40]
+ );
+ sum += (cg1(lmpml, 0, l, -1) * (summ - sume - suem + suee) +
+ cg1(lmpml, 0, l, 1) * (summ + sume + suem + suee)) * cofl;
+ }
+ }
+ }
+ }
+ gaps[i40] = sum.real() * cofs;
+ }
+}
+
+void cbf(int n, complex<double> z, int &nm, 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] = 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;
+ }
+ complex<double> csa = csj[0];
+ complex<double> csb = csj[1];
+ int m = msta1(a0, 200);
+ if (m < n) nm = m;
+ else m = msta2(a0, n, 15);
+ complex<double> cf0 = 0.0;
+ complex<double> cf1 = 1.0e-100;
+ 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 = abs(csa);
+ double abs_csb = abs(csb);
+ if (abs_csa > abs_csb) cs = csa / cf;
+ else cs = csb / cf0;
+ for (int k = 0; k <= nm; k++) {
+ csj[k] = cs * csj[k];
+ }
+}
+
+double cg1(int lmpml, int mu, int l, int m) {
+ double result = 0.0;
+ double xd, xn;
+ if (lmpml == -1) { // Interpreted as GOTO 30
+ xd = 2.0 * l * (2 * l - 1);
+ if (mu == -1) {
+ xn = 1.0 * (l - 1 - m) * (l - m);
+ } else if (mu == 0) {
+ xn = 2.0 * (l - m) * (l + m);
+ } else if (mu == 1) {
+ xn = 1.0 * (l - 1 + m) * (l + m);
+ } else {
+ throw 111; // Need an exception for unpredicted indices.
+ }
+ result = sqrt(xn / xd);
+ } else if (lmpml == 0) { // Interpreted as GOTO 5
+ bool goto10 = (m != 0) || (mu != 0);
+ if (!goto10) {
+ result = 0.0;
+ return result;
+ }
+ if (mu != 0) {
+ xd = 2.0 * l * (l + 1);
+ if (mu == -1) {
+ xn = 1.0 * (l - m) * (l + m + 1);
+ result = -sqrt(xn / xd);
+ } else if (mu == 1) { // mu > 0
+ xn = 1.0 * (l + m) * (l - m + 1);
+ result = sqrt(xn / xd);
+ } else {
+ throw 111; // Need an exception for unpredicted indices.
+ }
+ } else { // mu = 0
+ xd = 1.0 * l * (l + 1);
+ xn = -1.0 * m;
+ result = xn / sqrt(xd);
+ }
+ } else if (lmpml == 1) { // Interpreted as GOTO 60
+ xd = 2.0 * (l * 2 + 3) * (l + 1);
+ if (mu == -1) {
+ xn = 1.0 * (l + 1 + m) * (l + 2 + m);
+ result = sqrt(xn / xd);
+ } else if (mu == 0) {
+ xn = 2.0 * (l + 1 - m) * (l + 1 + m);
+ result = -sqrt(xn / xd);
+ } else if (mu == 1) {
+ xn = 1.0 * (l + 1 - m) * (l + 2 - m);
+ result = sqrt(xn / xd);
+ } else { // mu was not recognized.
+ throw 111; // Need an exception for unpredicted indices.
+ }
+ } else { // lmpml was not recognized
+ throw 111; // Need an exception for unpredicted indices.
+ }
+ return result;
+}
+
+complex<double> dconjg(complex<double> z) {
+ double zreal = z.real();
+ double zimag = z.imag();
+ return complex<double>(zreal, -zimag);
+}
+
+void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
+ const double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1];
+ const double half_step = 0.5 * dif / npntmo;
+ double rr = c1->rc[i - 1][ns - 1];
+ const complex<double> delta = c2->dc0[ic] - c2->dc0[ic - 1];
+ const int kpnt = npntmo + npntmo;
+ c2->ris[kpnt] = c2->dc0[ic];
+ c2->dlri[kpnt] = complex<double>(0.0, 0.0);
+ const int i90 = i - 1;
+ const int ns90 = ns - 1;
+ const int ic90 = ic - 1;
+ for (int np90 = 0; np90 < kpnt; np90++) {
+ double ff = (rr - c1->rc[i90][ns90]) / dif;
+ c2->ris[np90] = delta * ff * ff * (-2.0 * ff + 3.0) + c2->dc0[ic90];
+ c2->dlri[np90] = 3.0 * delta * ff * (1.0 - ff) / (dif * vk * c2->ris[np90]);
+ rr += half_step;
+ }
+}
+
+void dme(
+ int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
+ C1 *c1, C2 *c2, int &jer, int &lcalc, complex<double> &arg
+) {
+ const int lipo = li + 1;
+ const int lipt = li + 2;
+ double *rfj = new double[lipt];
+ double *rfn = new double[lipt];
+ complex<double> cfj[lipt], fbi[lipt], fb[lipt], fn[lipt];
+ complex<double> rmf[li], drmf[li], ref[li], dref[li];
+ complex<double> dfbi, dfb, dfn, ccna, ccnb, ccnc, ccnd;
+ complex<double> y1, dy1, y2, dy2, arin, cri, uim;
+ jer = 0;
+ uim = complex<double>(0.0, 1.0);
+ int nstp = npnt - 1;
+ int nstpts = npntts - 1;
+ 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;
+ delete[] rfj;
+ delete[] rfn;
+ 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;
+ delete[] rfj;
+ delete[] rfn;
+ 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;
+ delete[] rfj;
+ delete[] rfn;
+ return;
+ }
+ rnf(lipo, arex, lcalc, rfn);
+ if (lcalc < lipo) {
+ jer = 8;
+ delete[] rfj;
+ delete[] rfn;
+ return;
+ }
+ for (int j43 = 1; j43 <= lipt; j43++) {
+ fb[j43 - 1] = rfj[j43 - 1];
+ fn[j43 - 1] = rfn[j43 - 1];
+ }
+ 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);
+ 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);
+ }
+ } // nsh <= 1 ?
+ delete[] rfj;
+ delete[] rfn;
+ return;
+}
+
+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;
+}
+
+void mmulc(complex<double> *vint, double **cmullr, double **cmul) {
+ double sm2 = vint[0].real();
+ double s24 = vint[1].real();
+ double d24 = vint[1].imag();
+ double sm4 = vint[5].real();
+ double s23 = vint[8].real();
+ double d32 = vint[8].imag();
+ double s34 = vint[9].real();
+ double d34 = vint[9].imag();
+ double sm3 = vint[10].real();
+ double s31 = vint[11].real();
+ double d31 = vint[11].imag();
+ double s21 = vint[12].real();
+ double d12 = vint[12].imag();
+ double s41 = vint[13].real();
+ double d14 = vint[13].imag();
+ double sm1 = vint[15].real();
+ cmullr[0][0] = sm2;
+ cmullr[0][1] = sm3;
+ cmullr[0][2] = -s23;
+ cmullr[0][3] = -d32;
+ cmullr[1][0] = sm4;
+ cmullr[1][1] = sm1;
+ cmullr[1][2] = -s41;
+ cmullr[1][3] = -d14;
+ cmullr[2][0] = -s24 * 2.0;
+ cmullr[2][1] = -s31 * 2.0;
+ cmullr[2][2] = s21 + s34;
+ cmullr[2][3] = d34 + d12;
+ cmullr[3][0] = -d24 * 2.0;
+ cmullr[3][1] = -d31 * 2.0;
+ cmullr[3][2] = d34 - d12;
+ cmullr[3][3] = s21 - s34;
+ cmul[0][0] = (sm2 + sm3 + sm4 + sm1) * 0.5;
+ cmul[0][1] = (sm2 - sm3 + sm4 - sm1) * 0.5;
+ cmul[0][2] = -s23 - s41;
+ cmul[0][3] = -d32 - d14;
+ cmul[1][0] = (sm2 + sm3 - sm4 - sm1) * 0.5;
+ cmul[1][1] = (sm2 - sm3 - sm4 + sm1) * 0.5;
+ cmul[1][2] = -s23 + s41;
+ cmul[1][3] = -d32 + d14;
+ cmul[2][0] = -s24 - s31;
+ cmul[2][1] = -s24 + s31;
+ cmul[2][2] = s21 + s34;
+ cmul[2][3] = d34 + d12;
+ cmul[3][0] = -d24 - d31;
+ cmul[3][1] = -d24 + d31;
+ cmul[3][2] = d34 - d12;
+ cmul[3][3] = s21 - s34;
+}
+
+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;
+}
+
+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;
+}
+
+void orunve( double *u1, double *u2, double *u3, int iorth, double torth) {
+ if (iorth <= 0) {
+ double cp = u1[0] * u2[0] + u1[1] * u2[1] + u1[2] * u2[2];
+ double abs_cp = cp;
+ if (abs_cp < 0.0) abs_cp *= -1.0;
+ if (iorth == 0 || abs_cp >= torth) {
+ double fn = 1.0 / sqrt(1.0 - cp * cp);
+ u3[0] = (u1[1] * u2[2] - u1[2] * u2[1]) * fn;
+ u3[1] = (u1[2] * u2[0] - u1[0] * u2[2]) * fn;
+ u3[2] = (u1[0] * u2[1] - u1[1] * u2[0]) * fn;
+ return;
+ }
+ }
+ u3[0] = u1[1] * u2[2] - u1[2] * u2[1];
+ u3[1] = u1[2] * u2[0] - u1[0] * u2[2];
+ u3[2] = u1[0] * u2[1] - u1[1] * u2[0];
+}
+
+void pwma(
+ double *up, double *un, complex<double> *ylm, int inpol, int lw,
+ int isq, C1 *c1
+) {
+ const double four_pi = 8.0 * acos(0.0);
+ int is = isq;
+ if (isq == -1) is = 0;
+ int ispo = is + 1;
+ int ispt = is + 2;
+ int nlwm = lw * (lw + 2);
+ int nlwmt = nlwm + nlwm;
+ const double sqrtwi = 1.0 / sqrt(2.0);
+ const complex<double> uim(0.0, 1.0);
+ complex<double> cm1 = 0.5 * complex<double>(up[0], up[1]);
+ complex<double> cp1 = 0.5 * complex<double>(up[0], -up[1]);
+ double cz1 = up[2];
+ complex<double> cm2 = 0.5 * complex<double>(un[0], un[1]);
+ complex<double> cp2 = 0.5 * complex<double>(un[0], -un[1]);
+ double cz2 = un[2];
+ for (int l20 = 0; l20 < lw; l20++) {
+ int l = l20 + 1;
+ int lf = l + 1;
+ int lftl = lf * l;
+ double x = 1.0 * lftl;
+ complex<double> cl = complex<double>(four_pi / sqrt(x), 0.0);
+ for (int powi = 1; powi <= l; powi++) cl *= uim;
+ int mv = l + lf;
+ int m = -lf;
+ for (int mf20 = 0; mf20 < mv; mf20++) {
+ m += 1;
+ int k = lftl + m;
+ x = 1.0 * (lftl - m * (m + 1));
+ double cp = sqrt(x);
+ x = 1.0 * (lftl - m * (m - 1));
+ double cm = sqrt(x);
+ double cz = 1.0 * m;
+ c1->w[k - 1][ispo - 1] = dconjg(
+ cp1 * cp * ylm[k + 1] +
+ cm1 * cm * ylm[k - 1] +
+ cz1 * cz * ylm[k]
+ ) * cl;
+ c1->w[k - 1][ispt - 1] = dconjg(
+ cp2 * cp * ylm[k + 1] +
+ cm2 * cm * ylm[k - 1] +
+ cz2 * cz * ylm[k]
+ ) * cl;
+ }
+ }
+ for (int k30 = 0; k30 < nlwm; k30++) {
+ int i = k30 + nlwm;
+ c1->w[i][ispo - 1] = uim * c1->w[k30][ispt - 1];
+ c1->w[i][ispt - 1] = -uim * c1->w[k30][ispo - 1];
+ }
+ if (inpol != 0) {
+ for (int k40 = 0; k40 < nlwm; k40++) {
+ int i = k40 + nlwm;
+ complex<double> cc1 = sqrtwi * (c1->w[k40][ispo - 1] + uim * c1->w[k40][ispt - 1]);
+ complex<double> cc2 = sqrtwi * (c1->w[k40][ispo - 1] - uim * c1->w[k40][ispt - 1]);
+ c1->w[k40][ispo - 1] = cc2;
+ c1->w[i][ispo - 1] = -cc2;
+ c1->w[k40][ispt - 1] = cc1;
+ c1->w[i][ispt - 1] = cc1;
+ }
+ } else {
+ if (isq == 0) {
+ return;
+ }
+ }
+ if (isq != 0) {
+ for (int i50 = 0; i50 < 2; i50++) {
+ int ipt = i50 + 2;
+ int ipis = i50 + is;
+ for (int k50 = 0; k50 < nlwmt; k50++) {
+ c1->w[k50][ipt] = dconjg(c1->w[k50][ipis]);
+ }
+ }
+ }
+}
+
+void rabas(
+ int inpol, int li, int nsph, C1 *c1, double **tqse, complex<double> **tqspe,
+ double **tqss, complex<double> **tqsps
+) {
+ complex<double> cc0 = complex<double>(0.0, 0.0);
+ complex<double> uim = complex<double>(0.0, 1.0);
+ double two_pi = 4.0 * acos(0.0);
+ for (int i80 = 0; i80 < nsph; i80++) {
+ int i = i80 + 1;
+ if(c1->iog[i80] >= i) {
+ tqse[0][i80] = 0.0;
+ tqse[1][i80] = 0.0;
+ tqspe[0][i80] = cc0;
+ tqspe[1][i80] = cc0;
+ tqss[0][i80] = 0.0;
+ tqss[1][i80] = 0.0;
+ tqsps[0][i80] = cc0;
+ tqsps[1][i80] = cc0;
+ for (int l70 = 0; l70 < li; l70++) {
+ int l = l70 + 1;
+ double fl = 1.0 * (l + l + 1);
+ complex<double> rm = 1.0 / c1->rmi[l70][i80];
+ double rmm = (rm * dconjg(rm)).real();
+ complex<double> re = 1.0 / c1->rei[l70][i80];
+ double rem = (re * dconjg(re)).real();
+ if (inpol == 0) {
+ complex<double> pce = ((rm + re) * uim) * fl;
+ complex<double> pcs = ((rmm + rem) * fl) * uim;
+ tqspe[0][i80] -= pce;
+ tqspe[1][i80] += pce;
+ tqsps[0][i80] -= pcs;
+ tqsps[1][i80] += pcs;
+ } else {
+ double ce = (rm + re).real() * fl;
+ double cs = (rmm + rem) * fl;
+ tqse[0][i80] -= ce;
+ tqse[1][i80] += ce;
+ tqss[0][i80] -= cs;
+ tqss[1][i80] += cs;
+ }
+ }
+ if (inpol == 0) {
+ tqspe[0][i80] *= two_pi;
+ tqspe[1][i80] *= two_pi;
+ tqsps[0][i80] *= two_pi;
+ tqsps[1][i80] *= two_pi;
+ } else {
+ tqse[0][i80] *= two_pi;
+ tqse[1][i80] *= two_pi;
+ tqss[0][i80] *= two_pi;
+ tqss[1][i80] *= two_pi;
+ }
+ }
+ }
+}
+
+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 = 1; k <= n; k++)
+ sj[k] = 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 < 0.0) ? -sa : sa;
+ double abs_sb = (sb < 0.0) ? -sb : sb;
+ if (abs_sa > abs_sb) cs = sa / f;
+ else cs = sb / f0;
+ for (int k = 0; k <= nm; k++) {
+ sj[k] = cs * sj[k];
+ }
+}
+
+void rkc(
+ int npntmo, double step, complex<double> dcc, double &x, int lpo,
+ complex<double> &y1, complex<double> &y2, complex<double> &dy1,
+ complex<double> &dy2
+) {
+ complex<double> cy1, cdy1, c11, cy23, yc2, c12, c13;
+ 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 = 0; 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;
+ }
+}
+
+void rkt(
+ int npntmo, double step, double &x, int lpo, complex<double> &y1,
+ complex<double> &y2, complex<double> &dy1, complex<double> &dy2,
+ C2 *c2
+) {
+ complex<double> cy1, cdy1, c11, cy23, cdy23, yc2, c12, c13;
+ 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 = 0; ipnt60 < npntmo; ipnt60++) {
+ int ipnt = ipnt60 + 1;
+ int jpnt = ipnt + ipnt - 1;
+ int jpnt60 = jpnt - 1;
+ cy1 = cl / (x * x) - c2->ris[jpnt60];
+ 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[jpnt];
+ 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[jpntpo];
+ 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[jpnt60];
+ cdy1 += 2.0 * c2->dlri[jpnt60];
+ c21 = (cy1 * y2 + cdy1 * dy2) * step;
+ cy23 -= cdy23 * c2->dlri[jpnt];
+ cdy23 += 2.0 * c2->dlri[jpnt];
+ 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[jpntpo];
+ cdy4 += 2.0 * c2->dlri[jpntpo];
+ 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;
+ }
+}
+
+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;
+}
+
+void sphar(
+ double cosrth, double sinrth, double cosrph, double sinrph,
+ int ll, complex<double> *ylm
+) {
+ const int rmp_size = ll;
+ const int plegn_size = (ll + 1) * ll / 2 + ll + 1;
+ double sinrmp[rmp_size], cosrmp[rmp_size], plegn[plegn_size];
+ double four_pi = 8.0 * acos(0.0);
+ double pi4irs = 1.0 / sqrt(four_pi);
+ double x = cosrth;
+ double y = sinrth;
+ if (y < 0.0) y *= -1.0;
+ double cllmo = 3.0;
+ double cll = 1.5;
+ double ytol = y;
+ plegn[0] = 1.0;
+ plegn[1] = x * sqrt(cllmo);
+ plegn[2] = ytol * sqrt(cll);
+ sinrmp[0] = sinrph;
+ cosrmp[0] = cosrph;
+ if (ll >= 2) {
+ int k = 3;
+ for (int l20 = 2; l20 <= ll; l20++) {
+ int lmo = l20 - 1;
+ int ltpo = l20 + l20 + 1;
+ int ltmo = ltpo - 2;
+ int lts = ltpo * ltmo;
+ double cn = 1.0 * lts;
+ for (int mpo10 = 1; mpo10 <= lmo; mpo10++) {
+ int m = mpo10 - 1;
+ int mpopk = mpo10 + k;
+ int ls = (l20 + m) * (l20 - m);
+ double cd = 1.0 * ls;
+ double cnm = 1.0 * ltpo * (lmo + m) * (l20 - mpo10);
+ double cdm = 1.0 * ls * (ltmo - 2);
+ plegn[mpopk - 1] = plegn[mpopk - l20 - 1] * x * sqrt(cn / cd) -
+ plegn[mpopk - ltmo - 1] * sqrt(cnm / cdm);
+ }
+ int lpk = l20 + k;
+ double cltpo = 1.0 * ltpo;
+ plegn[lpk - 1] = plegn[k - 1] * x * sqrt(cltpo);
+ k = lpk + 1;
+ double clt = 1.0 * (ltpo - 1);
+ cll *= (cltpo / clt);
+ ytol *= y;
+ plegn[k - 1] = ytol * sqrt(cll);
+ sinrmp[l20 - 1] = sinrph * cosrmp[lmo - 1] + cosrph * sinrmp[lmo - 1];
+ cosrmp[l20 - 1] = cosrph * cosrmp[lmo - 1] - sinrph * sinrmp[lmo - 1];
+ } // end l20 loop
+ }
+ // label 30
+ int l = 0;
+ int m, k, l0y, l0p, lmy, lmp;
+ double save;
+ label40:
+ m = 0;
+ k = l * (l + 1);
+ l0y = k + 1;
+ l0p = k / 2 + 1;
+ ylm[l0y - 1] = pi4irs * plegn[l0p - 1];
+ goto label45;
+ label44:
+ lmp = l0p + m;
+ save = pi4irs * plegn[lmp - 1];
+ lmy = l0y + m;
+ ylm[lmy - 1] = save * complex<double>(cosrmp[m - 1], sinrmp[m - 1]);
+ if (m % 2 != 0) ylm[lmy - 1] *= -1.0;
+ lmy = l0y - m;
+ ylm[lmy - 1] = save * complex<double>(cosrmp[m - 1], -sinrmp[m - 1]);
+ label45:
+ if (m >= l) goto label47;
+ m += 1;
+ goto label44;
+ label47:
+ if (l >= ll) return;
+ l += 1;
+ goto label40;
+}
+
+void sscr0(complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
+ complex<double> sum21, rm, re, csam;
+ const complex<double> cc0 = complex<double>(0.0, 0.0);
+ const double exdc = exri * exri;
+ double ccs = 4.0 * acos(0.0) / (vk * vk);
+ double cccs = ccs / exdc;
+ csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ tfsas = cc0;
+ for (int i12 = 0; i12 < nsph; i12++) {
+ int i = i12 + 1;
+ int iogi = c1->iog[i12];
+ if (iogi >= i) {
+ double sums = 0.0;
+ complex<double> sum21 = cc0;
+ for (int l10 = 0; l10 < lm; l10++) {
+ int l = l10 + 1;
+ double fl = 1.0 + l + l;
+ rm = 1.0 / c1->rmi[l10][i12];
+ re = 1.0 / c1->rei[l10][i12];
+ complex<double> rm_cnjg = dconjg(rm);
+ complex<double> re_cnjg = dconjg(re);
+ 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] = scasec;
+ c1->sexs[i12] = extsec;
+ c1->sabs[i12] = abssec;
+ double gcss = c1->gcsv[i12];
+ c1->sqscs[i12] = scasec / gcss;
+ c1->sqexs[i12] = extsec / gcss;
+ c1->sqabs[i12] = abssec / gcss;
+ c1->fsas[i12] = sum21 * csam;
+ }
+ tfsas += c1->fsas[iogi - 1];
+ }
+}
+
+void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
+ complex<double> s11, s21, s12, s22, rm, re, csam, cc;
+ const complex<double> cc0(0.0, 0.0);
+ double ccs = 1.0 / (vk * vk);
+ csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
+ const double pigfsq = 64.0 * acos(0.0) * acos(0.0);
+ double cfsq = 4.0 / (pigfsq * ccs * ccs);
+ int nlmm = lm * (lm + 2);
+ for (int i14 = 0; i14 < nsph; i14++) {
+ int i = i14 + 1;
+ int iogi = c1->iog[i14];
+ if (iogi >= i) {
+ int k = 0;
+ s11 = cc0;
+ s21 = cc0;
+ s12 = cc0;
+ s22 = cc0;
+ for (int l10 = 0; l10 < lm; l10++) {
+ int l = l10 + 1;
+ rm = 1.0 / c1->rmi[l10][i14];
+ re = 1.0 / c1->rei[l10][i14];
+ int ltpo = l + l + 1;
+ for (int im10 = 0; im10 < ltpo; im10++) {
+ k += 1;
+ int ke = k + nlmm;
+ s11 = s11 - c1->w[k - 1][2] * c1->w[k - 1][0] * rm - c1->w[ke - 1][2] * c1->w[ke - 1][0] * re;
+ s21 = s21 - c1->w[k - 1][3] * c1->w[k - 1][0] * rm - c1->w[ke - 1][3] * c1->w[ke - 1][0] * re;
+ s12 = s12 - c1->w[k - 1][2] * c1->w[k - 1][1] * rm - c1->w[ke - 1][2] * c1->w[ke - 1][1] * re;
+ s22 = s22 - c1->w[k - 1][3] * c1->w[k - 1][1] * rm - c1->w[ke - 1][3] * c1->w[ke - 1][1] * re;
+ }
+ }
+ c1->sas[i14][0][0] = s11 * csam;
+ c1->sas[i14][1][0] = s21 * csam;
+ c1->sas[i14][0][1] = s12 * csam;
+ c1->sas[i14][1][1] = s22 * csam;
+ }
+ } // loop i14
+ for (int i24 = 0; i24 < nsph; i24++) {
+ int i = i24 + 1;
+ int iogi = c1->iog[i24];
+ if (iogi >= i) {
+ int j = 0;
+ for (int ipo1 = 0; ipo1 < 2; ipo1++) {
+ for (int jpo1 = 0; jpo1 < 2; jpo1++) {
+ complex<double> cc = dconjg(c1->sas[i24][jpo1][ipo1]);
+ for (int ipo2 = 0; ipo2 < 2; ipo2++) {
+ for (int jpo2 = 0; jpo2 < 2; jpo2++) {
+ c1->vints[i24][j++] = c1->sas[i24][jpo2][ipo2] * cc * cfsq;
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+void thdps(int lm, double ****zpv) {
+ for (int l15 = 0; l15 < lm; l15++) {
+ int l = l15 + 1;
+ double xd = 1.0 * l * (l + 1);
+ double zp = -1.0 / sqrt(xd);
+ zpv[l15][1][0][1] = zp;
+ zpv[l15][1][1][0] = zp;
+ }
+ if (lm != 1) {
+ for (int l20 = 1; l20 < lm; l20++) {
+ int l = l20 + 1;
+ double xn = 1.0 * (l - 1) * (l + 1);
+ double xd = 1.0 * l * (l + l + 1);
+ double zp = sqrt(xn / xd);
+ zpv[l20][0][0][0] = zp;
+ zpv[l20][0][1][1] = zp;
+ }
+ int lmmo = lm - 1;
+ for (int l25 = 0; l25 < lmmo; l25++) {
+ int l = l25 + 1;
+ double xn = 1.0 * l * (l + 2);
+ double xd = (l + 1) * (l + l + 1);
+ double zp = -1.0 * sqrt(xn / xd);
+ zpv[l25][2][0][0] = zp;
+ zpv[l25][2][1][1] = zp;
+ }
+ }
+}
+
+void upvmp(
+ double thd, double phd, int icspnv, double &cost, double &sint,
+ double &cosp, double &sinp, double *u, double *up, double *un
+) {
+ double half_pi = acos(0.0);
+ double rdr = half_pi / 90.0;
+ double th = thd * rdr;
+ double ph = phd * rdr;
+ cost = cos(th);
+ sint = sin(th);
+ cosp = cos(ph);
+ sinp = sin(ph);
+ u[0] = cosp * sint;
+ u[1] = sinp * sint;
+ u[2] = cost;
+ up[0] = cosp * cost;
+ up[1] = sinp * cost;
+ up[2] = -sint;
+ un[0] = -sinp;
+ un[1] = cosp;
+ un[2] = 0.0;
+ if (icspnv != 0) {
+ up[0] *= -1.0;
+ up[1] *= -1.0;
+ up[2] *= -1.0;
+ un[0] *= -1.0;
+ un[1] *= -1.0;
+ }
+}
+
+void upvsp(
+ double *u, double *upmp, double *unmp, double *us, double *upsmp, double *unsmp,
+ double *up, double *un, double *ups, double *uns, double *duk, int &isq,
+ int &ibf, double &scand, double &cfmp, double &sfmp, double &cfsp, double &sfsp
+) {
+ double rdr = acos(0.0) / 90.0;
+ double small = 1.0e-6;
+ isq = 0;
+ scand = u[0] * us[0] + u[1] * us[1] + u[2] * us[2];
+ double abs_scand = (scand >= 1.0) ? scand - 1.0 : 1.0 - scand;
+ if (abs_scand >= small) {
+ abs_scand = scand + 1.0;
+ if (abs_scand < 0.0) abs_scand *= -1.0;
+ if (abs_scand >= small) {
+ scand = acos(scand) / rdr;
+ duk[0] = u[0] - us[0];
+ duk[1] = u[1] - us[1];
+ duk[2] = u[2] - us[2];
+ ibf = 0;
+ } else { // label 15
+ scand = 180.0;
+ duk[0] = 2.0 * u[0];
+ duk[1] = 2.0 * u[1];
+ duk[2] = 2.0 * u[2];
+ ibf = 1;
+ ups[0] = -upsmp[0];
+ ups[1] = -upsmp[1];
+ ups[2] = -upsmp[2];
+ uns[0] = -unsmp[0];
+ uns[1] = -unsmp[1];
+ uns[2] = -unsmp[2];
+ }
+ } else { // label 10
+ scand = 0.0;
+ duk[0] = 0.0;
+ duk[1] = 0.0;
+ duk[2] = 0.0;
+ ibf = -1;
+ isq = -1;
+ ups[0] = upsmp[0];
+ ups[1] = upsmp[1];
+ ups[2] = upsmp[2];
+ uns[0] = unsmp[0];
+ uns[1] = unsmp[1];
+ uns[2] = unsmp[2];
+ }
+ if (ibf == -1 || ibf == 1) { // label 20
+ up[0] = upmp[0];
+ up[1] = upmp[1];
+ up[2] = upmp[2];
+ un[0] = unmp[0];
+ un[1] = unmp[1];
+ un[2] = unmp[2];
+ } else { // label 25
+ orunve(u, us, un, -1, small);
+ uns[0] = un[0];
+ uns[1] = un[1];
+ uns[2] = un[2];
+ orunve(un, u, up, 1, small);
+ orunve(uns, us, ups, 1, small);
+ }
+ // label 85
+ cfmp = upmp[0] * up[0] + upmp[1] * up[1] + upmp[2] * up[2];
+ sfmp = unmp[0] * up[0] + unmp[1] * up[1] + unmp[2] * up[2];
+ cfsp = ups[0] * upsmp[0] + ups[1] * upsmp[1] + ups[2] * upsmp[2];
+ sfsp = uns[0] * upsmp[0] + uns[1] * upsmp[1] + uns[2] * upsmp[2];
+}
+
+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
+) {
+ const int ylm_size = (lm + 1) * (lm + 1) + 1;
+ complex<double> *ylm = new complex<double>[ylm_size];
+ const int nlmp = lm * (lm + 2) + 2;
+ ylm[nlmp - 1] = complex<double>(0.0, 0.0);
+ if (idot != 0) {
+ if (idot != 1) {
+ for (int n40 = 0; n40 < nsph; n40++) {
+ arg[n40] = u[0] * c1->rxx[n40] + u[1] * c1->ryy[n40] + u[2] * c1->rzz[n40];
+ }
+ } else {
+ for (int n50 = 0; n50 < nsph; n50++) {
+ arg[n50] = c1->rzz[n50];
+ }
+ }
+ if (iis == 2) {
+ for (int n60 = 0; n60 < nsph; n60++) arg[n60] *= -1;
+ }
+ }
+ sphar(cost, sint, cosp, sinp, lm, ylm);
+ pwma(up, un, ylm, inpol, lm, iis, c1);
+ delete[] ylm;
+}
+
+void wmasp(
+ double cost, double sint, double cosp, double sinp, double costs, double sints,
+ double cosps, double sinps, double *u, double *up, double *un, double *us,
+ double *ups, double *uns, int isq, int ibf, int inpol, int lm, int idot,
+ int nsph, double *argi, double *args, C1 *c1
+) {
+ const int ylm_size = (lm + 1) * (lm + 1) + 1;
+ complex<double> *ylm = new complex<double>[ylm_size];
+ const int nlmp = lm * (lm + 2) + 2;
+ ylm[nlmp - 1] = complex<double>(0.0, 0.0);
+ if (idot != 0) {
+ if (idot != 1) {
+ for (int n40 = 0; n40 < nsph; n40++) {
+ argi[n40] = u[0] * c1->rxx[n40] + u[1] * c1->ryy[n40] + u[2] * c1->rzz[n40];
+ if (ibf != 0) {
+ args[n40] = argi[n40] * ibf;
+ } else {
+ args[n40] = -1.0 * (us[0] * c1->rxx[n40] + us[1] * c1->ryy[n40] + us[2] * c1->rzz[n40]);
+ }
+ }
+ } else { // label 50
+ for (int n60 = 0; n60 < nsph; n60++) {
+ argi[n60] = cost * c1->rzz[n60];
+ if (ibf != 0) {
+ args[n60] = argi[n60] * ibf;
+ } else {
+ args[n60] = -costs * c1->rzz[n60];
+ }
+ }
+ }
+ }
+ sphar(cost, sint, cosp, sinp, lm, ylm);
+ pwma(up, un, ylm, inpol, lm, isq, c1);
+ if (ibf >= 0) {
+ sphar(costs, sints, cosps, sinps, lm, ylm);
+ pwma(ups, uns, ylm, inpol, lm, 2, c1);
+ }
+ delete[] ylm;
+}
diff --git a/src/sphere/Makefile b/src/sphere/Makefile
index f70b833d..e89f6afd 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
- $(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_sphere $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o
+np_sphere: $(BUILDDIR)/np_sphere.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.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)/sph_subs.o $(BUILDDIR)/sphere.o
$(BUILDDIR)/np_sphere.o:
$(CXX) $(CXXFLAGS) -c np_sphere.cpp -o $(BUILDDIR)/np_sphere.o
@@ -29,6 +29,9 @@ $(BUILDDIR)/Configuration.o:
$(BUILDDIR)/Parsers.o:
$(CXX) $(CXXFLAGS) -c ../libnptm/Parsers.cpp -o $(BUILDDIR)/Parsers.o
+$(BUILDDIR)/sph_subs.o:
+ $(CXX) $(CXXFLAGS) -c ../libnptm/sph_subs.cpp -o $(BUILDDIR)/sph_subs.o
+
$(BUILDDIR)/sphere.o:
$(CXX) $(CXXFLAGS) -c sphere.cpp -o $(BUILDDIR)/sphere.o
--
GitLab