Port WAMFF subroutine to C++

#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?

 * \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?