diff --git a/src/include/tra_subs.h b/src/include/tra_subs.h
index f4f50ecd690474032a7e2e5795a9ef5349f9cfb8..6e1ebe05a8087bdedd7f56e2beca690508570333 100644
--- a/src/include/tra_subs.h
+++ b/src/include/tra_subs.h
@@ -1,10 +1,175 @@
#include <fstream>
+#include <cmath>
#include <complex>
#ifndef INCLUDE_TRA_SUBS_H_
#define INCLUDE_TRA_SUBS_H_
#endif
+//Externally defined subroutines
+extern void sphar(
+ double cosrth, double sinrth, double cosrph, double sinrph,
+ int ll, std::complex<double> *ylm
+);
+//End of externally defined subroutines
+
+/*! 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)
+ 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?
@@ -15,7 +180,7 @@
* \param tra: `double` QUESTION: definition?
* \param spd: `double` QUESTION: definition?
* \param rir: `double` QUESTION: definition?
- * \param fctn: `double` QUESTION: definition?
+ * \param ftcn: `double` QUESTION: definition?
* \param le: `int` QUESTION: definition?
* \param lmode: `int` QUESTION: definition?
* \param pmf: `double` QUESTION: definition?