diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index 081b13ac27a1f4d862576b26354d558840b2368e..49bd8dcc7d0a1bf02fcbc11b5098a33d64ba7f21 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -362,6 +362,188 @@ void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
}
}
+/*! \brief C++ porting of ORUNVE
+ *
+ * \param u1: `double *`
+ * \param u2: `double *`
+ * \param u3: `double *`
+ * \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];
+}
+
+/*! \brief C++ porting of PWMA
+ *
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param ylm: Vector of complex
+ * \param inpol: `int`
+ * \param lw: `int`
+ * \param isq: `int`
+ * \param c1: `C1 *`
+ */
+void pwma(
+ double *up, double *un, std::complex<double> *ylm, int inpol, int lw,
+ int isq, C1 *c1
+) {
+ //std::complex<double> cp1, cm1, cc1, cp2, c2, cc2, uim, cl;
+ 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;
+ double sqrtwi = 1.0 / sqrt(2.0);
+ 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 = 1; l20 <= lw; l20++) {
+ int lf = l20 + 1;
+ int lftl = lf * l20;
+ double x = 1.0 * lftl;
+ std::complex<double> cl = (four_pi / sqrt(x)) * std::pow(uim, l20);
+ int mv = l20 + lf;
+ int m = -lf;
+ for (int mf20 = 1; 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 = 1; k30 <= nlwm; k30++) {
+ int i = k30 + nlwm;
+ c1->w[i - 1][ispo - 1] = uim * c1->w[k30 - 1][ispt - 1];
+ c1->w[i - 1][ispt - 1] = -uim * c1->w[k30 - 1][ispo - 1];
+ }
+ if (inpol != 0) {
+ for (int k40 = 1; k40 <= nlwm; k40++) {
+ int i = k40 + nlwm;
+ std::complex<double> cc1 = sqrtwi * (c1->w[k40 - 1][ispo - 1] + uim * c1->w[k40 - 1][ispt - 1]);
+ std::complex<double> cc2 = sqrtwi * (c1->w[k40 - 1][ispo - 1] - uim * c1->w[k40 - 1][ispt - 1]);
+ c1->w[k40 - 1][ispo - 1] = cc2;
+ c1->w[i - 1][ispo - 1] = -cc2;
+ c1->w[k40 - 1][ispt - 1] = cc1;
+ c1->w[i - 1][ispt - 1] = cc1;
+ }
+ }
+ if (isq != 0) {
+ for (int i50 = 1; i50 <= 2; i50++) {
+ int ipt = i50 + 2;
+ int ipis = i50 + is;
+ for (int k50 = 1; k50 <= nlwmt; k50++)
+ c1->w[k50 - 1][ipt - 1] = dconjg(c1->w[k50 - 1][ipis - 1]);
+ }
+ }
+ }
+}
+
+/*! \brief C++ porting of RABAS
+ *
+ * \param inpol: `int`
+ * \param li: `int`
+ * \param nsph: `int`
+ * \param c1: `C1 *`
+ * \param TQSE: Matrix of complex.
+ * \param TQSPE: Matrix of complex.
+ * \param TQSS: Matrix of complex.
+ * \param TQSPS: Matrix of complex.
+ */
+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 = 1; i80 <= nsph; i80++) {
+ if(c1->iog[i80 - 1] >= i80) {
+ tqse[0][i80 - 1] = 0.0;
+ tqse[1][i80 - 1] = 0.0;
+ tqspe[0][i80 - 1] = cc0;
+ tqspe[1][i80 - 1] = cc0;
+ tqss[0][i80 - 1] = 0.0;
+ tqss[1][i80 - 1] = 0.0;
+ tqsps[0][i80 - 1] = cc0;
+ tqsps[1][i80 - 1] = cc0;
+ for (int l70 = 1; l70 <= li; l70++) {
+ double fl = 1.0 * l70 + l70 + 1;
+ std::complex<double> rm = 1.0 / c1->rmi[l70 - 1][i80 - 1];
+ double rmm = (rm * dconjg(rm)).real();
+ std::complex<double> re = 1.0 / c1->rei[l70 - 1][i80 - 1];
+ 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 - 1] -= pce;
+ tqspe[1][i80 - 1] += pce;
+ tqsps[0][i80 - 1] -= pcs;
+ tqsps[1][i80 - 1] += pcs;
+ } else {
+ double ce = (rm + re).real() * fl;
+ double cs = (rmm + rem) * fl;
+ tqse[0][i80 - 1] -= ce;
+ tqse[1][i80 - 1] += ce;
+ tqss[0][i80 - 1] -= cs;
+ tqss[1][i80 - 1] += cs;
+ }
+ }
+ if (inpol == 0) {
+ tqspe[0][i80 - 1] *= two_pi;
+ tqspe[1][i80 - 1] *= two_pi;
+ tqsps[0][i80 - 1] *= two_pi;
+ tqsps[1][i80 - 1] *= two_pi;
+ printf("DEBUG: TQSPE(1, %d ) = ( %lE, %lE)\n", i80, tqspe[0][i80 - 1].real(), tqspe[0][i80 - 1].imag());
+ printf("DEBUG: TQSPE(2, %d ) = ( %lE, %lE)\n", i80, tqspe[1][i80 - 1].real(), tqspe[1][i80 - 1].imag());
+ printf("DEBUG: TQSPS(1, %d ) = ( %lE, %lE)\n", i80, tqsps[0][i80 - 1].real(), tqsps[0][i80 - 1].imag());
+ printf("DEBUG: TQSPS(2, %d ) = ( %lE, %lE)\n", i80, tqsps[1][i80 - 1].real(), tqsps[1][i80 - 1].imag());
+ } else {
+ tqse[0][i80 - 1] *= two_pi;
+ tqse[1][i80 - 1] *= two_pi;
+ tqss[0][i80 - 1] *= two_pi;
+ tqss[1][i80 - 1] *= two_pi;
+ printf("DEBUG: TQSE(1, %d ) = %lE)\n", i80, tqse[0][i80 - 1]);
+ printf("DEBUG: TQSE(2, %d ) = %lE)\n", i80, tqse[1][i80 - 1]);
+ printf("DEBUG: TQSS(1, %d ) = %lE)\n", i80, tqss[0][i80 - 1]);
+ printf("DEBUG: TQSS(2, %d ) = %lE)\n", i80, tqss[1][i80 - 1]);
+ }
+ }
+ }
+}
+
/*! \brief C++ porting of RBF
*
* This is the Real Bessel Function.
@@ -639,6 +821,95 @@ void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri
}
}
+/*! \brief C++ porting of SPHAR
+ *
+ * \param cosrth: `double`
+ * \param sinrth: `double`
+ * \param cosrph: `double`
+ * \param sinrph: `double`
+ * \param ll: `int`
+ * \param ylm: Vector of complex.
+ */
+void sphar(
+ double cosrth, double sinrth, double cosrph, double sinrph,
+ int ll, std::complex<double> *ylm
+) {
+ double sinrmp[40], cosrmp[40], plegn[861];
+ double four_pi = 8.0 * acos(0.0);
+ double pi4irs = 1.0 / sqrt(four_pi);
+ double x = cosrth;
+ double y = sinrth;
+ //printf("DEBUG: X = %lE\n", x);
+ if (y < 0.0) y *= -1.0;
+ //printf("DEBUG: Y = %lE\n", y);
+ double cllmo = 3.0;
+ double cll = 1.5;
+ double ytol = y;
+ plegn[0] = 1.0;
+ //printf("DEBUG: PLEGN( %d ) = %lE\n", 1, plegn[0]);
+ plegn[1] = x * sqrt(cllmo);
+ //printf("DEBUG: PLEGN( %d ) = %lE\n", 2, plegn[1]);
+ plegn[2] = ytol * sqrt(cll);
+ //printf("DEBUG: PLEGN( %d ) = %lE\n", 3, plegn[2]);
+ sinrmp[0] = sinrph;
+ cosrmp[0] = cosrph;
+ int k;
+ if (ll >= 2) {
+ 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);
+ //printf("DEBUG: PLEGN( %d ) = %lE\n", mpopk, plegn[mpopk - 1]);
+ }
+ int lpk = l20 + k;
+ double cltpo = 1.0 * ltpo;
+ plegn[lpk - 1] = plegn[k - 1] * x * sqrt(cltpo);
+ //printf("DEBUG: PLEGN( %d ) = %lE\n", lpk, plegn[lpk - 1]);
+ k = lpk + 1;
+ double clt = 1.0 * (ltpo - 1);
+ cll *= (cltpo / clt);
+ ytol *= y;
+ plegn[k - 1] = ytol * sqrt(cll);
+ //printf("DEBUG: PLEGN( %d ) = %lE\n", k, plegn[k - 1]);
+ 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
+ }
+ //l = 0;
+ //bool goto40 = true;
+ for (int l = 0; l < ll; l++) {
+ //int m = 0;
+ k = l * (l + 1);
+ int l0y = k + 1;
+ int l0p = k / 2 + 1;
+ ylm[l0y - 1] = pi4irs * plegn[l0p - 1];
+ //printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", l0y, ylm[l0y - 1].real(), ylm[l0y - 1].imag());
+ for (int m = 0; m < l; m++) {
+ int lmp = l0p + m;
+ double save = pi4irs * plegn[lmp - 1];
+ int 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;
+ //printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", lmy, ylm[lmy - 1].real(), ylm[lmy - 1].imag());
+ lmy = l0y - m;
+ ylm[lmy - 1] = save * std::complex<double>(cosrmp[m - 1], -sinrmp[m - 1]);
+ //printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", lmy, ylm[lmy - 1].real(), ylm[lmy - 1].imag());
+ }
+ }
+}
+
/*! \brief C++ porting of THDPS
*
* \param lm: `int`
@@ -680,6 +951,255 @@ void thdps(int lm, double ****zpv) {
}
}
+/*! \brief C++ porting of UPVMP
+ *
+ * \param thd: `double`
+ * \param phd: `double`
+ * \param icspnv: `int`
+ * \param cost: `double`
+ * \param sint: `double`
+ * \param cosp: `double`
+ * \param sinp: `double`
+ * \param u: `double *`
+ * \param up: `double *`
+ * \param un: `double *`
+ */
+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);
+ //printf("DEBUG: cost = %lE\n", cost);
+ //printf("DEBUG: sint = %lE\n", sint);
+ //printf("DEBUG: cosp = %lE\n", cosp);
+ //printf("DEBUG: sinp = %lE\n", sinp);
+ 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 C++ porting of UPVSP
+ *
+ * \param u: `double *`
+ * \param upmp: `double *`
+ * \param unmp: `double *`
+ * \param us: `double *`
+ * \param upsmp: `double *`
+ * \param unsmp: `double *`
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param ups: `double *`
+ * \param uns: `double *`
+ * \param duk: `double *`
+ * \param isq: `int &`
+ * \param ibf: `int &`
+ * \param scand: `double &`
+ * \param cfmp: `double &`
+ * \param sfmp: `double &`
+ * \param cfsp: `double &`
+ * \param sfsp: `double &`
+ */
+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;
+ if (abs_scand < 0.0) abs_scand *= -1.0;
+ 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 C++ porting of WMAMP
+ *
+ * \param iis: `int`
+ * \param cost: `double`
+ * \param sint: `double`
+ * \param cosp: `double`
+ * \param sinp: `double`
+ * \param inpol: `int`
+ * \param lm: `int`
+ * \param idot: `int`
+ * \param nsph: `int`
+ * \param arg: `double *`
+ * \param u: `double *`
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param c1: `C1 *`
+ */
+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
+) {
+ std::complex<double> *ylm = new std::complex<double>[1682];
+ int nlmp = lm * (lm + 2) + 2;
+ ylm[nlmp] = 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 C++ porting of WMASP
+ *
+ * \param cost: `double`
+ * \param sint: `double`
+ * \param cosp: `double`
+ * \param sinp: `double`
+ * \param costs: `double`
+ * \param sints: `double`
+ * \param cosps: `double`
+ * \param sinps: `double`
+ * \param u: `double *`
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param us: `double *`
+ * \param ups: `double *`
+ * \param uns: `double *`
+ * \param isq: `int`
+ * \param ibf: `int`
+ * \param inpol: `int`
+ * \param lm: `int`
+ * \param idot: `int`
+ * \param nsph: `int`
+ * \param argi: `double *`
+ * \param args: `double *`
+ * \param c1: `C1 *`
+ */
+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
+) {
+ std::complex<double> *ylm = new std::complex<double>[1682];
+ 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 = 1; n40 <= nsph; n40++) {
+ argi[n40 - 1] = u[0] * c1->rxx[n40 - 1] + u[1] * c1->ryy[n40 - 1] + u[2] * c1->rzz[n40 - 1];
+ if (ibf != 0) {
+ args[n40 - 1] = argi[n40 - 1] * ibf;
+ } else {
+ args[n40 - 1] = -1.0 * (us[0] * c1->rxx[n40 - 1] + us[1] * c1->ryy[n40 - 1] + us[2] * c1->rzz[n40 - 1]);
+ }
+ }
+ } else { // label 50
+ for (int n60 = 1; n60 <= nsph; n60++) {
+ argi[n60 - 1] = cost * c1->rzz[n60 - 1];
+ if (ibf != 0) {
+ args[n60 - 1] = argi[n60 - 1] * ibf;
+ } else {
+ args[n60 - 1] = -costs * c1->rzz[n60 - 1];
+ }
+ }
+ }
+ }
+ 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 C++ porting of DME
*
* \param li: `int`
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index e422603dcbfe03776ae24f5a23e223a6831ffb00..d1c9ffb1ac742eb87428ab15ca6a5d17c6a56d6d 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -10,7 +10,10 @@ using namespace std;
//! \brief C++ implementation of SPH
void sphere() {
complex<double> arg, s0, tfsas;
- double gaps[2];
+ complex<double> **tqspe, **tqsps;
+ double **tqse, **tqss;
+ double *argi, *args, *gaps;
+ double th, ph;
printf("INFO: making legacy configuration ...\n");
ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/sphere/DEDFB");
conf->write_formatted("c_OEDFB");
@@ -34,6 +37,9 @@ void sphere() {
c1->rei[i][1] = complex<double>(0.0, 0.0);
}
C2 *c2 = new C2;
+ argi = new double[1];
+ args = new double[1];
+ gaps = new double[2];
FILE *output = fopen("c_OSPH", "w");
fprintf(output, " READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM\n");
fprintf(
@@ -265,8 +271,217 @@ void sphere() {
printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag());
double sqk = vk * vk * sconf->exdc;
aps(zpv, gconf->l_max, nsph, c1, sqk, gaps);
- // Would call RABAS
- } //jxi488
+ tqse = new double*[2];
+ tqss = new double*[2];
+ tqspe = new std::complex<double>*[2];
+ tqsps = new std::complex<double>*[2];
+ for (int ti = 0; ti < 2; ti++) {
+ tqse[ti] = new double[2];
+ tqss[ti] = new double[2];
+ tqspe[ti] = new std::complex<double>[2];
+ tqsps[ti] = new std::complex<double>[2];
+ }
+ rabas(gconf->in_pol, gconf->l_max, nsph, c1, tqse, tqspe, tqss, tqsps);
+ for (int i170 = 1; i170 <= nsph; i170++) {
+ if (c1->iog[i170 - 1] >= i170) {
+ double albeds = c1->sscs[i170 - 1] / c1->sexs[i170 - 1];
+ c1->sqscs[i170 - 1] *= sqsfi;
+ c1->sqabs[i170 - 1] *= sqsfi;
+ c1->sqexs[i170 - 1] *= sqsfi;
+ fprintf(output, " SPHERE %2d\n", i170);
+ if (c1->nshl[i170 - 1] != 1) {
+ fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i170 - 1]);
+ } else {
+ fprintf(
+ output,
+ " SIZE=%15.7lE, REFRACTIVE INDEX=(%15.7lE,%15.7lE)\n",
+ c2->vsz[i170 -1],
+ c2->vkt[i170 - 1].real(),
+ c2->vkt[i170 - 1].imag()
+ );
+ }
+ fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
+ fprintf(
+ output,
+ " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+ c1->sscs[i170 - 1], c1->sabs[i170 - 1],
+ c1->sexs[i170 - 1], albeds
+ );
+ fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
+ fprintf(
+ output,
+ " %14.7lE%15.7lE%15.7lE\n",
+ c1->sqscs[i170 - 1], c1->sqabs[i170 - 1],
+ c1->sqexs[i170 - 1]
+ );
+ fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i170 - 1].real(), c1->fsas[i170 - 1].imag());
+ double csch = 2.0 * vk * sqsfi / c1->gcsv[i170 -1];
+ s0 = c1->fsas[i170 - 1] * exri;
+ double qschu = csch * s0.imag();
+ double pschu = csch * s0.real();
+ double s0mag = cs0 * sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag()));
+ fprintf(
+ output,
+ " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ qschu, pschu, s0mag
+ );
+ double rapr = c1->sexs[i170 - 1] - gaps[i170 - 1];
+ double cosav = gaps[i170 - 1] / c1->sscs[i170 - 1];
+ fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
+ fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 1, tqse[0][i170 - 1], tqss[0][i170 - 1]);
+ fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170 - 1], tqss[1][i170 - 1]);
+ tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170 - 1])), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170 - 1])), sizeof(double));
+ double val = tqspe[0][i170 - 1].real();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqspe[0][i170 - 1].imag();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqsps[0][i170 - 1].real();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqsps[0][i170 - 1].imag();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&(tqse[1][i170 - 1])), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&(tqss[1][i170 - 1])), sizeof(double));
+ val = tqspe[1][i170 - 1].real();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqspe[1][i170 - 1].imag();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqsps[1][i170 - 1].real();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ val = tqsps[1][i170 - 1].imag();
+ tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
+ } // End if iog[i170 - 1] >= i170
+ } // i170 loop
+ if (nsph != 1) {
+ fprintf(output, " FSAT=(%15.7lE,%15.7lE)\n", tfsas.real(), tfsas.imag());
+ double csch = 2.0 * vk * sqsfi / gcs;
+ s0 = tfsas * exri;
+ double qschu = csch * s0.imag();
+ double pschu = csch * s0.real();
+ double s0mag = cs0 * sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag()));
+ fprintf(
+ output,
+ " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ qschu, pschu, s0mag
+ );
+ }
+ th = gconf->in_theta_start;
+ int isq, ibf;
+ double cost, sint, cosp, sinp;
+ double costs, sints, cosps, sinps;
+ double scan;
+ double *duk = new double[3];
+ double *u = new double[3];
+ double *us = new double[3];
+ double *un = new double[3];
+ double *uns = new double[3];
+ double *up = new double[3];
+ double *ups = new double[3];
+ double *upmp = new double[3];
+ double *upsmp = new double[3];
+ double *unmp = new double[3];
+ double *unsmp = new double[3];
+ double frx, fry, frz;
+ double cfmp, cfsp, sfmp, sfsp;
+ int jw;
+ for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP parallelizable section
+ ph = gconf->in_phi_start;
+ for (int jph484 = 1; jph484 <= nph; jph484++) {
+ bool goto182 = (nk == 1) && (jxi488 > 1);
+ if (!goto182) {
+ upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
+ }
+ if (gconf->meridional_type >= 0) {
+ wmamp(0, cost, sint, cosp, sinp, gconf->in_pol, gconf->l_max, 0, nsph, argi, u, upmp, unmp, c1);
+ for (int i183 = 0; i183 < nsph; i183++) {
+ double rapr = c1->sexs[i183] - gaps[i183];
+ frx = rapr * u[0];
+ fry = rapr * u[1];
+ frx = rapr * u[2];
+ }
+ jw = 1;
+ }
+ double thsl = ths1;
+ double phsph = 0.0;
+ for (int jths482 = 1; jths482 <= nths; jths482++) {
+ double ths = thsl;
+ int icspnv = 0;
+ if (gconf->meridional_type > 1) ths = th + thsca;
+ if (gconf->meridional_type >= 1) {
+ phsph = 0.0;
+ if ((ths < 0.0) || (ths > 180.0)) phsph = 180.0;
+ if (ths < 0.0) ths *= -1.0;
+ if (ths > 180.0) ths = 360.0 - ths;
+ if (phsph != 0.0) icspnv = 1;
+ }
+ double phs = phs1;
+ for (int jphs480 = 1; jphs480 <= nphs; jphs480++) {
+ if (gconf->meridional_type >= 1) {
+ phs = ph + phsph;
+ if (phs >= 360.0) phs -= 360.0;
+ }
+ bool goto190 = ((nks == 1) && (jxi488 > 1)) || (jth486 > 1) || (jph484 > 1);
+ if (!goto190) {
+ upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
+ if (gconf->meridional_type >= 0) {
+ wmamp(2, costs, sints, cosps, sinps, gconf->in_pol, gconf->l_max, 0, nsph, args, us, upsmp, unsmp, c1);
+ }
+ }
+ if (nkks != 0 || jxi488 == 1) {
+ upvsp(u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns, duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp);
+ if (gconf->meridional_type < 0) {
+ wmasp(
+ cost, sint, cosp, sinp, costs, sints, cosps, sinps,
+ u, up, un, us, ups, uns, isq, ibf, gconf->in_pol,
+ gconf->l_max, 0, nsph, argi, args, c1
+ );
+ }
+ for (int i193 = 1; i193 <= 3; i193++) {
+ un[i193 - 1] = unmp[i193 - 1];
+ uns[i193 - 1] = unsmp[i193 - 1];
+ }
+ }
+ if (gconf->meridional_type < 0) jw = 1;
+ tppoan.write(reinterpret_cast<char *>(&th), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&ph), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&ths), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&phs), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&scan), sizeof(double));
+ if (jw != 0) {
+ jw = 0;
+ tppoan.write(reinterpret_cast<char *>(&(u[0])), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&(u[1])), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&(u[2])), sizeof(double));
+ }
+ fprintf(
+ output,
+ "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n",
+ jth486, jph484, jths482, jphs480
+ );
+ fprintf(
+ output,
+ " TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n",
+ th, ph, ths, phs
+ );
+ fprintf(output, " SCAND=%10.3lE\n", scan);
+ fprintf(output, " CFMP=%15.7lE, SFMP=%15.7lE\n", cfmp, sfmp);
+ fprintf(output, " CFSP=%15.7lE, SFSP=%15.7lE\n", cfsp, sfsp);
+ if (gconf->meridional_type >= 0) {
+ fprintf(output, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
+ fprintf(output, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n", uns[0], uns[1], uns[2]);
+ } else {
+ fprintf(output, " UN=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
+ }
+
+ if (gconf->meridional_type < 1) phs += gconf->sc_phi_step;
+ } // jphs480 loop
+ if (gconf->meridional_type <= 1) thsl += gconf->sc_theta_step;
+ } // jths482 loop
+ ph += gconf->in_phi_step;
+ } // jph484 loop on elevation
+ th += gconf->in_theta_step;
+ } // jth486 loop on azimuth
+ } //jxi488 loop on scales
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
printf("ERROR: failed to open TPPOAN file.\n");