diff --git a/src/cluster/Makefile b/src/cluster/Makefile
index 181647521ebfdeadc696de00957a8d471153e6c3..ae2539bf7cdf45073f9af17ea6c78bf0cbbdef37 100644
--- a/src/cluster/Makefile
+++ b/src/cluster/Makefile
@@ -2,7 +2,6 @@ BUILDDIR=../../build/cluster
FC=gfortran
FCFLAGS=-std=legacy -O3
LFLAGS=
-LFLAGS=
CXX=g++
CXXFLAGS=-O2 -ggdb -pg -coverage
CXXLFLAGS=
diff --git a/src/include/Parsers.h b/src/include/Parsers.h
index a4523de6d03c78d8efe2467e5e37b067ddcbe1e3..94dd68e891b175c7263ca45585625e96122e0bf1 100644
--- a/src/include/Parsers.h
+++ b/src/include/Parsers.h
@@ -1,14 +1,14 @@
/*! \file Parsers.h
*/
-#ifndef INCLUDE_PARSERS_H_
-#define INCLUDE_PARSERS_H_
-
#ifndef FILE_NOT_FOUND_ERROR
//! Error code if a file is not found.
#define FILE_NOT_FOUND_ERROR 21
#endif
+#ifndef INCLUDE_PARSERS_H_
+#define INCLUDE_PARSERS_H_
+
/*! \brief Load a text file as a sequence of strings in memory.
*
* The configuration of the field expansion code in FORTRAN uses
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 243918001a6d821e7ecfd9ab81b90b769556d760..9351e6eabaf81dc463d4697587d633ed770a2b03 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -2,7 +2,6 @@
*
* \brief C++ porting of CLU functions and subroutines.
*
- *
* This library includes a collection of functions that are used to solve the
* scattering problem in the case of a cluster of spheres. The functions that
* were generalized from the case of the single sphere are imported the `sph_subs.h`
diff --git a/src/include/file_io.h b/src/include/file_io.h
deleted file mode 100644
index 24f2c5f0cff1165afacbea0b5d2ac443fd0f6df7..0000000000000000000000000000000000000000
--- a/src/include/file_io.h
+++ /dev/null
@@ -1,43 +0,0 @@
-/*! \file file_io.h
- *
- * C++ wrapper of FORTRAN I/O operations with files.
- */
-
-/*! \brief Open a file for subsequent access.
- *
- * \param uid: `int*` Pointer to the unit ID to be associated to the file.
- * \param name: `char*` C-string for file name (max. length of 63).
- * \param sta: `char*` C-string for file status (max. length of 7).
- * \param mode: `char*` C-string for access mode (max. length of 11).
- */
-extern "C" void open_file_(int* uid, const char* name, const char* sta, const char* mode);
-/*! \brief Close a previously opened file.
- *
- * \param uid: `int*` Pointer to the unit ID of the file.
- */
-extern "C" void close_file_(int* uid);
-/*! \brief Read an integer value from a file.
- *
- * \param uid: `int*` Pointer to the unit ID of the file.
- * \param value: `int*` Pointer of the variable to be updated.
- */
-extern "C" void read_int_(int* uid, int* value);
-/*! \brief Write a complex value to a file.
- *
- * \param uid: `int*` Pointer to the unit ID of the file.
- * \param real: `double*` Pointer to the real part of the value.
- * \param imag: `double*` Pointer to the imaginary part of the value.
- */
-extern "C" void write_complex_(int* uid, double* real, double* imag);
-/*! \brief Write a double precision float value to a file.
- *
- * \param uid: `int*` Pointer to the unit ID of the file.
- * \param value: `double*` Pointer to the variable to be written.
- */
-extern "C" void write_double_(int* uid, double* value);
-/*! \brief Write an integer value to a file.
- *
- * \param uid: `int*` Pointer to the unit ID of the file.
- * \param value: `int*` Pointer to the variable to be written.
- */
-extern "C" void write_int_(int* uid, int* value);
diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index 803670c18a2f66a2dca4c7fa2352c75fc76f2580..1478812e9b7394c4c92adab4d9057dd054012616 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -1,12 +1,30 @@
-#include <fstream>
-#include <cmath>
-#include <complex>
+/*! \file tra_subs.h
+ *
+ * \brief C++ porting of TRAPPING functions and subroutines.
+ *
+ * This library includes a collection of functions that are used to solve the
+ * trapping problem. The functions that were generalized from the case of the
+ * single sphere are imported the `sph_subs.h` library. As it occurs with the
+ * single sphere case functions, in most cases, the results of calculations do
+ * not fall back to fundamental data types. They are rather multi-component
+ * structures. In order to manage access to such variety of return values, most
+ * functions are declared as `void` and they operate on output arguments passed
+ * by reference.
+ */
#ifndef INCLUDE_TRA_SUBS_H_
#define INCLUDE_TRA_SUBS_H_
#endif
-// Structures for TRAPPING (that will probably move to Commons)
+#ifndef INCLUDE_SPH_SUBS_H_
+#include "../include/sph_subs.h"
+#endif
+
+#include <fstream>
+#include <cmath>
+#include <complex>
+
+// Structures for TRAPPING
/*! \brief CIL data structure.
*
* A structure containing field expansion order configuration.
@@ -36,278 +54,6 @@ struct CCR {
};
//End of TRAPPING structures
-//Externally defined subroutines
-extern double cg1(int lmpml, int mu, int l, int m);
-extern std::complex<double> dconjg(std::complex<double> z);
-extern void sphar(
- double cosrth, double sinrth, double cosrph, double sinrph,
- int ll, std::complex<double> *ylm
-);
-extern void thdps(int lm, double ****zpv);
-//End of externally defined subroutines
-
-/*! C++ porting of PWMALP
- *
- * \param w: Matrix of complex. QUESTION: definition?
- * \param up: `double *`
- * \param un: `double *`
- * \param ylm: Vector of complex
- * \param lw: `int`
- */
-void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw) {
- std::complex<double> cp1, cm1, cp2, cm2, cl;
- const std::complex<double> uim(0.0, 1.0);
- const double four_pi = acos(0.0) * 8.0;
- const int nlwm = lw * (lw + 2);
- cm1 = 0.5 * std::complex<double>(up[0], up[1]);
- cp1 = 0.5 * std::complex<double>(up[0], -up[1]);
- double cz1 = up[2];
- cm2 = 0.5 * std::complex<double>(un[0], un[1]);
- 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 = std::pow((four_pi / sqrt(x)) * uim, 1.0 * l20);
- int mv = l20 + lf;
- int m = -lf;
- for (int mf20 = 1; mf20 <= mv; mf20++) {
- m++;
- 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;
- w[k - 1][0] = dconjg(cp1 * cp * ylm[k + 1] + cm1 * cm * ylm[k - 1] + cz1 * cz * ylm[k]) * cl;
- w[k - 1][1] = dconjg(cp2 * cp * ylm[k + 1] + cm2 * cm * ylm[k - 1] + cz2 * cz * ylm[k]) * cl;
- } // mf20 loop
- } // l20 loop
- for (int k30 = 0; k30 < nlwm; k30++) {
- int i = k30 + nlwm;
- w[i][0] = uim * w[k30][1];
- w[i][1] = -uim * w[k30][0];
- } // k30 loop
-}
-
-/*! C++ porting of WAMFF
- *
- * \param wk: Vector of complex. QUESTION: definition?
- * \param x: `double`
- * \param y: `double`
- * \param z: `double`
- * \param lm: `int`
- * \param apfafa: `double` QUESTION: definition?
- * \param tra: `double` QUESTION: definition?
- * \param spd: `double` QUESTION: definition?
- * \param rir: `double` QUESTION: definition?
- * \param ftcn: `double` QUESTION: definition?
- * \param lmode: `int` QUESTION: definition?
- * \param pmf: `double` QUESTION: definition?
- */
-void wamff(
- std::complex<double> *wk, double x, double y, double z, int lm, double apfafa,
- double tra, double spd, double rir, double ftcn, int lmode, double pmf
-) {
- const int nlmm = lm * (lm + 2);
- const int nlmmt = 2 * nlmm;
- const int nlmp = nlmm + 2;
- std::complex<double> **w, *ylm;
- const std::complex<double> cc0(0.0, 0.0);
- const std::complex<double> uim(0.0, 1.0);
- std::complex<double> cfam, cf1, cf2;
- double rho, cph, sph, cth, sth, r;
- double ftc1, ftc2;
- double *up = new double[3];
- double *un = new double[3];
- w = new std::complex<double>*[nlmmt];
- for (int wi = 0; wi < nlmmt; wi++) w[wi] = new std::complex<double>[2]();
- ylm = new std::complex<double>[nlmp]();
- bool onx = (y == 0.0);
- bool ony = (x == 0.0);
- bool onz = (onx && ony);
- if (!(onz && onx && ony)) {
- rho = sqrt(x * x + y * y);
- cph = x / rho;
- sph = y = rho;
- } else {
- if (onz) { // label 10
- cph = 1.0;
- sph = 0.0;
- } else {
- if (onx) { // label 12
- rho = sqrt(x * x);
- cph = (x < 0.0)? -1.0 : 1.0;
- sph = 0.0;
- } else {
- if (ony) { // label 13
- rho = sqrt(y * y);
- cph = 0.0;
- sph = (y < 0.0)? -1.0: 1.0;
- }
- }
- }
- }
- // label 15
- if (z == 0.0) {
- if (!onz) {
- r = rho;
- cth = 0.0;
- sth = 1.0;
- } else { // label 17
- r = 0.0;
- cth = 1.0;
- sth = 0.0;
- }
- } else { // label 18
- if (!onz) {
- r = sqrt(rho * rho + z * z);
- cth = z / r;
- sth = rho / r;
- } else { // label 20
- r = sqrt(z * z);
- cth = (z < 0.0)? -1.0: 1.0;
- sth = 0.0;
- }
- }
- if (lmode == 0 || sth != 0.0) { // label 25
- ylm[nlmp - 1] = cc0;
- sphar(cth, sth, cph, sph, lm, ylm);
- up[0] = cph * cth;
- up[1] = sph * cth;
- up[2] = -sth;
- un[0] = -sph;
- un[1] = cph;
- un[2] = 0.0;
- // Would call PWMALP(W,UP,UN,YLM,LM)
- pwmalp(w, up, un, ylm, lm);
- double apfa = sth * apfafa;
- if (spd > 0.0) {
- double sthl = sth * rir;
- double cthl = sqrt(1.0 - sthl * sthl);
- double arg = spd * (z - (r / rir) * cthl);
- cfam = (tra * std::exp(-apfa * apfa) / sqrt(cthl)) * std::exp(uim * arg);
- if (lmode == 0) { // CHECK: Computed GO TO, not sure what happens with LMODE = 0
- if (sth == 0.0) { // label 45
- ftc1 = ftcn;
- ftc2 = ftcn;
- // goes to 48
- }
- } else if (lmode == 1) { // label 46
- cfam *= ((cph + uim * sph) * sth * pmf);
- ftc1 = 2.0 * cthl / (cthl * rir + cth);
- ftc2 = 2.0 * cthl / (cthl + cth * rir);
- // follows on to 48
- } else if (lmode == 2) { // label 50
- ftc1 = 2.0 * cthl / (cthl * rir + cth);
- cfam *= (sth * pmf * ftc1);
- for (int i52 = 0; i52 < nlmmt; i52++) wk[i52] = w[i52][0] * cfam;
- // returns
- } else if (lmode == 3) { // label 53
- ftc2 = 2.0 * cthl / (cthl + cth * rir);
- cfam *= (sth * pmf * ftc2);
- for (int i55 = 0; i55 < nlmmt; i55++) wk[i55] = w[i55][1] * cfam;
- // returns
- }
- if (lmode == 0 || lmode == 1) { //label 48
- cf1 = cph * ftc1 * cfam;
- cf2 = -sph * ftc2 * cfam;
- // goes to 62
- }
- } else { // label 57
- double fam = tra * std::exp(-apfa * apfa) / sqrt(cth);
- if (lmode == 0) {
- double f1 = cph * fam;
- double f2 = -sph * fam;
- for (int i58 = 0; i58 < nlmmt; i58++) wk[i58] = w[i58][0] * f1 + w[i58][1] * f2;
- // returns
- } else if (lmode == 1) { // label 60
- cfam = (pmf * sth * fam) * (cph * uim * sph);
- cf1 = cph * cfam;
- cf2 = -sph * cfam;
- // follows on to 62
- } else if (lmode == 2) { // label 65
- fam *= (pmf * sth);
- for (int i67 = 0; i67 < nlmmt; i67++) wk[i67] = w[i67][0] * fam;
- // returns
- } else if (lmode == 3) { // label 68
- fam *= (pmf * sth);
- for (int i70 = 0; i70 < nlmmt; i70++) wk[i70] = w[i70][1] * fam;
- // returns
- }
- }
- if (lmode == 0 || lmode == 1) { // label 62
- for (int i63 = 0; i63 < nlmmt; i63++) wk[i63] = w[i63][0] * cf1 + w[i63][1] * cf2;
- }
- }
- // Clean up memory
- delete[] up;
- delete[] un;
- for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
- delete[] w;
- delete[] ylm;
-}
-
-/*! C++ porting of FRFMER
- *
- * \param nkv: `int` QUESTION: definition?
- * \param vkm: `double` QUESTION: definition?
- * \param vkv: `double *` QUESTION: definition?
- * \param vknmx: `double` QUESTION: definition?
- * \param apfafa: `double` QUESTION: definition?
- * \param tra: `double` QUESTION: definition?
- * \param spd: `double` QUESTION: definition?
- * \param rir: `double` QUESTION: definition?
- * \param ftcn: `double` QUESTION: definition?
- * \param le: `int` QUESTION: definition?
- * \param lmode: `int` QUESTION: definition?
- * \param pmf: `double` QUESTION: definition?
- * \param tt1: `fstream &` Handle to first temporary binary file.
- * \param tt2: `fstream &` Handle to second temporary binary file.
- */
-void frfmer(
- int nkv, double vkm, double *vkv, double vknmx, double apfafa, double tra,
- double spd, double rir, double ftcn, int le, int lmode, double pmf,
- std::fstream &tt1, std::fstream &tt2
-) {
- const int nlemt = le * (le + 2) * 2;
- const std::complex<double> cc0(0.0, 0.0);
- std::complex<double> *wk = new std::complex<double>[nlemt]();
- double sq = vkm * vkm;
- for (int jy90 = 0; jy90 < nkv; jy90++) {
- double vky = vkv[jy90];
- double sqy = vky * vky;
- for (int jx80 = 0; jx80 < nkv; jx80++) {
- double vkx = vkv[jx80];
- double sqx = vkx * vkx;
- double sqn = sqx + sqy;
- double vkn = sqrt(sqn);
- if (vkn <= vknmx) {
- double vkz = sqrt(sq - sqn);
- wamff(wk, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf);
- for (int j = 0; j < nlemt; j++) {
- double vreal = wk[j].real();
- double vimag = wk[j].imag();
- tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- }
- tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
- } else { // label 50
- for (int j = 0; j < nlemt; j++) {
- double vreal = 0.0;
- double vimag = 0.0;
- tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
- tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
- }
- double vkz = 0.0;
- tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
- }
- } // jx80 loop
- } // jy90 loop
- delete[] wk;
-}
-
/*! C++ porting of CAMP
*
* \param ac: Vector of complex. QUESTION: definition?
@@ -320,13 +66,7 @@ void frfmer(
void camp(
std::complex<double> *ac, std::complex<double> **am0m, std::complex<double> *ws,
CIL *cil
-) {
- for (int j = 0; j < cil->nlemt; j++) {
- for (int i = 0; i < cil->nlemt; i++) {
- ac[j] += (am0m[j][i] * ws[i]);
- } // i loop
- } // j loop
-}
+);
/*! C++ porting of CZAMP
*
@@ -341,31 +81,7 @@ void camp(
void czamp(
std::complex<double> *ac, std::complex<double> **amd, int **indam,
std::complex<double> *ws, CIL *cil
-) {
- const std::complex<double> cc0(0.0, 0.0);
- const std::complex<double> uim(0.0, 1.0);
- std::complex<double> summ, sume;
- for (int im20 = 1; im20 <= cil->mxim; im20++) {
- int m = im20 - cil->mxmpo;
- int abs_m = (m < 0) ? -m : m;
- int lmn = (abs_m > 1) ? abs_m : 1;
- for (int l20 = lmn; l20 <= cil->le; l20++) {
- int i = m + l20 * (l20 + 1);
- int ie = i + cil->nlem;
- summ = cc0;
- sume = cc0;
- for (int ls15 = lmn; ls15 <= cil->le; ls15++) {
- int is = m + ls15 * (ls15 + 15) - 1;
- int ise = is + cil->nlem;
- int num = indam[l20 - 1][ls15 - 1] + m - 1;
- summ += (amd[num][0] * ws[is] + amd[num][1] * ws[ise]);
- sume += (amd[num][2] * ws[is] + amd[num][3] * ws[ise]);
- } // ls15 loop
- ac[i - 1] = summ;
- ac[ie - 1] = sume;
- } // l20 loop
- } // im20 loop
-}
+);
/*! C++ porting of FFRF
*
@@ -380,120 +96,7 @@ void czamp(
void ffrf(
double ****zpv, std::complex<double> *ac, std::complex<double> *ws, double *fffe,
double *fffs, CIL *cil, CCR *ccr
-) {
- const std::complex<double> cc0(0.0, 0.0);
- const std::complex<double> uim(0.0, 1.0);
- //const std::complex<double> sq2iti = uim / sqrt(2.0);
- std::complex<double> uimmp, summ, sume, suem, suee;
- std::complex<double> *gap = new std::complex<double>[3]();
-
- for (int imu50 = 1; imu50 <= 3; imu50++) {
- int mu = imu50 - 2;
- gap[imu50 - 1] = cc0;
- for (int l40 = 1; l40 <= cil->le; l40++) {
- int lpo = l40 + 1;
- int ltpo = lpo + l40;
- int imm = l40 * lpo;
- for (int ilmp40 = 1; ilmp40 <= 3; ilmp40++) {
- if ((l40 == 1 && ilmp40 == 1) || (l40 == cil->le && ilmp40 == 3)) continue; // ilmp40 loop
- int lmpml = ilmp40 - 2;
- int lmp = l40 + lmpml;
- uimmp = uim * (-1.0 * lmpml);
- int impmmmp = lmp * (lmp + 1);
- for (int im30 = 1; im30 <= ltpo; im30++) {
- int m = im30 - lpo;
- int mmp = m - mu;
- int abs_mmp = (mmp < 0) ? -mmp : mmp;
- if (abs_mmp <= lmp) {
- int i = imm + m;
- int ie = i + cil->nlem;
- int imp = impmmmp + mmp;
- int impe = imp + cil->nlem;
- double cgc = cg1(lmpml, mu, l40, m);
- summ = dconjg(ws[i - 1]) * ac[imp - 1];
- sume = dconjg(ws[i - 1]) * ac[impe - 1];
- suem = dconjg(ws[ie - 1]) * ac[imp - 1];
- suee = dconjg(ws[ie - 1]) * ac[impe - 1];
- if (lmpml != 0) {
- summ *= uimmp;
- sume *= uimmp;
- suem *= uimmp;
- suee *= uimmp;
- }
- // label 25
- gap[imu50 - 1] += (cgc * (
- summ * zpv[l40 - 1][ilmp40 - 1][0][0] +
- sume * zpv[l40 - 1][ilmp40 - 1][0][1] +
- suem * zpv[l40 - 1][ilmp40 - 1][1][0] +
- suee * zpv[l40 - 1][ilmp40 - 1][1][1]
- )
- );
- }
- } // im30 loop
- } // ilmp40
- } // l40 loop
- } // imu50 loop
- sume = -gap[0] * ccr->cimu;
- suee = -gap[1] * ccr->cof;
- suem = -gap[2] * ccr->cimu;
- fffe[0] = (sume - suem).real();
- fffe[1] = ((sume + suem) * uim).real();
- fffe[2] = suee.real();
-
- for (int imu90 = 1; imu90 <= 3; imu90++) {
- int mu = imu90 - 2;
- gap[imu90 - 1] = cc0;
- for (int l80 = 1; l80 <= cil->le; l80++) {
- int lpo = l80 + 1;
- int ltpo = lpo + l80;
- int imm = l80 * lpo;
- for (int ilmp80 = 1; ilmp80 <= 3; ilmp80++) {
- if ((l80 == 1 && ilmp80 == 1) || (l80 == cil->le && ilmp80 == 3)) continue; // ilmp80 loop
- int lmpml = ilmp80 - 2;
- int lmp = l80 + lmpml;
- uimmp = uim * (-1.0 * lmpml);
- int impmmmp = lmp * (lmp + 1);
- for (int im70 = 1; im70 <= ltpo; im70++) {
- int m = im70 - lpo;
- int mmp = m - mu;
- int abs_mmp = (mmp < 0) ? -mmp : mmp;
- if (abs_mmp <= lmp) {
- int i = imm + m;
- int ie = i + cil->nlem;
- int imp = impmmmp + mmp;
- int impe = imp + cil->nlem;
- double cgc = cg1(lmpml, mu, l80, m);
- summ = dconjg(ac[i - 1]) * ac[imp - 1];
- sume = dconjg(ac[i - 1]) * ac[impe - 1];
- suem = dconjg(ac[ie - 1]) * ac[imp - 1];
- suee = dconjg(ac[ie - 1]) * ac[impe - 1];
- if (lmpml != 0) {
- summ *= uimmp;
- sume *= uimmp;
- suem *= uimmp;
- suee *= uimmp;
- }
- // label 65
- gap[imu90 - 1] += (cgc * (
- summ * zpv[l80 - 1][ilmp80 - 1][0][0] +
- sume * zpv[l80 - 1][ilmp80 - 1][0][1] +
- suem * zpv[l80 - 1][ilmp80 - 1][1][0] +
- suee * zpv[l80 - 1][ilmp80 - 1][1][1]
- )
- );
- }
- } // im70 loop
- } // ilmp80 loop
- } // l80 loop
- } // imu90 loop
- sume = gap[0] * ccr->cimu;
- suee = gap[1] * ccr->cof;
- suem = gap[2] * ccr->cimu;
- fffs[0] = (sume - suem).real();
- fffs[1] = ((sume + suem) * uim).real();
- fffs[2] = suee.real();
- delete[] gap;
-}
+);
/*! C++ porting of FFRT
*
@@ -507,65 +110,40 @@ void ffrf(
void ffrt(
std::complex<double> *ac, std::complex<double> *ws, double *ffte, double *ffts,
CIL *cil
-) {
- const std::complex<double> cc0(0.0, 0.0);
- const std::complex<double> uim(0.0, 1.0);
- const double sq2i = 1.0 / sqrt(2.0);
- const std::complex<double> sq2iti = uim * sq2i;
- std::complex<double> aca, acw;
- std::complex<double> *ctqce, *ctqcs;
+);
- ctqce = new std::complex<double>[3]();
- ctqcs = new std::complex<double>[3]();
- for (int l60 = 1; l60 < cil->le; l60++) {
- int lpo = l60 + 1;
- int il = l60 * lpo;
- int ltpo = l60 + lpo;
- for (int im60 = 1; im60 <= ltpo; im60++) {
- double rmu;
- int m = im60 - lpo;
- int i = m + il;
- int ie = i + cil->nlem;
- int mmmu = m + 1;
- int mmmmu = (mmmu < 0) ? -mmmu: mmmu;
- if (mmmmu <= l60) {
- int immu = mmmu + il;
- int immue = immu + cil->nlem;
- rmu = -sqrt(1.0 * (l60 + mmmu) * (l60 - m)) * sq2i;
- acw = dconjg(ac[i - 1]) * ws[immu - 1] + dconjg(ac[ie - 1]) * ws[immue - 1];
- aca = dconjg(ac[i - 1]) * ac[immu - 1] + dconjg(ac[ie - 1]) * ac[immue - 1];
- ctqce[0] += (acw * rmu);
- ctqcs[0] += (aca * rmu);
- }
- // label 30
- rmu = -1.0 * m;
- acw = dconjg(ac[i - 1]) * ws[i - 1] + dconjg(ac[ie - 1]) * ws[ie - 1];
- aca = dconjg(ac[i - 1]) * ac[i - 1] + dconjg(ac[ie - 1]) * ac[ie - 1];
- ctqce[1] += (acw * rmu);
- ctqcs[1] += (aca * rmu);
- mmmu = m - 1;
- mmmmu = (mmmu < 0) ? -mmmu: mmmu;
- if (mmmmu <= l60) {
- int immu = mmmu + il;
- int immue = immu + cil->nlem;
- rmu = sqrt(1.0 * (l60 - mmmu) * (l60 + m)) * sq2i;
- acw = dconjg(ac[i - 1]) * ws[immu - 1] + dconjg(ac[ie - 1]) * ws[immue - 1];
- aca = dconjg(ac[i - 1]) * ac[immu - 1] + dconjg(ac[ie - 1]) * ac[immue - 1];
- ctqce[2] += (acw * rmu);
- ctqcs[2] += (aca * rmu);
- }
- } // im60 loop
- } // l60 loop
- ffte[0] = (ctqce[0] - ctqce[2]).real() * sq2i;
- ffte[1] = (sq2iti * (ctqce[0] + ctqce[2])).real();
- ffte[2] = ctqce[1].real();
- ffts[0] = -sq2i * (ctqcs[0] - ctqcs[2]).real();
- ffts[1] = -1.0 * (sq2iti * (ctqcs[0] + ctqcs[2])).real();
- ffts[2] = -1.0 * ctqcs[1].real();
+/*! C++ porting of FRFMER
+ *
+ * \param nkv: `int` QUESTION: definition?
+ * \param vkm: `double` QUESTION: definition?
+ * \param vkv: `double *` QUESTION: definition?
+ * \param vknmx: `double` QUESTION: definition?
+ * \param apfafa: `double` QUESTION: definition?
+ * \param tra: `double` QUESTION: definition?
+ * \param spd: `double` QUESTION: definition?
+ * \param rir: `double` QUESTION: definition?
+ * \param ftcn: `double` QUESTION: definition?
+ * \param le: `int` QUESTION: definition?
+ * \param lmode: `int` QUESTION: definition?
+ * \param pmf: `double` QUESTION: definition?
+ * \param tt1: `fstream &` Handle to first temporary binary file.
+ * \param tt2: `fstream &` Handle to second temporary binary file.
+ */
+void frfmer(
+ int nkv, double vkm, double *vkv, double vknmx, double apfafa, double tra,
+ double spd, double rir, double ftcn, int le, int lmode, double pmf,
+ std::fstream &tt1, std::fstream &tt2
+);
- delete[] ctqce;
- delete[] ctqcs;
-}
+/*! C++ porting of PWMALP
+ *
+ * \param w: Matrix of complex. QUESTION: definition?
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param ylm: Vector of complex
+ * \param lw: `int`
+ */
+void pwmalp(std::complex<double> **w, double *up, double *un, std::complex<double> *ylm, int lw);
/*! C++ porting of SAMP
*
@@ -580,19 +158,7 @@ void ffrt(
void samp(
std::complex<double> *ac, std::complex<double> *tmsm, std::complex<double> *tmse,
std::complex<double> *ws, CIL *cil
-) {
- int i = 0;
- for (int l20 = 0; l20 < cil->le; l20++) {
- int l = l20 + 1;
- int ltpo = l + l + 1;
- for (int im20 = 0; im20 < ltpo; im20++) {
- int ie = i + cil->nlem;
- ac[i] = tmsm[l20] * ws[i];
- ac[ie] = tmse[l20] * ws[ie];
- i++;
- } // im20 loop
- } // l20 loop
-}
+);
/*! C++ porting of SAMPOA
*
@@ -606,26 +172,24 @@ void samp(
void sampoa(
std::complex<double> *ac, std::complex<double> **tms, std::complex<double> *ws,
CIL *cil
-) {
- std::complex<double> **tm = new std::complex<double>*[2];
- tm[0] = new std::complex<double>[2]();
- tm[1] = new std::complex<double>[2]();
- int i = 0;
- for (int l20 = 0; l20 < cil->le; l20++) {
- tm[0][0] = tms[l20][0];
- tm[0][1] = tms[l20][1];
- tm[1][1] = tms[l20][2];
- tm[1][0] = tm[0][1];
- int l = l20 + 1;
- int ltpo = l + l + 1;
- for (int im20 = 0; im20 < ltpo; im20++) {
- int ie = i + cil->nlem;
- ac[i] = tm[0][0] * ws[i] + tm[0][1] * ws[ie];
- ac[ie] = tm[1][0] * ws[i] + tm[1][1] * ws[ie];
- i++;
- } // im20 loop
- } // l20 loop
- delete[] tm[1];
- delete[] tm[0];
- delete[] tm;
-}
+);
+
+/*! C++ porting of WAMFF
+ *
+ * \param wk: Vector of complex. QUESTION: definition?
+ * \param x: `double`
+ * \param y: `double`
+ * \param z: `double`
+ * \param lm: `int`
+ * \param apfafa: `double` QUESTION: definition?
+ * \param tra: `double` QUESTION: definition?
+ * \param spd: `double` QUESTION: definition?
+ * \param rir: `double` QUESTION: definition?
+ * \param ftcn: `double` QUESTION: definition?
+ * \param lmode: `int` QUESTION: definition?
+ * \param pmf: `double` QUESTION: definition?
+ */
+void wamff(
+ std::complex<double> *wk, double x, double y, double z, int lm, double apfafa,
+ double tra, double spd, double rir, double ftcn, int lmode, double pmf
+);
diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..73d3809866fa8541e1f2720d95b1bdabf8891a6f
--- /dev/null
+++ b/src/libnptm/tra_subs.cpp
@@ -0,0 +1,494 @@
+/*! \file tra_subs.cpp
+ *
+ * \brief C++ implementation of TRAPPING subroutines.
+ */
+
+#include "../include/tra_subs.h"
+
+using namespace std;
+
+void camp(
+ complex<double> *ac, complex<double> **am0m, complex<double> *ws,
+ CIL *cil
+) {
+ for (int j = 0; j < cil->nlemt; j++) {
+ for (int i = 0; i < cil->nlemt; i++) {
+ ac[j] += (am0m[j][i] * ws[i]);
+ } // i loop
+ } // j loop
+}
+
+void czamp(
+ complex<double> *ac, complex<double> **amd, int **indam,
+ complex<double> *ws, CIL *cil
+) {
+ const complex<double> cc0(0.0, 0.0);
+ const complex<double> uim(0.0, 1.0);
+ complex<double> summ, sume;
+ for (int im20 = 1; im20 <= cil->mxim; im20++) {
+ int m = im20 - cil->mxmpo;
+ int abs_m = (m < 0) ? -m : m;
+ int lmn = (abs_m > 1) ? abs_m : 1;
+ for (int l20 = lmn; l20 <= cil->le; l20++) {
+ int i = m + l20 * (l20 + 1);
+ int ie = i + cil->nlem;
+ summ = cc0;
+ sume = cc0;
+ for (int ls15 = lmn; ls15 <= cil->le; ls15++) {
+ int is = m + ls15 * (ls15 + 15) - 1;
+ int ise = is + cil->nlem;
+ int num = indam[l20 - 1][ls15 - 1] + m - 1;
+ summ += (amd[num][0] * ws[is] + amd[num][1] * ws[ise]);
+ sume += (amd[num][2] * ws[is] + amd[num][3] * ws[ise]);
+ } // ls15 loop
+ ac[i - 1] = summ;
+ ac[ie - 1] = sume;
+ } // l20 loop
+ } // im20 loop
+}
+
+void ffrf(
+ double ****zpv, complex<double> *ac, complex<double> *ws, double *fffe,
+ double *fffs, CIL *cil, CCR *ccr
+) {
+ const complex<double> cc0(0.0, 0.0);
+ const complex<double> uim(0.0, 1.0);
+ complex<double> uimmp, summ, sume, suem, suee;
+ complex<double> *gap = new complex<double>[3]();
+
+ for (int imu50 = 1; imu50 <= 3; imu50++) {
+ int mu = imu50 - 2;
+ gap[imu50 - 1] = cc0;
+ for (int l40 = 1; l40 <= cil->le; l40++) {
+ int lpo = l40 + 1;
+ int ltpo = lpo + l40;
+ int imm = l40 * lpo;
+ for (int ilmp40 = 1; ilmp40 <= 3; ilmp40++) {
+ if ((l40 == 1 && ilmp40 == 1) || (l40 == cil->le && ilmp40 == 3)) continue; // ilmp40 loop
+ int lmpml = ilmp40 - 2;
+ int lmp = l40 + lmpml;
+ uimmp = uim * (-1.0 * lmpml);
+ int impmmmp = lmp * (lmp + 1);
+ for (int im30 = 1; im30 <= ltpo; im30++) {
+ int m = im30 - lpo;
+ int mmp = m - mu;
+ int abs_mmp = (mmp < 0) ? -mmp : mmp;
+ if (abs_mmp <= lmp) {
+ int i = imm + m;
+ int ie = i + cil->nlem;
+ int imp = impmmmp + mmp;
+ int impe = imp + cil->nlem;
+ double cgc = cg1(lmpml, mu, l40, m);
+ summ = dconjg(ws[i - 1]) * ac[imp - 1];
+ sume = dconjg(ws[i - 1]) * ac[impe - 1];
+ suem = dconjg(ws[ie - 1]) * ac[imp - 1];
+ suee = dconjg(ws[ie - 1]) * ac[impe - 1];
+ if (lmpml != 0) {
+ summ *= uimmp;
+ sume *= uimmp;
+ suem *= uimmp;
+ suee *= uimmp;
+ }
+ // label 25
+ gap[imu50 - 1] += (cgc * (
+ summ * zpv[l40 - 1][ilmp40 - 1][0][0] +
+ sume * zpv[l40 - 1][ilmp40 - 1][0][1] +
+ suem * zpv[l40 - 1][ilmp40 - 1][1][0] +
+ suee * zpv[l40 - 1][ilmp40 - 1][1][1]
+ )
+ );
+ }
+ } // im30 loop
+ } // ilmp40
+ } // l40 loop
+ } // imu50 loop
+ sume = -gap[0] * ccr->cimu;
+ suee = -gap[1] * ccr->cof;
+ suem = -gap[2] * ccr->cimu;
+ fffe[0] = (sume - suem).real();
+ fffe[1] = ((sume + suem) * uim).real();
+ fffe[2] = suee.real();
+
+ for (int imu90 = 1; imu90 <= 3; imu90++) {
+ int mu = imu90 - 2;
+ gap[imu90 - 1] = cc0;
+ for (int l80 = 1; l80 <= cil->le; l80++) {
+ int lpo = l80 + 1;
+ int ltpo = lpo + l80;
+ int imm = l80 * lpo;
+ for (int ilmp80 = 1; ilmp80 <= 3; ilmp80++) {
+ if ((l80 == 1 && ilmp80 == 1) || (l80 == cil->le && ilmp80 == 3)) continue; // ilmp80 loop
+ int lmpml = ilmp80 - 2;
+ int lmp = l80 + lmpml;
+ uimmp = uim * (-1.0 * lmpml);
+ int impmmmp = lmp * (lmp + 1);
+ for (int im70 = 1; im70 <= ltpo; im70++) {
+ int m = im70 - lpo;
+ int mmp = m - mu;
+ int abs_mmp = (mmp < 0) ? -mmp : mmp;
+ if (abs_mmp <= lmp) {
+ int i = imm + m;
+ int ie = i + cil->nlem;
+ int imp = impmmmp + mmp;
+ int impe = imp + cil->nlem;
+ double cgc = cg1(lmpml, mu, l80, m);
+ summ = dconjg(ac[i - 1]) * ac[imp - 1];
+ sume = dconjg(ac[i - 1]) * ac[impe - 1];
+ suem = dconjg(ac[ie - 1]) * ac[imp - 1];
+ suee = dconjg(ac[ie - 1]) * ac[impe - 1];
+ if (lmpml != 0) {
+ summ *= uimmp;
+ sume *= uimmp;
+ suem *= uimmp;
+ suee *= uimmp;
+ }
+ // label 65
+ gap[imu90 - 1] += (cgc * (
+ summ * zpv[l80 - 1][ilmp80 - 1][0][0] +
+ sume * zpv[l80 - 1][ilmp80 - 1][0][1] +
+ suem * zpv[l80 - 1][ilmp80 - 1][1][0] +
+ suee * zpv[l80 - 1][ilmp80 - 1][1][1]
+ )
+ );
+ }
+ } // im70 loop
+ } // ilmp80 loop
+ } // l80 loop
+ } // imu90 loop
+ sume = gap[0] * ccr->cimu;
+ suee = gap[1] * ccr->cof;
+ suem = gap[2] * ccr->cimu;
+ fffs[0] = (sume - suem).real();
+ fffs[1] = ((sume + suem) * uim).real();
+ fffs[2] = suee.real();
+ delete[] gap;
+}
+
+void ffrt(
+ complex<double> *ac, complex<double> *ws, double *ffte, double *ffts,
+ CIL *cil
+) {
+ const complex<double> cc0(0.0, 0.0);
+ const complex<double> uim(0.0, 1.0);
+ const double sq2i = 1.0 / sqrt(2.0);
+ const complex<double> sq2iti = uim * sq2i;
+ complex<double> aca, acw;
+ complex<double> *ctqce, *ctqcs;
+
+ ctqce = new complex<double>[3]();
+ ctqcs = new complex<double>[3]();
+ for (int l60 = 1; l60 < cil->le; l60++) {
+ int lpo = l60 + 1;
+ int il = l60 * lpo;
+ int ltpo = l60 + lpo;
+ for (int im60 = 1; im60 <= ltpo; im60++) {
+ double rmu;
+ int m = im60 - lpo;
+ int i = m + il;
+ int ie = i + cil->nlem;
+ int mmmu = m + 1;
+ int mmmmu = (mmmu < 0) ? -mmmu: mmmu;
+ if (mmmmu <= l60) {
+ int immu = mmmu + il;
+ int immue = immu + cil->nlem;
+ rmu = -sqrt(1.0 * (l60 + mmmu) * (l60 - m)) * sq2i;
+ acw = dconjg(ac[i - 1]) * ws[immu - 1] + dconjg(ac[ie - 1]) * ws[immue - 1];
+ aca = dconjg(ac[i - 1]) * ac[immu - 1] + dconjg(ac[ie - 1]) * ac[immue - 1];
+ ctqce[0] += (acw * rmu);
+ ctqcs[0] += (aca * rmu);
+ }
+ // label 30
+ rmu = -1.0 * m;
+ acw = dconjg(ac[i - 1]) * ws[i - 1] + dconjg(ac[ie - 1]) * ws[ie - 1];
+ aca = dconjg(ac[i - 1]) * ac[i - 1] + dconjg(ac[ie - 1]) * ac[ie - 1];
+ ctqce[1] += (acw * rmu);
+ ctqcs[1] += (aca * rmu);
+ mmmu = m - 1;
+ mmmmu = (mmmu < 0) ? -mmmu: mmmu;
+ if (mmmmu <= l60) {
+ int immu = mmmu + il;
+ int immue = immu + cil->nlem;
+ rmu = sqrt(1.0 * (l60 - mmmu) * (l60 + m)) * sq2i;
+ acw = dconjg(ac[i - 1]) * ws[immu - 1] + dconjg(ac[ie - 1]) * ws[immue - 1];
+ aca = dconjg(ac[i - 1]) * ac[immu - 1] + dconjg(ac[ie - 1]) * ac[immue - 1];
+ ctqce[2] += (acw * rmu);
+ ctqcs[2] += (aca * rmu);
+ }
+ } // im60 loop
+ } // l60 loop
+ ffte[0] = (ctqce[0] - ctqce[2]).real() * sq2i;
+ ffte[1] = (sq2iti * (ctqce[0] + ctqce[2])).real();
+ ffte[2] = ctqce[1].real();
+ ffts[0] = -sq2i * (ctqcs[0] - ctqcs[2]).real();
+ ffts[1] = -1.0 * (sq2iti * (ctqcs[0] + ctqcs[2])).real();
+ ffts[2] = -1.0 * ctqcs[1].real();
+
+ delete[] ctqce;
+ delete[] ctqcs;
+}
+
+void frfmer(
+ int nkv, double vkm, double *vkv, double vknmx, double apfafa, double tra,
+ double spd, double rir, double ftcn, int le, int lmode, double pmf,
+ std::fstream &tt1, std::fstream &tt2
+) {
+ const int nlemt = le * (le + 2) * 2;
+ const complex<double> cc0(0.0, 0.0);
+ complex<double> *wk = new complex<double>[nlemt]();
+ double sq = vkm * vkm;
+ for (int jy90 = 0; jy90 < nkv; jy90++) {
+ double vky = vkv[jy90];
+ double sqy = vky * vky;
+ for (int jx80 = 0; jx80 < nkv; jx80++) {
+ double vkx = vkv[jx80];
+ double sqx = vkx * vkx;
+ double sqn = sqx + sqy;
+ double vkn = sqrt(sqn);
+ if (vkn <= vknmx) {
+ double vkz = sqrt(sq - sqn);
+ wamff(wk, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf);
+ for (int j = 0; j < nlemt; j++) {
+ double vreal = wk[j].real();
+ double vimag = wk[j].imag();
+ tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
+ tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
+ }
+ tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
+ } else { // label 50
+ for (int j = 0; j < nlemt; j++) {
+ double vreal = 0.0;
+ double vimag = 0.0;
+ tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
+ tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
+ }
+ double vkz = 0.0;
+ tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
+ }
+ } // jx80 loop
+ } // jy90 loop
+ delete[] wk;
+}
+
+void pwmalp(complex<double> **w, double *up, double *un, complex<double> *ylm, int lw) {
+ complex<double> cp1, cm1, cp2, cm2, cl;
+ const complex<double> uim(0.0, 1.0);
+ const double four_pi = acos(0.0) * 8.0;
+ const int nlwm = lw * (lw + 2);
+ cm1 = 0.5 * complex<double>(up[0], up[1]);
+ cp1 = 0.5 * complex<double>(up[0], -up[1]);
+ double cz1 = up[2];
+ cm2 = 0.5 * complex<double>(un[0], un[1]);
+ cp2 = 0.5 * 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;
+ complex<double> cl = std::pow((four_pi / sqrt(x)) * uim, 1.0 * l20);
+ int mv = l20 + lf;
+ int m = -lf;
+ for (int mf20 = 1; mf20 <= mv; mf20++) {
+ m++;
+ 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;
+ w[k - 1][0] = dconjg(cp1 * cp * ylm[k + 1] + cm1 * cm * ylm[k - 1] + cz1 * cz * ylm[k]) * cl;
+ w[k - 1][1] = dconjg(cp2 * cp * ylm[k + 1] + cm2 * cm * ylm[k - 1] + cz2 * cz * ylm[k]) * cl;
+ } // mf20 loop
+ } // l20 loop
+ for (int k30 = 0; k30 < nlwm; k30++) {
+ int i = k30 + nlwm;
+ w[i][0] = uim * w[k30][1];
+ w[i][1] = -uim * w[k30][0];
+ } // k30 loop
+}
+
+void samp(
+ complex<double> *ac, complex<double> *tmsm, complex<double> *tmse,
+ complex<double> *ws, CIL *cil
+) {
+ int i = 0;
+ for (int l20 = 0; l20 < cil->le; l20++) {
+ int l = l20 + 1;
+ int ltpo = l + l + 1;
+ for (int im20 = 0; im20 < ltpo; im20++) {
+ int ie = i + cil->nlem;
+ ac[i] = tmsm[l20] * ws[i];
+ ac[ie] = tmse[l20] * ws[ie];
+ i++;
+ } // im20 loop
+ } // l20 loop
+}
+
+void sampoa(
+ complex<double> *ac, complex<double> **tms, complex<double> *ws,
+ CIL *cil
+) {
+ complex<double> **tm = new complex<double>*[2];
+ tm[0] = new complex<double>[2]();
+ tm[1] = new complex<double>[2]();
+ int i = 0;
+ for (int l20 = 0; l20 < cil->le; l20++) {
+ tm[0][0] = tms[l20][0];
+ tm[0][1] = tms[l20][1];
+ tm[1][1] = tms[l20][2];
+ tm[1][0] = tm[0][1];
+ int l = l20 + 1;
+ int ltpo = l + l + 1;
+ for (int im20 = 0; im20 < ltpo; im20++) {
+ int ie = i + cil->nlem;
+ ac[i] = tm[0][0] * ws[i] + tm[0][1] * ws[ie];
+ ac[ie] = tm[1][0] * ws[i] + tm[1][1] * ws[ie];
+ i++;
+ } // im20 loop
+ } // l20 loop
+ delete[] tm[1];
+ delete[] tm[0];
+ delete[] tm;
+}
+
+void wamff(
+ complex<double> *wk, double x, double y, double z, int lm, double apfafa,
+ double tra, double spd, double rir, double ftcn, int lmode, double pmf
+) {
+ const int nlmm = lm * (lm + 2);
+ const int nlmmt = 2 * nlmm;
+ const int nlmp = nlmm + 2;
+ complex<double> **w, *ylm;
+ const complex<double> cc0(0.0, 0.0);
+ const complex<double> uim(0.0, 1.0);
+ complex<double> cfam, cf1, cf2;
+ double rho, cph, sph, cth, sth, r;
+ double ftc1, ftc2;
+ double *up = new double[3];
+ double *un = new double[3];
+ w = new complex<double>*[nlmmt];
+ for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[2]();
+ ylm = new complex<double>[nlmp]();
+ bool onx = (y == 0.0);
+ bool ony = (x == 0.0);
+ bool onz = (onx && ony);
+ if (!(onz && onx && ony)) {
+ rho = sqrt(x * x + y * y);
+ cph = x / rho;
+ sph = y = rho;
+ } else {
+ if (onz) { // label 10
+ cph = 1.0;
+ sph = 0.0;
+ } else {
+ if (onx) { // label 12
+ rho = sqrt(x * x);
+ cph = (x < 0.0)? -1.0 : 1.0;
+ sph = 0.0;
+ } else {
+ if (ony) { // label 13
+ rho = sqrt(y * y);
+ cph = 0.0;
+ sph = (y < 0.0)? -1.0: 1.0;
+ }
+ }
+ }
+ }
+ // label 15
+ if (z == 0.0) {
+ if (!onz) {
+ r = rho;
+ cth = 0.0;
+ sth = 1.0;
+ } else { // label 17
+ r = 0.0;
+ cth = 1.0;
+ sth = 0.0;
+ }
+ } else { // label 18
+ if (!onz) {
+ r = sqrt(rho * rho + z * z);
+ cth = z / r;
+ sth = rho / r;
+ } else { // label 20
+ r = sqrt(z * z);
+ cth = (z < 0.0)? -1.0: 1.0;
+ sth = 0.0;
+ }
+ }
+ if (lmode == 0 || sth != 0.0) { // label 25
+ ylm[nlmp - 1] = cc0;
+ sphar(cth, sth, cph, sph, lm, ylm);
+ up[0] = cph * cth;
+ up[1] = sph * cth;
+ up[2] = -sth;
+ un[0] = -sph;
+ un[1] = cph;
+ un[2] = 0.0;
+ // Would call PWMALP(W,UP,UN,YLM,LM)
+ pwmalp(w, up, un, ylm, lm);
+ double apfa = sth * apfafa;
+ if (spd > 0.0) {
+ double sthl = sth * rir;
+ double cthl = sqrt(1.0 - sthl * sthl);
+ double arg = spd * (z - (r / rir) * cthl);
+ cfam = (tra * std::exp(-apfa * apfa) / sqrt(cthl)) * std::exp(uim * arg);
+ if (lmode == 0) { // CHECK: Computed GO TO, not sure what happens with LMODE = 0
+ if (sth == 0.0) { // label 45
+ ftc1 = ftcn;
+ ftc2 = ftcn;
+ // goes to 48
+ }
+ } else if (lmode == 1) { // label 46
+ cfam *= ((cph + uim * sph) * sth * pmf);
+ ftc1 = 2.0 * cthl / (cthl * rir + cth);
+ ftc2 = 2.0 * cthl / (cthl + cth * rir);
+ // follows on to 48
+ } else if (lmode == 2) { // label 50
+ ftc1 = 2.0 * cthl / (cthl * rir + cth);
+ cfam *= (sth * pmf * ftc1);
+ for (int i52 = 0; i52 < nlmmt; i52++) wk[i52] = w[i52][0] * cfam;
+ // returns
+ } else if (lmode == 3) { // label 53
+ ftc2 = 2.0 * cthl / (cthl + cth * rir);
+ cfam *= (sth * pmf * ftc2);
+ for (int i55 = 0; i55 < nlmmt; i55++) wk[i55] = w[i55][1] * cfam;
+ // returns
+ }
+ if (lmode == 0 || lmode == 1) { //label 48
+ cf1 = cph * ftc1 * cfam;
+ cf2 = -sph * ftc2 * cfam;
+ // goes to 62
+ }
+ } else { // label 57
+ double fam = tra * std::exp(-apfa * apfa) / sqrt(cth);
+ if (lmode == 0) {
+ double f1 = cph * fam;
+ double f2 = -sph * fam;
+ for (int i58 = 0; i58 < nlmmt; i58++) wk[i58] = w[i58][0] * f1 + w[i58][1] * f2;
+ // returns
+ } else if (lmode == 1) { // label 60
+ cfam = (pmf * sth * fam) * (cph * uim * sph);
+ cf1 = cph * cfam;
+ cf2 = -sph * cfam;
+ // follows on to 62
+ } else if (lmode == 2) { // label 65
+ fam *= (pmf * sth);
+ for (int i67 = 0; i67 < nlmmt; i67++) wk[i67] = w[i67][0] * fam;
+ // returns
+ } else if (lmode == 3) { // label 68
+ fam *= (pmf * sth);
+ for (int i70 = 0; i70 < nlmmt; i70++) wk[i70] = w[i70][1] * fam;
+ // returns
+ }
+ }
+ if (lmode == 0 || lmode == 1) { // label 62
+ for (int i63 = 0; i63 < nlmmt; i63++) wk[i63] = w[i63][0] * cf1 + w[i63][1] * cf2;
+ }
+ }
+ // Clean up memory
+ delete[] up;
+ delete[] un;
+ for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
+ delete[] w;
+ delete[] ylm;
+}
diff --git a/src/trapping/Makefile b/src/trapping/Makefile
index bfad43bddcd843a6d71010b2a39c6e36d09e403d..635ccaa9b37ecae2b527c3539334658135d92d30 100644
--- a/src/trapping/Makefile
+++ b/src/trapping/Makefile
@@ -2,8 +2,11 @@ BUILDDIR=../../build/trapping
FC=gfortran
FCFLAGS=-std=legacy -O3
LFLAGS=
+CXX=g++
+CXXFLAGS=-O2 -ggdb -pg -coverage
+CXXLFLAGS=
-all: frfme lffft
+all: frfme lffft np_frfme np_lffft
frfme: frfme.o
$(FC) $(FCFLAGS) -o $(BUILDDIR)/frfme $(BUILDDIR)/frfme.o $(LFLAGS)
@@ -11,6 +14,33 @@ frfme: frfme.o
lffft: lffft.o
$(FC) $(FCFLAGS) -o $(BUILDDIR)/lffft $(BUILDDIR)/lffft.o $(LFLAGS)
+np_frfme: $(BUILDDIR)/cfrfme.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/tra_subs.o
+ $(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_frfme $(BUILDDIR)/cfrfme.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/tra_subs.o
+
+np_lffft: $(BUILDDIR)/clffft.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/tra_subs.o
+ $(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_lffft $(BUILDDIR)/clffft.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/tra_subs.o
+
+$(BUILDDIR)/cfrfme.o:
+ $(CXX) $(CXXFLAGS) frfme.cpp -c -o $(BUILDDIR)/cfrfme.o
+
+$(BUILDDIR)/clffft.o:
+ $(CXX) $(CXXFLAGS) lffft.cpp -c -o $(BUILDDIR)/clffft.o
+
+$(BUILDDIR)/Commons.o:
+ $(CXX) $(CXXFLAGS) ../libnptm/Commons.cpp -c -o $(BUILDDIR)/Commons.o
+
+$(BUILDDIR)/Configuration.o:
+ $(CXX) $(CXXFLAGS) ../libnptm/Configuration.cpp -c -o $(BUILDDIR)/Configuration.o
+
+$(BUILDDIR)/Parsers.o:
+ $(CXX) $(CXXFLAGS) ../libnptm/Parsers.cpp -c -o $(BUILDDIR)/Parsers.o
+
+$(BUILDDIR)/sph_subs.o:
+ $(CXX) $(CXXFLAGS) ../libnptm/sph_subs.cpp -c -o $(BUILDDIR)/sph_subs.o
+
+ $(BUILDDIR)/tra_subs.o:
+ $(CXX) $(CXXFLAGS) ../libnptm/tra_subs.cpp -c -o $(BUILDDIR)/tra_subs.o
+
clean:
rm -f $(BUILDDIR)/*.o
diff --git a/src/trapping/frfme.cpp b/src/trapping/frfme.cpp
index 9adbda7bf64ef1ae787130e779e9a16638c2498f..c0084283158bab19fdbbf6f303098828d23b3556 100644
--- a/src/trapping/frfme.cpp
+++ b/src/trapping/frfme.cpp
@@ -1,8 +1,11 @@
+/*! \file frfme.cpp
+ *
+ * \brief C++ implementation of FRFME.
+ */
+
#include <cstdio>
-#include <fstream>
#include <regex>
#include <string>
-#include <complex>
#ifndef INCLUDE_PARSERS_H_
#include "../include/Parsers.h"
#endif