/*! \file clu_subs.cpp
 *
 * \brief C++ implementation of CLUSTER subroutines.
 */
#include <complex>

#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
#endif

#ifndef INCLUDE_SPH_SUBS_H_
#include "../include/sph_subs.h"
#endif

#ifndef INCLUDE_CLU_SUBS_H_
#include "../include/clu_subs.h"
#endif

using namespace std;

void apc(
	 double ****zpv, int le, complex<double> **am0m, complex<double> **w,
	 double sqk, double **gapr, complex<double> **gapp
) {
  complex<double> **ac, **gap;
  const complex<double> cc0(0.0, 0.0);
  const complex<double> uim(0.0, 1.0);
  complex<double> uimmp, summ, sume, suem, suee, summp, sumep;
  complex<double> suemp, sueep;
  double cof = 1.0 / sqk;
  double cimu = cof / sqrt(2.0);
  int nlem = le * (le + 2);
  const int nlemt = nlem + nlem;
  ac = new complex<double>*[nlemt];
  gap = new complex<double>*[3];
  for (int ai = 0; ai < nlemt; ai++) ac[ai] = new complex<double>[2]();
  for (int gi = 0; gi < 3; gi++) gap[gi] = new complex<double>[2]();
  for (int j45 = 1; j45 <= nlemt; j45++) {
    int j = j45 - 1;
    ac[j][0] = cc0;
    ac[j][1] = cc0;
    for (int i45 = 1; i45 <= nlemt; i45++) {
      int i = i45 - 1;
      ac[j][0] += (am0m[j][i] * w[i][0]);
      ac[j][1] += (am0m[j][i] * w[i][1]);
    } //i45 loop
  } //j45 loop
  for (int imu90 = 1; imu90 <=3; imu90++) {
    int mu = imu90 - 2;
    gap[imu90 - 1][0] = cc0;
    gap[imu90 - 1][1] = cc0;
    gapp[imu90 - 1][0] = cc0;
    gapp[imu90 - 1][1] = cc0;
    for (int l80 =1; l80 <= le; l80++) {
      int lpo = l80 + 1;
      int ltpo = lpo + l80;
      int imm = l80 * lpo;
      for (int ilmp = 1; ilmp <= 3; ilmp++) {
	if ((l80 == 1 && ilmp == 1) || (l80 == le && ilmp == 3)) continue; // ilmp loop
	int lmpml = ilmp - 2;
	int lmp = l80 + lmpml;
	uimmp = (-1.0 * lmpml) * uim;
	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 + nlem;
	    int imp = impmmmp + mmp;
	    int impe = imp + nlem;
	    double cgc = cg1(lmpml, mu, l80, m);
	    int jpo = 2;
	    for (int ipo = 1; ipo <= 2; ipo++) {
	      if (ipo == 2) jpo = 1;
	      summ = dconjg(ac[i - 1][ipo - 1]) * ac[imp - 1][ipo - 1];
	      sume = dconjg(ac[i - 1][ipo - 1]) * ac[impe - 1][ipo - 1];
	      suem = dconjg(ac[ie - 1][ipo - 1]) * ac[imp - 1][ipo - 1];
	      suee = dconjg(ac[ie - 1][ipo - 1]) * ac[impe - 1][ipo - 1];
	      summp = dconjg(ac[i - 1][jpo - 1]) * ac[imp - 1][ipo - 1];
	      sumep = dconjg(ac[i - 1][jpo - 1]) * ac[impe - 1][ipo - 1];
	      suemp = dconjg(ac[ie - 1][jpo - 1]) * ac[imp - 1][ipo - 1];
	      sueep = dconjg(ac[ie - 1][jpo - 1]) * ac[impe - 1][ipo - 1];
	      if (lmpml != 0) {
		summ *= uimmp;
		sume *= uimmp;
		suem *= uimmp;
		suee *= uimmp;
		summp *= uimmp;
		sumep *= uimmp;
		suemp *= uimmp;
		sueep *= uimmp;
	      }
	      // label 55
	      gap[imu90 - 1][ipo - 1] += (
					  (
					   summ * zpv[l80 - 1][ilmp - 1][0][0]
					   + sume * zpv[l80 - 1][ilmp - 1][0][1]
					   + suem * zpv[l80 - 1][ilmp - 1][1][0]
					   + suee * zpv[l80 - 1][ilmp - 1][1][1]
					   ) * cgc
					  );
	      gapp[imu90 - 1][ipo - 1] += (
					   (
					    summp * zpv[l80 - 1][ilmp - 1][0][0]
					    + sumep * zpv[l80 - 1][ilmp - 1][0][1]
					    + suemp * zpv[l80 - 1][ilmp - 1][1][0]
					    + sueep * zpv[l80 - 1][ilmp - 1][1][1]
					    ) * cgc
					   );
	    } // ipo loop
	  } // ends im70 loop
	} // im70 loop
      } // ilmp loop
    } // l80 loop
  } // imu90 loop
  for (int ipo95 = 1; ipo95 <= 2; ipo95++) {
    sume = gap[0][ipo95 - 1] * cimu;
    suee = gap[1][ipo95 - 1] * cof;
    suem = gap[2][ipo95 - 1] * cimu;
    gapr[0][ipo95 - 1] = (sume - suem).real();
    gapr[1][ipo95 - 1] = ((sume + suem) * uim).real();
    gapr[2][ipo95 - 1] = suee.real();
    sumep = gapp[0][ipo95 - 1] * cimu;
    sueep = gapp[1][ipo95 - 1] * cof;
    suemp = gapp[2][ipo95 - 1] * cimu;
    gapp[0][ipo95 - 1] = sumep - suemp;
    gapp[1][ipo95 - 1] = (sumep + suemp) * uim;
    gapp[2][ipo95 - 1] = sueep;
  } // ipo95 loop
  // Clean memory
  for (int ai = nlemt - 1; ai > -1; ai--) delete[] ac[ai];
  for (int gi = 2; gi > -1; gi--) delete[] gap[gi];
  delete[] ac;
  delete[] gap;
}

void apcra(
	   double ****zpv, const int le, complex<double> **am0m, int inpol, double sqk,
	   double **gaprm, complex<double> **gappm
) {
  const complex<double> cc0(0.0, 0.0);
  const complex<double> uim(0.0, 1.0);
  complex<double> uimtl, uimtls, ca11, ca12, ca21, ca22;
  complex<double> a11, a12, a21, a22, sum1, sum2, fc;
  double ****svw = new double***[le];
  complex<double> ****svs = new complex<double>***[le];
  for (int i = 0; i < le; i++) {
    svw[i] = new double**[3];
    svs[i] = new complex<double>**[3];
    for (int j = 0; j < 3; j++) {
      svw[i][j] = new double*[2];
      svs[i][j] = new complex<double>*[2];
      for (int k = 0; k < 2; k++) {
	svw[i][j][k] = new double[2]();
	svs[i][j][k] = new complex<double>[2]();
      }
    }
  }
  int nlem = le * (le + 2);
  for (int l28 = 1; l28 <= le; l28++) {
    int lpo = l28 + 1;
    int ltpo = lpo + l28;
    double fl = sqrt(1.0 * ltpo);
    for (int ilmp = 1; ilmp <= 3; ilmp++) {
      if ((l28 == 1 && ilmp == 1) || (l28 == le && ilmp == 3)) continue; // ilmp loop
      int lmpml = ilmp - 2;
      int lmp = l28 + lmpml;
      double flmp = sqrt(1.0 * (lmp + lmp + 1));
      double fllmp = flmp / fl;
      double cgmmo = fllmp * cg1(lmpml, 0, l28, 1);
      double cgmpo = fllmp * cg1(lmpml, 0, l28, -1);
      if (inpol == 0) {
	double cgs = cgmpo + cgmmo;
	double cgd = cgmpo - cgmmo;
	svw[l28 - 1][ilmp - 1][0][0] = cgs;
	svw[l28 - 1][ilmp - 1][0][1] = cgd;
	svw[l28 - 1][ilmp - 1][1][0] = cgd;
	svw[l28 - 1][ilmp - 1][1][1] = cgs;
      } else { // label 22
	svw[l28 - 1][ilmp - 1][0][0] = cgmpo;
	svw[l28 - 1][ilmp - 1][1][0] = cgmpo;
	svw[l28 - 1][ilmp - 1][0][1] = -cgmmo;
	svw[l28 - 1][ilmp - 1][1][1] = cgmmo;
      }
      // label 26
    } // ilmp loop
  } // l28 loop
  for (int l30 = 1; l30 <= le; l30++) { // 0-init: can be omitted
    for (int ilmp = 1; ilmp <= 3; ilmp++) {
      for (int ipa = 1; ipa <= 2; ipa++) {
	for (int ipamp = 1; ipamp <= 2; ipamp++) {
	  svs[l30 - 1][ilmp - 1][ipa - 1][ipamp - 1] = cc0;
	}
      } // ipa loop
    } // ilmp loop
  } // l30 loop
  for (int l58 = 1; l58 <= le; l58 ++) {
    int lpo = l58 + 1;
    int ltpo = l58 + lpo;
    int imm = l58 * lpo;
    for (int ilmp = 1; ilmp <= 3; ilmp++) {
      if ((l58 == 1 && ilmp == 1) || (l58 == le && ilmp == 3)) continue; // ilmp loop
      int lmpml = ilmp - 2;
      int lmp = l58 + lmpml;
      int impmm = lmp * (lmp + 1);
      uimtl = uim * (1.0 * lmpml);
      if (lmpml == 0) uimtl = complex<double>(1.0, 0.0);
      for (int im54 = 1; im54 <= ltpo; im54++) {
	int m = im54 - lpo;
	int i = imm + m;
	int ie = i + nlem;
	for (int imu52 = 1; imu52 <= 3; imu52++) {
	  int mu = imu52 - 2;
	  int mmp = m - mu;
	  int abs_mmp = (mmp > 0) ? mmp : -mmp;
	  if (abs_mmp <= lmp) {
	    int imp = impmm + mmp;
	    int impe = imp + nlem;
	    double cgc = cg1(lmpml, -mu, l58, -m);
	    for (int ls = 1; ls <= le; ls++) {
	      int lspo = ls + 1;
	      int lstpo = ls + lspo;
	      int ismm = ls * lspo;
	      for (int ilsmp = 1; ilsmp <= 3; ilsmp++) {
		if ((ls == 1 && ilsmp == 1) || (ls == le && ilsmp == 3)) continue; // ilsmp loop
		int lsmpml = ilsmp - 2;
		int lsmp = ls + lsmpml;
		int ismpmm = lsmp * (lsmp + 1);
		uimtls = -uim * (1.0 * lsmpml);
		if (lsmpml == 0) uimtls = complex<double>(1.0, 0.0);
		for (int ims = 1; ims <= lstpo; ims++) {
		  int ms = ims - lspo;
		  int msmp = ms - mu;
		  int abs_msmp = (msmp > 0) ? msmp : -msmp;
		  if (abs_msmp <= lsmp) {
		    int is = ismm + ms;
		    int ise = is + nlem;
		    int ismp = ismpmm + msmp;
		    int ismpe = ismp + nlem;
		    double cgcs = cg1(lsmpml, mu, ls, ms);
		    fc = (uimtl * uimtls) * (cgc * cgcs);
		    ca11 = dconjg(am0m[is - 1][i - 1]);
		    ca12 = dconjg(am0m[is - 1][ie - 1]);
		    ca21 = dconjg(am0m[ise - 1][i - 1]);
		    ca22 = dconjg(am0m[ise - 1][ie - 1]);
		    a11 = am0m[ismp - 1][imp - 1];
		    a12 = am0m[ismp - 1][impe - 1];
		    a21 = am0m[ismpe - 1][imp - 1];
		    a22 = am0m[ismpe - 1][impe - 1];
		    double z11 = zpv[ls - 1][ilsmp - 1][0][0];
		    double z12 = zpv[ls - 1][ilsmp - 1][0][1];
		    double z21 = zpv[ls - 1][ilsmp - 1][1][0];
		    double z22 = zpv[ls - 1][ilsmp - 1][1][1];
		    svs[l58 - 1][ilmp - 1][0][0] += ((ca11 * a11 * z11
						      + ca11 * a21 * z12
						      + ca21 * a11 * z21
						      + ca21 * a21 * z22) * fc);
		    svs[l58 - 1][ilmp - 1][0][1] += ((ca11 * a12 * z11
						      + ca11 * a22 * z12
						      + ca21 * a12 * z21
						      + ca21 * a22 * z22) * fc);
		    svs[l58 - 1][ilmp - 1][1][0] += ((ca12 * a11 * z11
						      + ca12 * a21 * z12
						      + ca22 * a11 * z21
						      + ca22 * a21 * z22) * fc);
		    svs[l58 - 1][ilmp - 1][1][1] += ((ca12 * a12 * z11
						      + ca12 * a22 * z12
						      + ca22 * a12 * z21
						      + ca22 * a22 * z22) * fc);
		  } // ends ims loop
		} // ims loop
	      } // ilsmp loop
	    } // ls loop
	  } // ends imu52 loop
	} // imu52 loop
      } // im54 loop
    } // ilmp loop
  } // l58 loop
  sum1 = cc0;
  sum2 = cc0;
  for (int l68 = 1; l68 <= le; l68++) {
    //int lpo = l68 + 1;
    //int ltpo = l68 + lpo;
    //int imm = l68 * lpo;
    for (int ilmp = 1; ilmp <= 3; ilmp++) {
      if ((l68 == 1 && ilmp == 1) || (l68 == le && ilmp == 3)) continue; // ilmp loop
      if (inpol == 0) {
	sum1 += (
		 svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][0][0]
		 + svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][0][1]
		 + svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][1][0]
		 + svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][1][1]
		 );
	sum2 += (
		 svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][0][0]
		 + svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][0][1]
		 + svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][1][0]
		 + svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][1][1]
		 );
      } else { // label 62
	sum1 += (
		 svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][0][0]
		 + svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][0][1]
		 + svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][1][0]
		 + svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][1][1]
		 );
	sum2 += (
		 svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][0][0]
		 + svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][0][1]
		 + svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][1][0]
		 + svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][1][1]
		 );
      } // label 66, ends ilmp loop
    } // ilmp loop
  } // l68 loop
  const double half_pi = acos(0.0);
  double cofs = half_pi * 2.0 / sqk;
  gaprm[0][0] = 0.0;
  gaprm[0][1] = 0.0;
  gaprm[1][0] = 0.0;
  gaprm[1][1] = 0.0;
  gappm[0][0] = cc0;
  gappm[0][1] = cc0;
  gappm[1][0] = cc0;
  gappm[1][1] = cc0;
  if (inpol == 0) {
    sum1 *= cofs;
    sum2 *= cofs;
    gaprm[2][0] = sum1.real();
    gaprm[2][1] = sum1.real();
    gappm[2][0] = sum2 * uim;
    gappm[2][1] = -gappm[2][0];
  } else { // label 72
    cofs *= 2.0;
    gaprm[2][0] = sum1.real() * cofs;
    gaprm[2][1] = sum2.real() * cofs;
    gappm[2][0] = cc0;
    gappm[2][1] = cc0;
  }
  
  // Clean memory
  for (int i = le - 1; i > -1; i--) {
    for (int j = 2; j > -1; j--) {
      for (int k = 1; k > -1; k--) {
	delete[] svw[i][j][k];
	delete[] svs[i][j][k];
      }
      delete[] svw[i][j];
      delete[] svs[i][j];
    }
    delete[] svw[i];
    delete[] svs[i];
  }
  delete[] svw;
  delete[] svs;
}

complex<double> cdtp(
		     complex<double> z, complex<double> **am, int i, int jf,
		     int k, int nj
) {
  /* NOTE: the original FORTRAN code treats the AM matrix as a
   * vector. This is not directly allowed in C++ and it requires
   * accounting for the different dimensions.
   */
  complex<double> result = z;
  if (nj > 0) {
    int jl = jf + nj - 1;
    for (int j = jf; j <= jl; j++) {
      result += (am[i - 1][j - 1] * am[j - 1][k - 1]);
    }
  }
  return result;
}

double cgev(int ipamo, int mu, int l, int m) {
  double result = 0.0;
  double xd = 0.0, xn = 0.0;
  if (ipamo == 0) {
    if (m != 0 || mu != 0) { // label 10
      if (mu != 0) {
	xd = 2.0 * l * (l + 1);
	if (mu <= 0) {
	  xn = 1.0 * (l + m) * (l - m + 1);
	  result = sqrt(xn / xd);
	} else { // label 15
	  xn = 1.0 * (l - m) * (l + m + 1);
	  result = -sqrt(xn / xd);
	}
      } else { // label 20
	xd = 1.0 * (l + 1) * l;
	xn = -1.0 * m;
	result = xn / sqrt(xd);
      }
    }
  } else { // label 30
    xd = 2.0 * l * (l * 2 - 1);
    if (mu < 0) { // label 35
      xn = 1.0 * (l - 1 + m) * (l + m);
    } else if (mu == 0) { // label 40
      xn = 2.0 * (l - m) * (l + m);
    } else { // mu > 0, label 45
      xn = 1.0 * (l - 1 - m) * (l - m);
    }
    result = sqrt(xn / xd);
  }
  return result;
}

void cms(complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
  complex<double> dm, de, cgh, cgk;
  const complex<double> cc0(0.0, 0.0);
  int ndi = c4->nsph * c4->nlim;
  int nbl = 0;
  int nsphmo = c4->nsph - 1;
  for (int n1 = 1; n1 <= nsphmo; n1++) { // GPU portable?
    int in1 = (n1 - 1) * c4->nlim;
    int n1po = n1 + 1;
    for (int n2 = n1po; n2 <= c4->nsph; n2++) {
      int in2 = (n2 - 1) * c4->nlim;
      nbl++;
      for (int l1 = 1; l1 <= c4->li; l1++) {
	int l1po = l1 + 1;
	int il1 = l1po * l1;
	int l1tpo = l1po + l1;
	for (int im1 = 1; im1 <= l1tpo; im1++) {
	  int m1 = im1 - l1po;
	  int ilm1 = il1 + m1;
	  int ilm1e = ilm1 + ndi;
	  int i1 = in1 + ilm1;
	  int i1e = in1 + ilm1e;
	  int j1 = in2 + ilm1;
	  int j1e = in2 + ilm1e;
	  for (int l2 = 1; l2 <= c4->li; l2++) {
	    int l2po = l2 + 1;
	    int il2 = l2po * l2;
	    int l2tpo = l2po + l2;
	    int ish = ((l2 + l1) % 2 == 0) ? 1 : -1;
	    int isk = -ish;
	    for (int im2 = 1; im2 <= l2tpo; im2++) {
	      int m2 = im2 - l2po;
	      int ilm2 = il2 + m2;
	      int ilm2e = ilm2 + ndi;
	      int i2 = in2 + ilm2;
	      int i2e = in2 + ilm2e;
	      int j2 = in1 + ilm2;
	      int j2e = in1 + ilm2e;
	      cgh = ghit(0, 0, nbl, l1, m1, l2, m2, c1, c1ao, c4, c6);
	      cgk = ghit(0, 1, nbl, l1, m1, l2, m2, c1, c1ao, c4, c6);
	      am[i1 - 1][i2 - 1] = cgh;
	      am[i1 - 1][i2e - 1] = cgk;
	      am[i1e - 1][i2 - 1] = cgk;
	      am[i1e - 1][i2e - 1] = cgh;
	      am[j1 - 1][j2 - 1] = cgh * (1.0 * ish);
	      am[j1 - 1][j2e - 1] = cgk * (1.0 * isk);
	      am[j1e - 1][j2 - 1] = cgk * (1.0 * isk);
	      am[j1e - 1][j2e - 1] = cgh * (1.0 * ish);
	    }
	  }
	} // im1 loop
      } // l1 loop
    } // n2 loop
  } // n1 loop
  for (int n1 = 1; n1 <= c4->nsph; n1++) { // GPU portable?
    int in1 = (n1 - 1) * c4->nlim;
    for (int l1 = 1; l1 <= c4->li; l1++) {
      dm = c1->rmi[l1 - 1][n1 - 1];
      de = c1->rei[l1 - 1][n1 - 1];
      int l1po = l1 + 1;
      int il1 = l1po * l1;
      int l1tpo = l1po + l1;
      for (int im1 = 1; im1 <= l1tpo; im1++) {
	int m1 = im1 - l1po;
	int ilm1 = il1 + m1;
	int i1 = in1 + ilm1;
	int i1e = i1 + ndi;
	for (int ilm2 = 1; ilm2 <= c4->nlim; ilm2++) {
	  int i2 = in1 + ilm2;
	  int i2e = i2 + ndi;
	  am[i1 - 1][i2 - 1] = cc0;
	  am[i1 - 1][i2e - 1] = cc0;
	  am[i1e - 1][i2 - 1] = cc0;
	  am[i1e - 1][i2e - 1] = cc0;
	}
	am[i1 - 1][i1 - 1] = dm;
	am[i1e - 1][i1e - 1] = de;
      } // im1 loop
    } // l1 loop
  } // n1 loop
}

void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
  complex<double> ***svf, ***svw, **svs;
  const complex<double> cc0(0.0, 0.0);
  complex<double> cam(0.0, 0.0);
  const int le4po = 4 * c4->le + 1;
  svf = new complex<double>**[le4po];
  svw = new complex<double>**[le4po];
  svs = new complex<double>*[le4po];
  for (int si = 0; si < le4po; si++) {
    svf[si] = new complex<double>*[le4po];
    svw[si] = new complex<double>*[4];
    svs[si] = new complex<double>[4]();
    for (int sj = 0; sj < le4po; sj++) svf[si][sj] = new complex<double>[4]();
    for (int sj = 0; sj < 4; sj++) svw[si][sj] = new complex<double>[4]();
  }
  double exdc = exri * exri;
  double ccs = 1.0 / (vk * vk);
  const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
  double cint = ccs / (pi4sq * exdc);
  int letpo = c4->le + c4->le + 1;
  for (int i20 = 0; i20 < 16; i20++) c1ao->vintm[i20] = cc0; // 0-init: can be omitted
  for (int lpo40 = 1; lpo40 <= letpo; lpo40++) {
    int l = lpo40 - 1;
    int ltpo = lpo40 + l;
    int immn = letpo - l;
    int immx = letpo + l;
    for (int imf = immn; imf <= immx; imf++) { // 0-init: can be omitted
      for (int ims = immn; ims <= immx; ims++) {
	for (int ipo = 1; ipo <= 4; ipo++) {
	  svf[imf - 1][ims - 1][ipo - 1] = cc0;
	} // ipo loop
      } // ims loop
    } // imf loop
    for (int l1 = 1; l1 <= c4->le; l1++) {
      int il1 = l1 * (l1 + 1);
      for (int l2 = 1; l2 <= c4->le; l2++) {
	int abs_l2ml1 = (l2 > l1) ? l2 - l1 : l1 - l2;
	if (l < abs_l2ml1 || l > l2 + l1) continue; // l2 loop
	int il2 = l2 * (l2 + 1);
	for (int im = immn; im >= immx; im++) { // 0-init: can be omitted
	  for (int ipa = 1; ipa <= 4; ipa++) {
	    svs[im - 1][ipa - 1] = cc0;
	    for (int ipo = 1; ipo <= 4; ipo++) {
	      svw[im - 1][ipa - 1][ipo - 1] = cc0;
	    } // ipo loop
	  } // ipa loop
	} // im loop
	for (int im = immn; im <= immx; im++) {
	  int m = im - letpo;
	  r3jmr(l, l1, l2, m, c6);
	  int m1mnmo = (-l1 > -l2 - m) ? -(l1 + 1) : -(l2 + m + 1);
	  int nm1 = (l1 < l2 - m) ? (l1 - m1mnmo) : (l2 - m - m1mnmo);
	  for (int im1 = 1; im1 <= nm1; im1++) {
	    int m1 = -im1 - m1mnmo;
	    int isn = 1;
	    if (m1 % 2 != 0) isn = -1;
	    double cg3j = c6->rac3j[im1 - 1] * isn;
	    int ilm1 = il1 + m1;
	    int ilm2 = il2 + m1 - m;
	    int ipa = 0;
	    for (int ipa1 = 1; ipa1 <= 2; ipa1++) {
	      int i1 = ilm1;
	      if (ipa1 == 2) i1 = ilm1 + c4->nlem;
	      for (int ipa2 = 1; ipa2 <= 2; ipa2++) {
		int i2 = ilm2;
		if (ipa2 == 2) i2 = ilm2 + c4->nlem;
		ipa++;
		svs[im - 1][ipa - 1] += (c1ao->am0m[i1 - 1][i2 - 1] * cg3j);
		int ipo = 0;
		for (int ipo2 = 1; ipo2 <= 2; ipo2++) {
		  for (int ipo1 = 3; ipo1 <= 4; ipo1++) {
		    ipo++;
		    svw[im - 1][ipa - 1][ipo - 1] += (c1->w[i1 - 1][ipo1 - 1] * c1->w[i2 - 1][ipo2 - 1] * cg3j);
		  } // ipo1 loop
		} // ipo2 loop
	      } // ipa2 loop
	    } // ipa1 loop
	  } // im1 loop
	  // label 32 loops
	  for (int imf = immn; imf <= immx; imf++) {
	    for (int ims = immn; ims <= immx; ims++) {
	      for (int ipo = 1; ipo <= 4; ipo++) {
		for (int ipa = 1; ipa <= 4; ipa++) {
		  svf[imf - 1][ims - 1][ipo - 1] += (svw[imf - 1][ipa - 1][ipo - 1] * svs[ims - 1][ipa - 1]);
		} // ipa loop
	      } // ipo loop
	    } // ims loop
	  } // imf loop
	  // ends loop level 34, which are l2 loop and l1 loop
	} // im loop
      } // l2 loop
    } // l1 loop
    for (int imf = immn; imf <= immx; imf++) {
      for (int ims = immn; ims <= immx; ims++) {
	int i = 0;
	for (int ipo1 = 1; ipo1 <= 4; ipo1++) {
	  cam = dconjg(svf[imf - 1][ims - 1][ipo1 - 1]);
	  for (int ipo2 = 1; ipo2 <= 4; ipo2++) {
	    i++;
	    c1ao->vintm[i - 1] += (svf[imf - 1][ims - 1][ipo2 - 1] * cam * (1.0 * ltpo));
	  }
	} // ipo1 loop
      } // ims loop
    } // imf loop
  } // lpo40 loop
  for (int i42 = 0; i42 < 16; i42++) c1ao->vintm[i42] *= cint;

  // Clean memory
  for (int si = le4po - 1; si > -1; si--) {
    for (int sj = le4po - 1; sj > -1; sj--) delete[] svf[si][sj];
    for (int sj = 3; sj > -1; sj--) delete[] svw[si][sj];
    delete[] svf[si];
    delete[] svw[si];
    delete[] svs[si];
  }
  delete[] svf;
  delete[] svw;
  delete[] svs;
}

complex<double> ghit(
		     int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
		     C1_AddOns *c1ao, C4 *c4, C6 *c6
) {
  /* NBL identifies transfer vector going from N2 to N1;
   * IHI=0 for Hankel, IHI=1 for Bessel, IHI=2 for Bessel from origin;
   * depending on IHI, IPAM=0 gives H or I, IPAM= 1 gives K or L. */
  const complex<double> cc0(0.0, 0.0);
  const complex<double> uim(0.0, 1.0);
  complex<double> csum(0.0, 0.0), cfun(0.0, 0.0);
  complex<double> result = cc0;

  if (ihi == 2) {
    if (c1->rxx[nbl - 1] == 0.0 && c1->ryy[nbl - 1] == 0.0 && c1->rzz[nbl - 1] == 0.0) {
      if (ipamo == 0) {
	if (l1 == l2 && m1 == m2) result = complex(1.0, 0.0);
      }
      return result;
    }
  }
  // label 10
  int l1mp = l1 - ipamo;
  int l1po = l1 + 1;
  int m1mm2 = m1 - m2;
  int m1mm2m = (m1mm2 > 0) ? m1mm2 + 1 : 1 - m1mm2;
  int lminpo = (l2 - l1mp > 0) ? l2 - l1mp + 1 : l1mp - l2 + 1;
  int lmaxpo = l2 + l1mp + 1;
  int i3j0in = c1ao->ind3j[l1mp][l2 - 1];
  int ilin = -1;
  if (m1mm2m > lminpo && (m1mm2m - lminpo) % 2 != 0) ilin = 0;
  int isn = 1;
  if (m1 % 2 != 0) isn *= -1;
  if (lminpo % 2 == 0) {
    isn *= -1;
    if (l2 > l1mp) isn *= -1;
  }
  // label 12
  int nblmo = nbl - 1;
  if (ihi != 2) {
    int nbhj = nblmo * c4->litpo;
    int nby = nblmo * c4->litpos;
    if (ihi != 1) {
      for (int jm24 = 1; jm24 <= 3; jm24++) {
	csum = cc0;
	int mu = jm24 - 2;
	int mupm1 = mu + m1;
	int mupm2 = mu + m2;
	if (mupm1 >= -l1mp && mupm1 <= l1mp && mupm2 >= - l2 && mupm2 <= l2) {
	  int jsn = -isn;
	  if (mu == 0) jsn = isn;
	  double cr = cgev(ipamo, mu, l1, m1) * cgev(0, mu, l2, m2);
	  int i3j0 = i3j0in;
	  if (mupm1 == 0 && mupm2 == 0) {
	    int lt14 = lminpo;
	    while (lt14 <= lmaxpo) {
	      i3j0++;
	      int l3 = lt14 - 1;
	      int ny = l3 * l3 + lt14;
	      double aors = 1.0 * (l3 + lt14);
	      double f3j = (c1ao->v3j0[i3j0 - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)) * jsn;
	      cfun = (c1ao->vh[nbhj + lt14 - 1] * c1ao->vyhj[nby + ny - 1]) * f3j;
	      csum += cfun;
	      jsn *= -1;
	      lt14 += 2;
	    }
	    // goes to 22
	  } else { // label 16
	    r3jjr(l1mp, l2, -mupm1, mupm2, c6);
	    int il = ilin;
	    int lt20 = lminpo;
	    while (lt20 <= lmaxpo) {
	      i3j0++;
	      if (m1mm2m <= lt20) {
		il += 2;
		int l3 = lt20 - 1;
		int ny = l3 * l3  + lt20 + m1mm2;
		double aors = 1.0 * (l3 + lt20);
		double f3j = (c6->rac3j[il - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)) * jsn;
		cfun = (c1ao->vh[nbhj + lt20 - 1] * c1ao->vyhj[nby + ny - 1]) * f3j;
		csum += cfun;
	      }
	      // label 20
	      jsn *= -1;
	      lt20 += 2;
	    }
	  }
	  // label 22
	  csum *= cr;
	  result += csum;
	}
	// Otherwise there is nothing to add
      } // jm24 loop. Should go to 70
    } else { // label 30, IHI == 1
      for (int jm44 = 1; jm44 <= 3; jm44++) {
	csum = cc0;
	int mu = jm44 - 2;
	int mupm1 = mu + m1;
	int mupm2 = mu + m2;
	if (mupm1 >= -l1mp && mupm1 <= l1mp && mupm2 >= - l2 && mupm2 <= l2) {
	  int jsn = - isn;
	  if (mu == 0) jsn = isn;
	  double cr = cgev(ipamo, mu, l1, m1) * cgev(0, mu, l2, m2);
	  int i3j0 = i3j0in;
	  if (mupm1 == 0 && mupm2 == 0) {
	    int lt34 = lminpo;
	    while (lt34 <= lmaxpo) {
	      i3j0++;
	      int l3 = lt34 - 1;
	      int ny = l3 * l3 + lt34;
	      double aors = 1.0 * (l3 + lt34);
	      double f3j = (c1ao->v3j0[i3j0 - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)) * jsn;
	      cfun = (c1ao->vh[nbhj + lt34 - 1] * c1ao->vyhj[nby + ny - 1]) * f3j;
	      csum += cfun;
	      jsn *= -1;
	      lt34 += 2;
	    }
	    // goes to 42
	  } else { // label 36
	    r3jjr(l1mp, l2, -mupm1, mupm2, c6);
	    int il = ilin;
	    int lt40 = lminpo;
	    while (lt40 <= lmaxpo) {
	      i3j0++;
	      if (m1mm2m <= lt40) {
		il += 2;
		int l3 = lt40 - 1;
		int ny = l3 * l3  + lt40 + m1mm2;
		double aors = 1.0 * (l3 + lt40);
		double f3j = (c6->rac3j[il - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)) * jsn;
		cfun = (c1ao->vh[nbhj + lt40 - 1] * c1ao->vyhj[nby + ny - 1]) * f3j;
		csum += cfun;
	      }
	      // label 40
	      jsn *= -1;
	      lt40 += 2;
	    }
	  }
	  // label 42
	  csum *= cr;
	  result += csum;
	}
	// Otherwise there is nothing to add
      } // jm44 loop. Should go to 70
    }
    // goes to 70
  } else { // label 50, IHI == 2
    int nbhj = nblmo * c4->lmtpo;
    int nby = nblmo * c4->lmtpos;
    for (int jm64 = 1; jm64 <= 3; jm64++) {
      csum = cc0;
      int mu = jm64 - 2;
      int mupm1 = mu + m1;
      int mupm2 = mu + m2;
      if (mupm1 >= -l1mp && mupm1 <= l1mp && mupm2 >= - l2 && mupm2 <= l2) {
	int jsn = -isn;
	if (mu == 0) jsn = isn;
	double cr = cgev(ipamo, mu, l1, m1) * cgev(0, mu, l2, m2);
	int i3j0 = i3j0in;
	if (mupm1 == 0 && mupm2 == 0) {
	  int lt54 = lminpo;
	  while (lt54 <= lmaxpo) {
	    i3j0++;
	    int l3 = lt54 - 1;
	    int ny = l3 * l3 + lt54;
	    double aors = 1.0 * (l3 + lt54);
	    double f3j = (c1ao->v3j0[i3j0 - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)) * jsn;
	    cfun = (c1ao->vj0[nbhj + lt54 - 1] * c1ao->vyj0[nby + ny - 1]) * f3j;
	    csum += cfun;
	    jsn *= -1;
	    lt54 += 2;
	  }
	  // goes to 62
	} else { // label 56
	  r3jjr(l1mp, l2, -mupm1, mupm2, c6);
	  int il = ilin;
	  int lt60 = lminpo;
	  while (lt60 <= lmaxpo) {
	    i3j0++;
	    if (m1mm2m <= lt60) {
	      il += 2;
	      int l3 = lt60 - 1;
	      int ny = l3 * l3  + lt60 + m1mm2;
	      double aors = 1.0 * (l3 + lt60);
	      double f3j = (c6->rac3j[il - 1] * c1ao->v3j0[i3j0 - 1] * sqrt(aors)) * jsn;
	      cfun = (c1ao->vj0[nbhj + lt60 - 1] * c1ao->vyj0[nby + ny - 1]) * f3j;
	      csum += cfun;
	    }
	    // label 60
	    jsn *= -1;
	    lt60 += 2;
	  }
	}
	// label 62
	csum *= cr;
	result += csum;
      }
      // Otherwise there is nothing to add
    } // jm64 loop. Should go to 70
  }
  // label 70
  const double four_pi = acos(0.0) * 8.0;
  if (ipamo != 1) {
    double cr = sqrt(four_pi * (l1 + l1po) * (l2 + l2 + 1));
    result *= cr;
  } else {
    double cr = sqrt(four_pi * (l1 + l1mp) * (l1 + l1po) * (l2 + l2 + 1) / l1po);
    result *= (cr * uim);
  }
  return result;
}

void hjv(
	 double exri, double vk, int &jer, int &lcalc, complex<double> &arg,
	 C1 *c1, C1_AddOns *c1ao, C4 *c4
) {
  int nsphmo = c4->nsph - 1;
  int lit = c4->li + c4->li;
  int lmt = c4->li + c4->le;
  const int rfj_size = (lit > lmt) ? lit : lmt;
  const int rfn_size = c4->litpo;
  double *rfj, *rfn;
  rfj = new double[rfj_size]();
  rfn = new double[rfn_size]();
  jer = 0;
  int ivhb = 0;
  for (int nf40 = 1; nf40 <= nsphmo; nf40++) { // GPU portable?
    int nfpo = nf40 + 1;
    for (int ns40 = nfpo; ns40 <= c4->nsph; ns40++) {
      double rx = c1->rxx[nf40 - 1] - c1->rxx[ns40 - 1];
      double ry = c1->ryy[nf40 - 1] - c1->ryy[ns40 - 1];
      double rz = c1->rzz[nf40 - 1] - c1->rzz[ns40 - 1];
      double rr = sqrt(rx * rx + ry * ry + rz * rz);
      double rarg = rr * vk * exri;
      arg = complex<double>(rarg, 0.0);
      rbf(lit, rarg, lcalc, rfj);
      if (lcalc < lit) {
	jer = 1;
	delete[] rfj;
	delete[] rfn;
	return;
      }
      rnf(lit, rarg, lcalc, rfn);
      if (lcalc < lit) {
	jer = 2;
	delete[] rfj;
	delete[] rfn;
	return;
      }
      for (int lpo38 = 1; lpo38 <= c4->litpo; lpo38++) {
	double rpart = rfj[lpo38 - 1];
	double ipart = rfn[lpo38 - 1];
	c1ao->vh[lpo38 + ivhb - 1] = complex<double>(rpart, ipart);
      }
      ivhb += c4->litpo;
    } // ns40 loop
  } // nf40 loop
  ivhb = 0;
  for (int nf50 = 1; nf50 <= c4->nsph; nf50++) {
    double rx = c1->rxx[nf50 - 1];
    double ry = c1->ryy[nf50 - 1];
    double rz = c1->rzz[nf50 - 1];
    if (!(rx == 0.0 && ry == 0.0 && rz == 0.0)) {
      double rr = sqrt(rx * rx + ry * ry + rz * rz);
      double rarg = rr * vk * exri;
      rbf(lmt, rarg, lcalc, rfj);
      if (lcalc < lmt) {
	jer = 3;
	delete[] rfj;
	delete[] rfn;
	return;
      }
      for (int lpo47 = 1; lpo47 <= c4->lmtpo; lpo47++) {
	c1ao->vj0[lpo47 + ivhb - 1] = rfj[lpo47 - 1];
      }
    }
    ivhb += c4->lmtpo;
  } // nf50 loop
  delete[] rfj;
  delete[] rfn;
}

void lucin(complex<double> **am, const int nddmst, int n, int &ier) {
  /* NDDMST  FIRST DIMENSION OF AM AS DECLARED IN DIMENSION
   *         STATEMENT.
   * N       NUMBER OF ROWS IN AM.
   * IER     IS REPLACED BY 1 FOR SINGULARITY.
   */
  double *v = new double[nddmst];
  complex<double> ctemp, cfun;
  complex<double> cc0 = complex<double>(0.0, 0.0);
  ier = 0;
  int nminus = n - 1;
  for (int i = 1; i <= n; i++) {
    double sum = 0.0;
    for (int j = 1; j <= n; j++) {
      sum += (
	      am[i - 1][j - 1].real() * am[i - 1][j - 1].real()
	      + am[i - 1][j - 1].imag() * am[i - 1][j - 1].imag()
	      );
    } // j1319 loop
    v[i - 1] = 1.0 / sum;
  } // i1309 loop
  // 2. REPLACE AM BY TRIANGULAR MATRICES (L,U) WHERE AM=L*U.
  //    REPLACE L(I,I) BY 1/L(I,I), READY FOR SECTION 4.
  //    (ROW INTERCHANGES TAKE PLACE, AND THE INDICES OF THE PIVOTAL ROWS
  //    ARE PLACED IN V.)
  /* >>> THERE APPEARS TO BE A BUG IN THE FOLLOWING LOOP <<< */
  for (int k = 1; k <= n; k++) {
    int kplus = k + 1;
    int kminus = k - 1;
    int l = k;
    double psqmax = 0.0;
    for (int i = k; i <= n; i++) {
      cfun = cdtp(-am[i - 1][k - 1], am, i, 1, k, kminus);
      ctemp = -cfun;
      am[i - 1][k - 1] = ctemp;
      double psq = v[i - 1] * (ctemp.real() * ctemp.real() + ctemp.imag() * ctemp.imag());
      if (psq > psqmax) {
	psqmax = psq;
	l = i;
      }
    } // i2029 loop
    if (l != k) {
      for (int j = 1; j <= n; j++) {
	ctemp = am[k - 1][j - 1];
	am[k - 1][j - 1] = am[l - 1][j - 1];
	am[l - 1][j - 1] = ctemp;
      } // j2049 loop
      v[l - 1] = v[k - 1];
    }
    // label 2011
    v[k - 1] = 1.0 * l;
    if (psqmax == 0.0) {
      ier = 1;
      delete[] v;
      return;
    }
    ctemp = 1.0 / am[k - 1][k - 1];
    am[k - 1][k - 1] = ctemp;
    if (kplus <= n) {
      for (int j = kplus; j <= n; j++) {
	cfun = cdtp(-am[k - 1][j - 1], am, k, 1, j, kminus);
	am[k - 1][j - 1] = -ctemp * cfun;
      } // j2059 loop
    }
  } // k2019 loop
  // 4.  REPLACE AM BY ITS INVERSE AMINV.
  // 4.1 REPLACE L AND U BY THEIR INVERSE LINV AND UINV.
  for (int k = 1; k <= nminus; k++) {
    int kplus = k + 1;
    for (int i = kplus; i <= n; i++) {
      cfun = cdtp(cc0, am, i, k, k, i - k);
      am[i - 1][k - 1] = -am[i - 1][i - 1] * cfun;
      cfun = cdtp(am[k - 1][i - 1], am, k, kplus, i, i - k - 1);
      am[k - 1][i - 1] = -cfun;
    } // i4119 loop
  } // k4109 loop
  // 4.2 FORM AMINV=UINV*LINV.
  for (int k = 1; k <= n; k++) {
    for (int i = 1; i <= n; i++) {
      if (i < k) {
	cfun = cdtp(cc0, am, i, k, k, n - k + 1);
	am[i - 1][k -1] = cfun;
      }
      else {
	cfun = cdtp(am[i - 1][k - 1], am, i, i + 1, k, n - i);
	am[i - 1][k - 1] = cfun;
      }
    } // i4119 loop
  } // k4209 loop
  // 4.3 INTERCHANGE COLUMNS OF AMINV AS SPECIFIED BY V, BUT IN REVERSE
  //     ORDER.
  for (int l = 1; l <= n; l++) {
    int k = n - l + 1;
    int kcol = (int)(v[k - 1]);
    if (kcol != k) {
      for (int i = 1; i <= n; i++) {
	ctemp = am[i - 1][k - 1];
	am[i - 1][k - 1] = am[i - 1][kcol - 1];
	am[i - 1][kcol - 1] = ctemp;
      } // i4319 loop
    }
  } // l4309 loop
  delete[] v;
}

void mextc(double vk, double exri, complex<double> **fsac, double **cextlr, double **cext) {
  double fa11r = fsac[0][0].real();
  double fa11i = fsac[0][0].imag();
  double fa21r = fsac[1][0].real();
  double fa21i = fsac[1][0].imag();
  double fa12r = fsac[0][1].real();
  double fa12i = fsac[0][1].imag();
  double fa22r = fsac[1][1].real();
  double fa22i = fsac[1][1].imag();
  cextlr[0][0] = fa11i * 2.0;
  cextlr[0][1] = 0.0;
  cextlr[0][2] = -fa12i;
  cextlr[0][3] = -fa12r;
  cextlr[1][0] = 0.0;
  cextlr[1][1] = fa22i * 2.0;
  cextlr[1][2] = -fa21i;
  cextlr[1][3] = fa21r;
  cextlr[2][0] = -fa21i * 2.0;
  cextlr[2][1] = -fa12i * 2.0;
  cextlr[2][2] = fa11i + fa22i;
  cextlr[2][3] = fa22r - fa11r;
  cextlr[3][0] = fa21r * 2.0;
  cextlr[3][1] = -fa12r * 2.0;
  cextlr[3][2] = fa11r - fa22r;
  cextlr[3][3] = cextlr[2][2];
  cext[0][0] = cextlr[3][3];
  cext[1][1] = cextlr[3][3];
  cext[2][2] = cextlr[3][3];
  cext[2][3] = cextlr[2][3];
  cext[3][2] = cextlr[3][2];
  cext[3][3] = cextlr[3][3];
  cext[0][1] = fa11i - fa22i;
  cext[0][2] = -fa12i - fa21i;
  cext[0][3] = fa21r - fa12r;
  cext[1][0] = cext[0][1];
  cext[1][2] = fa21i - fa12i;
  cext[3][1] = fa12r + fa21r;
  cext[1][3] = -cext[3][1];
  cext[2][0] = cext[0][2];
  cext[2][1] = -cext[1][2];
  cext[3][0] = cext[1][3];
  double ckm = vk / exri;
  for (int i10 = 0; i10 < 4; i10++) {
    for (int j10 = 0; j10 < 4; j10++) {
      cextlr[i10][j10] *= ckm;
      cext[i10][j10] *= ckm;
    }
  }
}

void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
  const complex<double> cc0(0.0, 0.0);
  complex<double> sump, sum1, sum2, sum3, sum4, am, amp, cc, csam;
  const double exdc = exri * exri;
  double ccs = 1.0 / (vk * vk);
  double cccs = ccs / exdc;
  csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
  const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
  double cfsq = 4.0 / (pi4sq *ccs * ccs);
  const int nlemt = c4->nlem + c4->nlem;
  int jpo = 2;
  for (int ipo18 = 1; ipo18 <= 2; ipo18++) {
    if (ipo18 == 2) jpo = 1;
    int ipopt = ipo18 + 2;
    int jpopt = jpo + 2;
    double sum = 0.0;
    sump = cc0;
    sum1 = cc0;
    sum2 = cc0;
    sum3 = cc0;
    sum4 = cc0;
    for (int i12 = 1; i12 <= nlemt; i12++) {
      int i = i12 - 1;
      am = cc0;
      amp = cc0;
      for (int j10 = 1; j10 <= nlemt; j10++) {
	int j = j10 - 1;
	am += (c1ao->am0m[i][j] * c1->w[j][ipo18 - 1]);
	amp += (c1ao->am0m[i][j] * c1->w[j][jpo - 1]);
      } // j10 loop
      sum += (dconjg(am) * am).real();
      sump += (dconjg(amp) * am);
      sum1 += (dconjg(c1->w[i][ipo18 - 1]) * am);
      sum2 += (dconjg(c1->w[i][jpo - 1]) * am);
      sum3 += (c1->w[i][ipopt - 1] * am);
      sum4 += (c1->w[i][jpopt - 1] * am);
    } // i12 loop
    c1ao->scsc[ipo18 - 1] = cccs * sum;
    c1ao->scscp[ipo18 - 1] = cccs * sump;
    c1ao->ecsc[ipo18 - 1] = -cccs * sum1.real();
    c1ao->ecscp[ipo18 - 1] = -cccs * sum2;
    c1ao->fsac[ipo18 - 1][ipo18 - 1] = csam * sum1;
    c1ao->fsac[jpo - 1][ipo18 - 1] = csam * sum2;
    c1ao->sac[ipo18 - 1][ipo18 - 1] = csam * sum3;
    c1ao->sac[jpo - 1][ipo18 - 1] = csam * sum4;
  } // ipo18 loop
  int i = 0;
  for (int ipo1 = 1; ipo1 <= 2; ipo1++) {
    for (int jpo1 = 1; jpo1 <= 2; jpo1++) {
      cc = dconjg(c1ao->sac[jpo1 - 1][ipo1 - 1]);
      for (int ipo2 = 1; ipo2 <= 2; ipo2 ++) {
	for (int jpo2 = 1; jpo2 <= 2; jpo2++) {
	  c1ao->vint[i++] = c1ao->sac[jpo2 - 1][ipo2 - 1] * cc * cfsq;
	} // jpo2 loop
      } // ipo2 loop
    } // jpo1 loop
  } // ipo1 loop
}

void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
  const complex<double> cc0(0.0, 0.0);
  const complex<double> uim(0.0, 1.0);
  complex<double> sum1, sum2, sum3, sum4, sumpd;
  complex<double> sums1, sums2, sums3, sums4, csam;
  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);
  sum2 = cc0;
  sum3 = cc0;
  for (int i14 = 1; i14 <= c4->nlem; i14++) { // GPU portable?
    int ie = i14 + c4->nlem;
    sum2 += (c1ao->am0m[i14 - 1][i14 - 1] + c1ao->am0m[ie - 1][ie - 1]);
    sum3 += (c1ao->am0m[i14 - 1][ie - 1] + c1ao->am0m[ie - 1][i14 - 1]);
  } // i14 loop
  double sumpi = 0.0;
  sumpd = cc0;
  int nlemt = c4->nlem + c4->nlem;
  for (int i16 = 1; i16 <= nlemt; i16++) {
    for (int j16 = 1; j16 <= c4->nlem; j16++) {
      int je = j16 + c4->nlem;
      double rvalue = (
		       dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][j16 - 1]
		       + dconjg(c1ao->am0m[i16 - 1][je - 1]) * c1ao->am0m[i16 - 1][je - 1]
		       ).real();
      sumpi += rvalue;
      sumpd += (
		dconjg(c1ao->am0m[i16 - 1][j16 - 1]) * c1ao->am0m[i16 - 1][je - 1]
		+ dconjg(c1ao->am0m[i16 - 1][je - 1]) * c1ao->am0m[i16 - 1][j16 - 1]
		);
    } // j16 loop
  } // i16 loop
  if (inpol == 0) {
    sum1 = sum2;
    sum4 = sum3 * uim;
    sum3 = -sum4;
    sums1 = sumpi;
    sums2 = sumpi;
    sums3 = sumpd * uim;
    sums4 = -sums3;
  } else { // label 18
    sum1 = sum2 + sum3;
    sum2 = sum2 - sum3;
    sum3 = cc0;
    sum4 = cc0;
    sums1 = sumpi - sumpd;
    sums2 = sumpi + sumpd;
    sums3 = cc0;
    sums4 = cc0;
  }
  // label 20
  c1ao->ecscm[0] = -cccs * sum2.real();
  c1ao->ecscm[1] = -cccs * sum1.real();
  c1ao->ecscpm[0] = -cccs * sum4;
  c1ao->ecscpm[1] = -cccs * sum3;
  c1ao->fsacm[0][0] = csam * sum2;
  c1ao->fsacm[1][0] = csam * sum4;
  c1ao->fsacm[1][1] = csam * sum1;
  c1ao->fsacm[0][1] = csam * sum3;
  c1ao->scscm[0] = cccs * sums1.real();
  c1ao->scscm[1] = cccs * sums2.real();
  c1ao->scscpm[0] = cccs * sums3;
  c1ao->scscpm[1] = cccs * sums4;
}

void polar(
	   double x, double y, double z, double &r, double &cth, double &sth,
	   double &cph, double &sph
) {
  bool onx = (y == 0.0);
  bool ony = (x == 0.0);
  bool onz = (onx && ony);
  double rho = 0.0;
  if (!onz) {
    if (!onx) {
      if (!ony) {
	rho = sqrt(x * x + y * y);
	cph = x / rho;
	sph = y / rho;
	// goes to 25
      } else { // label 20
	rho = (y > 0.0) ? y : -y;
	cph = 0.0;
	sph = (y > 0.0) ? 1.0 : -1.0;
	// goes to 25
      }
    } else { // label 15
      rho = (x > 0.0) ? x : -x;
      cph = (x > 0.0) ? 1.0 : -1.0;
      sph = 0.0;
      // goes to 25
    }
  } else { // label 10
    cph = 1.0;
    sph = 0.0;
    // goes to 25
  }
  // label 25
  if (z == 0.0) {
    if (!onz) {
      r = rho;
      cth = 0.0;
      sth = 1.0;
      // returns
    } else { // label 30
      r = 0.0;
      cth = 1.0;
      sth = 0.0;
      // returns
    }
  } else { // label 35
    if (!onz) {
      r = sqrt(rho * rho + z * z);
      cth = z / r;
      sth = rho / r;
      // returns
    } else { // label 40
      r = (z > 0.0) ? z : -z;
      cth = (z > 0.0) ? 1.0 : -1.0;
      sth = 0.0;
      // returns
    }
  }
}

void r3j000(int j2, int j3, C6 *c6) {
  int jmx = j3 + j2;
  if (jmx <= 0) {
    c6->rac3j[0] = 1.0;
    return;
  }
  int jmn = j3 - j2;
  if (jmn < 0) jmn *= -1;
  int njmo = (jmx - jmn) / 2;
  int jf = jmx + jmx + 1;
  int isn = 1;
  if (jmn % 2 != 0) isn = -1;
  if (njmo <= 0) {
    double sj = 1.0 * jf;
    double cnr = (1 / sqrt(sj)) * isn;
    c6->rac3j[0] = cnr;
    return;
  }
  double sjr = 1.0 * jf;
  int jmxpos = (jmx + 1) * (jmx + 1);
  int jmns = jmn * jmn;
  int j1mo = jmx - 1;
  int j1s = (j1mo + 1) * (j1mo + 1);
  double cj = sqrt(1.0 * (jmxpos - j1s) * (j1s - jmns));
  int j1mos = j1mo * j1mo;
  double cjmo = sqrt(1.0 * (jmxpos - j1mos) * (j1mos - jmns));
  if (njmo <= 1) {
    c6->rac3j[0] = -cj / cjmo;
    double sj = sjr + (c6->rac3j[0] * c6->rac3j[0]) * (jf - 4);
    double cnr = (1.0 / sqrt(sj)) * isn;
    c6->rac3j[1] = cnr;
    c6->rac3j[0] *= cnr;
    return;
  }
  int nj = njmo + 1;
  int nmat = (nj + 1) / 2;
  c6->rac3j[nj - 1] = 1.0;
  c6->rac3j[njmo - 1] = -cj / cjmo;
  if (nmat != njmo) {
    int nbr = njmo - nmat;
    for (int ibr45 = 1; ibr45 <= nbr; ibr45++) {
      int irr = nj - ibr45;
      jf -= 4;
      j1mo -= 2;
      j1s = (j1mo + 1) * (j1mo + 1);
      cj = sqrt(1.0 * (jmxpos - j1s) * (j1s - jmns));
      j1mos = j1mo * j1mo;
      cjmo = sqrt(1.0 * (jmxpos - j1mos) * (j1mos - jmns));
      c6->rac3j[irr - 2] = c6->rac3j[irr - 1] * (-cj / cjmo);
      sjr = sjr + (c6->rac3j[irr - 1] * c6->rac3j[irr - 1]) * jf;
    }
  }
  // label 50
  double racmat = c6->rac3j[nmat - 1];
  sjr = sjr + (racmat * racmat) * (jf - 4);
  c6->rac3j[0] = 1.0;
  jf = jmn + jmn + 1;
  double sjl = 1.0 * jf;
  int j1pt = jmn + 2;
  int j1pos = (j1pt - 1) * (j1pt - 1);
  double cjpo = sqrt(1.0 * (jmxpos - j1pos) * (j1pos - jmns));
  int j1pts = j1pt * j1pt;
  double cjpt = sqrt(1.0 * (jmxpos - j1pts) * (j1pts - jmns));
  c6->rac3j[1] = -cjpo / cjpt;
  int nmatmo = nmat - 1;
  if (nmatmo >= 2) {
    for (int irl70 = 2; irl70 <= nmatmo; irl70++) {
      jf += 4;
      j1pt += 2;
      j1pos = (j1pt - 1) * (j1pt - 1);
      cjpo = sqrt(1.0 * (jmxpos - j1pos) * (j1pos - jmns));
      j1pts = j1pt * j1pt;
      cjpt = sqrt(1.0 * (jmxpos - j1pts) * (j1pts - jmns));
      c6->rac3j[irl70] = c6->rac3j[irl70 - 1] * (-cjpo / cjpt);
      sjl = sjl + (c6->rac3j[irl70 - 1] * c6->rac3j[irl70 - 1]) * jf;
    }
  }
  // label 75
  double ratrac = racmat / c6->rac3j[nmat - 1];
  double rats = ratrac * ratrac;
  double sj = sjr + sjl * rats;
  c6->rac3j[nmat - 1] = racmat;
  double cnr = (1.0 / sqrt(sj)) * isn;
  for (int irr80 = nmat; irr80 <= nj; irr80++) {
    c6->rac3j[irr80 - 1] *= cnr;
  }
  double cnl = cnr * ratrac;
  for (int irl85 = 1; irl85 <= nmatmo; irl85++) {
    c6->rac3j[irl85 - 1] *= cnl;
  }
}

void r3jjr(int j2, int j3, int m2, int m3, C6 *c6) {
  int jmx = j3 + j2;
  int jdf = j3 - j2;
  int m1 = -m2 - m3;
  int abs_jdf = (jdf >= 0) ? jdf : -jdf;
  int abs_m1 = (m1 >= 0) ? m1 : -m1;
  int jmn = (abs_jdf > abs_m1) ? abs_jdf : abs_m1;
  int njmo = jmx - jmn;
  int jf = jmx + jmx + 1;
  int isn = 1;
  if ((jdf + m1) % 2 != 0) isn = -1;
  if (njmo <= 0) {
    double sj = 1.0 * jf;
    double cnr = (1.0 / sqrt(sj)) * isn;
    c6->rac3j[0] = cnr;
  } else { // label 15
    double sjt = 1.0;
    double sjr = 1.0 * jf;
    int jsmpos = (jmx + 1) * (jmx + 1);
    int jdfs = jdf * jdf;
    int m1s = m1 * m1;
    int mdf = m3 - m2;
    int idjc = m1 * (j3 * (j3 + 1) - j2 * (j2 +1));
    int j1 = jmx;
    int j1s = j1 * j1;
    int j1po = j1 + 1;
    double ccj = 1.0 * (j1s - jdfs) * (j1s - m1s);
    double cj = sqrt(ccj * (jsmpos - j1s));
    double dj = 1.0 * jf * (j1 * j1po * mdf + idjc);
    if (njmo <= 1) {
      c6->rac3j[0] = -dj / (cj * j1po);
      double sj = sjr + (c6->rac3j[0] * c6->rac3j[0]) * (jf - 2);
      double cnr = (1.0 / sqrt(sj)) * isn;
      c6->rac3j[1] = cnr;
      c6->rac3j[0] *= cnr;
    } else { // label 20
      double cjp = 0.0;
      int nj = njmo + 1;
      int nmat = (nj + 1) / 2;
      c6->rac3j[nj - 1] = 1.0;
      c6->rac3j[njmo - 1] = -dj / (cj * j1po);
      if (nmat != njmo) {
	int nbr = njmo - nmat;
	for (int ibr45 = 1; ibr45 <= nbr; ibr45++) {
	  int irr = nj - ibr45;
	  jf -= 2;
	  j1--;
	  j1s = j1 * j1;
	  j1po = j1 + 1;
	  cjp = cj;
	  ccj = 1.0 * (j1s - jdfs) * (j1s - m1s);
	  cj = sqrt(ccj * (jsmpos - j1s));
	  sjt = c6->rac3j[irr - 1] * c6->rac3j[irr - 1];
	  dj = 1.0 * jf * (j1 * j1po * mdf + idjc);
	  c6->rac3j[irr - 2] = -(c6->rac3j[irr - 1] * dj
				 + c6->rac3j[irr] * cjp * j1) / (cj * j1po);
	  sjr += (sjt * jf);
	} // ibr45 loop
      }
      // label 50
      double osjt = sjt;
      sjt = c6->rac3j[nmat - 1] * c6->rac3j[nmat - 1];
      if (sjt >= osjt) {
	sjr += (sjt * (jf - 2));
      } else { // label 55
	nmat++;
      }
      // label 60
      double racmat = c6->rac3j[nmat - 1];
      c6->rac3j[0] = 1.0;
      jf = jmn + jmn + 1;
      double sjl = 1.0 * jf;
      j1 = jmn;
      if (j1 != 0) {
	j1po = j1 + 1;
	int j1pos = j1po * j1po;
	double ccjp = 1.0 * (j1pos - jdfs) * (j1pos - m1s);
	cjp = sqrt(ccjp * (jsmpos - j1pos));
	dj = 1.0 * jf * (j1 * j1po * mdf + idjc);
	c6->rac3j[1] = - dj / (cjp * j1);
      } else { // label 62
	cjp = sqrt(1.0 * (jsmpos - 1));
	dj = 1.0 * mdf;
	c6->rac3j[1] = -dj / cjp;
      }
      // label 63
      int nmatmo = nmat - 1;
      if (nmatmo >= 2) {
	for (int irl70 = 2; irl70 <= nmatmo; irl70++) {
	  jf += 2;
	  j1++;
	  j1po = j1 + 1;
	  int j1pos = j1po * j1po;
	  cj = cjp;
	  double ccjp = 1.0 * (j1pos - jdfs) * (j1pos - m1s);
	  cjp = sqrt(ccjp * (jsmpos - j1pos));
	  sjt = c6->rac3j[irl70 - 1] * c6->rac3j[irl70 - 1];
	  dj = 1.0 * jf * (j1 * j1po * mdf + idjc);
	  c6->rac3j[irl70] = -(
			       c6->rac3j[irl70 - 1] * dj
			       + c6->rac3j[irl70 - 2] * cj * j1po
			       ) / (cjp * j1);
	  sjl += (sjt * jf);
	}
      }
      // label 75
      double ratrac = racmat / c6->rac3j[nmat - 1];
      double rats = ratrac * ratrac;
      double sj = sjr + sjl * rats;
      c6->rac3j[nmat - 1] = racmat;
      double cnr = (1.0 / sqrt(sj)) * isn;
      for (int irr80 = nmat; irr80 <= nj; irr80++) c6->rac3j[irr80 - 1] *= cnr;
      double cnl = cnr * ratrac;
      for (int irl85 = 1; irl85 <= nmatmo; irl85++) c6->rac3j[irl85 - 1] *= cnl;
    }
  }
}

void r3jmr(int j1, int j2, int j3, int m1, C6 *c6) {
  int mmx = (j2 < j3 - m1) ? j2 : j3 - m1;
  int mmn = (-j2 > -(j3 + m1)) ? -j2 : -(j3 + m1);
  int nmmo = mmx - mmn;
  int j1po = j1 + 1;
  int j1tpo = j1po + j1;
  int isn = 1;
  if ((j2 - j3 - m1) % 2 != 0) isn = -1;
  if (nmmo <= 0) {
    double sj = 1.0 * j1tpo;
    double cnr = (1.0 / sqrt(sj)) * isn;
    c6->rac3j[0] = cnr;
    // returns
  } else { // label 15
    int j1s = j1 * j1po;
    int j2po = j2 + 1;
    int j2s = j2 * j2po;
    int j3po = j3 + 1;
    int j3s = j3 * j3po;
    int id = j1s - j2s - j3s;
    int m2 = mmx;
    int m3 = m1 + m2;
    double cm = sqrt(1.0 * (j2po - m2) * (j2 + m2) * (j3po - m3) * (j3 + m3));
    double dm = 1.0 * (id + m2 * m3 * 2);
    if (nmmo <= 1) {
      c6->rac3j[0] = dm / cm;
      double sj = (1.0 + c6->rac3j[0] * c6->rac3j[0]) * j1tpo;
      double cnr = 1.0 / sqrt(sj) * isn;
      c6->rac3j[1] = cnr;
      c6->rac3j[0] *= cnr;
      // returns
    } else { // label 20
      int nm = nmmo + 1;
      int nmat = (nm + 1) / 2;
      c6->rac3j[nm - 1] = 1.0;
      c6->rac3j[nmmo - 1] = dm / cm;
      double sjt = 1.0;
      double sjr = 1.0;
      if (nmat != nmmo) {
	int nbr = nmmo - nmat;
	for (int ibr45 = 1; ibr45 <= nbr; ibr45++) {
	  int irr = nm - ibr45;
	  m2--;
	  m3 = m1 + m2;
	  double cmp = cm;
	  cm = sqrt(1.0 * (j2po - m2) * (j2 + m2) * (j3po - m3) * (j3 + m3));
	  sjt = c6->rac3j[irr - 1] * c6->rac3j[irr - 1];
	  dm = 1.0 * (id + m2 * m3 * 2);
	  c6->rac3j[irr - 1] *= ((dm - c6->rac3j[irr] * cmp) / cm);
	  sjr += sjt;
	} // ibr45 loop
      }
      // label 50
      double osjt = sjt;
      sjt = c6->rac3j[nmat - 1] * c6->rac3j[nmat - 1];
      if (sjt >= osjt) {
	sjr += sjt;
      } else { // label 55
	nmat++;
      }
      // label 60
      double racmat = c6->rac3j[nmat - 1];
      c6->rac3j[0] = 1.0;
      m2 = mmn;
      m3 = m1 + m2;
      double cmp = sqrt(1.0 * (j2 - m2) * (j2po + m2) * (j3 - m3) * (j3po + m3));
      dm = 1.0 * (id + m2 * m3 * 2);
      c6->rac3j[1] = dm / cmp;
      double sjl = 1.0;
      int nmatmo = nmat - 1;
      if (nmatmo > 1) {
	for (int irl70 = 2; irl70 <= nmatmo; irl70++) {
	  m2++;
	  m3 = m1 + m2;
	  cm = cmp;
	  cmp = sqrt(1.0 * (j2 - m2) * (j2po + m2) * (j3 - m3) * (j3po + m3));
	  sjt = c6->rac3j[irl70 - 1] * c6->rac3j[irl70 - 1];
	  dm = 1.0 * (id + m2 * m3 * 2);
	  c6->rac3j[irl70] = (c6->rac3j[irl70 - 1] * dm - c6->rac3j[irl70 - 2] * cm) / cmp;
	  sjl += sjt;
	}
      } // label 75
      double ratrac = racmat / c6->rac3j[nmat - 1];
      double rats = ratrac * ratrac;
      double sj = (sjr + sjl * rats) * j1tpo;
      c6->rac3j[nmat - 1] = racmat;
      double cnr = 1.0 / sqrt(sj) * isn;
      for (int irr80 = nmat; irr80 <= nm; irr80++) c6->rac3j[irr80 - 1] *= cnr;
      double cnl = cnr * ratrac;
      for (int irl85 = 1; irl85 <= nmatmo; irl85++) c6->rac3j[irl85 - 1] *= cnl;
      // returns
    }
  }
}

void raba(
	  int le, complex<double> **am0m, complex<double> **w, double **tqce,
	  complex<double> **tqcpe, double **tqcs, complex<double> **tqcps
) {
  complex<double> **a, **ctqce, **ctqcs;
  complex<double> acw, acwp, aca, acap, c1, c2, c3;
  const complex<double> cc0(0.0, 0.0);
  const complex<double> uim(0.0, 1.0);
  const double sq2i = 1.0 / sqrt(2.0);
  int nlem = le * (le + 2);
  const int nlemt = nlem + nlem;
  a = new complex<double>*[nlemt];
  ctqce = new complex<double>*[2];
  ctqcs = new complex<double>*[2];
  for (int ai = 0; ai < nlemt; ai++) a[ai] = new complex<double>[2]();
  for (int ci = 0; ci < 2; ci++) {
    ctqce[ci] = new complex<double>[3]();
    ctqcs[ci] = new complex<double>[3]();
  }
  for (int i20 = 1; i20 <= nlemt; i20++) {
    int i = i20 - 1;
    c1 = cc0;
    c2 = cc0;
    for (int j10 = 1; j10 <= nlemt; j10++) {
      int j = j10 - 1;
      c1 += (am0m[i][j] * w[j][0]);
      c2 += (am0m[i][j] * w[j][1]);
    } // j10 loop
    a[i][0] = c1;
    a[i][1] = c2;
  } //i20 loop
  int jpo = 2;
  for (int ipo70 = 1; ipo70 <= 2; ipo70++) {
    if (ipo70 == 2) jpo = 1;
    int ipo = ipo70 - 1;
    ctqce[ipo][0] = cc0;
    ctqce[ipo][1] = cc0;
    ctqce[ipo][2] = cc0;
    tqcpe[ipo][0] = cc0;
    tqcpe[ipo][1] = cc0;
    tqcpe[ipo][2] = cc0;
    ctqcs[ipo][0] = cc0;
    ctqcs[ipo][1] = cc0;
    ctqcs[ipo][2] = cc0;
    tqcps[ipo][0] = cc0;
    tqcps[ipo][1] = cc0;
    tqcps[ipo][2] = cc0;
    for (int l60 = 1; l60 <= le; l60 ++) {
      int lpo = l60 + 1;
      int il = l60 * lpo;
      int ltpo = l60 + lpo;
      for (int im60 = 1; im60 <= ltpo; im60++) {
	int m = im60 - lpo;
	int i = m + il;
	int ie = i + nlem;
	int mmmu = m + 1;
	int mmmmu = (mmmu > 0) ? mmmu : -mmmu;
	double  rmu = 0.0;
	if (mmmmu <= l60) {
	  int immu = mmmu + il;
	  int immue = immu + nlem;
	  rmu = -sqrt(1.0 * (l60 + mmmu) * (l60 - m)) * sq2i;
	  acw = dconjg(a[i - 1][ipo]) * w[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * w[immue - 1][ipo];
	  acwp = dconjg(a[i - 1][ipo]) * w[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * w[immue - 1][jpo - 1];
	  aca = dconjg(a[i - 1][ipo]) * a[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * a[immue - 1][ipo];
	  acap = dconjg(a[i - 1][ipo]) * a[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * a[immue - 1][jpo - 1];
	  ctqce[ipo][0] += (acw * rmu);
	  tqcpe[ipo][0] += (acwp * rmu);
	  ctqcs[ipo][0] += (aca * rmu);
	  tqcps[ipo][0] += (acap * rmu);
	}
	// label 30
	rmu = -1.0 * m;
	acw = dconjg(a[i - 1][ipo]) * w[i - 1][ipo] + dconjg(a[ie - 1][ipo]) * w[ie - 1][ipo];
	acwp = dconjg(a[i - 1][ipo]) * w[i - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * w[ie - 1][jpo - 1];
	aca = dconjg(a[i - 1][ipo]) * a[i - 1][ipo] + dconjg(a[ie - 1][ipo]) * a[ie - 1][ipo];
	acap = dconjg(a[i - 1][ipo]) * a[i - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * a[ie - 1][jpo - 1];
	ctqce[ipo][1] += (acw * rmu);
	tqcpe[ipo][1] += (acwp * rmu);
	ctqcs[ipo][1] += (aca * rmu);
	tqcps[ipo][1] += (acap * rmu);
	mmmu = m - 1;
	mmmmu = (mmmu > 0) ? mmmu : -mmmu;
	if (mmmmu <= l60) {
	  int immu = mmmu + il;
	  int immue = immu + nlem;
	  rmu = sqrt(1.0 * (l60 - mmmu) * (l60 + m)) * sq2i;
	  acw = dconjg(a[i - 1][ipo]) * w[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * w[immue - 1][ipo];
	  acwp = dconjg(a[i - 1][ipo]) * w[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * w[immue - 1][jpo - 1];
	  aca = dconjg(a[i - 1][ipo]) * a[immu - 1][ipo] + dconjg(a[ie - 1][ipo]) * a[immue - 1][ipo];
	  acap = dconjg(a[i - 1][ipo]) * a[immu - 1][jpo - 1] + dconjg(a[ie - 1][ipo]) * a[immue - 1][jpo - 1];
	  ctqce[ipo][2] += (acw * rmu);
	  tqcpe[ipo][2] += (acwp * rmu);
	  ctqcs[ipo][2] += (aca * rmu);
	  tqcps[ipo][2] += (acap * rmu);
	} // ends im60 loop
      } // im60 loop
    } // l60 loop
  } // ipo70 loop
  for (int ipo78 = 1; ipo78 <= 2; ipo78++) {
    int ipo = ipo78 - 1;
    tqce[ipo][0] = (ctqce[ipo][0] - ctqce[ipo][2]).real() * sq2i;
    tqce[ipo][1] = ((ctqce[ipo][0] + ctqce[ipo][2]) * uim).real() * sq2i;
    tqce[ipo][2] = ctqce[ipo][1].real();
    c1 = tqcpe[ipo][0];
    c2 = tqcpe[ipo][1];
    c3 = tqcpe[ipo][2];
    tqcpe[ipo][0] = (c1 - c3) * sq2i;
    tqcpe[ipo][1] = (c1 + c3) * (uim * sq2i);
    tqcpe[ipo][2] = c2;
    tqcs[ipo][0] = -sq2i * (ctqcs[ipo][0] - ctqcs[ipo][2]).real();
    tqcs[ipo][1] = -sq2i * ((ctqcs[ipo][0] + ctqcs[ipo][2]) * uim).real();
    tqcs[ipo][2] = -1.0 * ctqcs[ipo][1].real();
    c1 = tqcps[ipo][0];
    c2 = tqcps[ipo][1];
    c3 = tqcps[ipo][2];
    tqcps[ipo][0] = -(c1 - c3) * sq2i;
    tqcps[ipo][1] = -(c1 + c3) * (uim * sq2i);
    tqcps[ipo][2] = -c2;
  } // ipo78 loop
  // Clean memory
  for (int ai = 0; ai < nlemt; ai++) delete[] a[ai];
  for (int ci = 0; ci < 2; ci++) {
    delete[] ctqce[ci];
    delete[] ctqcs[ci];
  }
  delete[] a;
  delete[] ctqce;
  delete[] ctqcs;
}

void rftr(
	  double *u, double *up, double *un, double *gapv, double extins, double scatts,
	  double &rapr, double &cosav, double &fp, double &fn, double &fk, double &fx,
	  double &fy, double &fz
) {
  fk = u[0] * gapv[0] + u[1] * gapv[1] + u[2] * gapv[2];
  rapr = extins - fk;
  cosav = fk / scatts;
  fp = -(up[0] * gapv[0] + up[1] * gapv[1] + up[2] * gapv[2]);
  fn = -(un[0] * gapv[0] + un[1] * gapv[1] + un[2] * gapv[2]);
  fk = rapr;
  fx = u[0] * extins - gapv[0];
  fy = u[1] * extins - gapv[1];
  fz = u[2] * extins - gapv[2];
}

void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
  const complex<double> cc0(0.0, 0.0);
  double exdc = exri * exri;
  double ccs = 4.0 * acos(0.0) / (vk * vk);
  double cccs = ccs / exdc;
  complex<double> sum21, rm, re, csam;
  csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
  //double scs = 0.0, ecs = 0.0, acs = 0.0;
  c3->scs = 0.0;
  c3->ecs = 0.0;
  c3->acs = 0.0;
  c3->tfsas = cc0;
  for (int i14 = 1; i14 <= c4->nsph; i14++) {
    int iogi = c1->iog[i14 - 1];
    if (iogi >= i14) {
      double sums = 0.0;
      sum21 = cc0;
      for (int l10 = 1; l10 <= c4->li; l10++) {
	double fl = 1.0 * (l10 + l10 + 1);
	rm = 1.0 / c1->rmi[l10 - 1][i14 - 1];
	re = 1.0 / c1->rei[l10 - 1][i14 - 1];
	double rvalue = (dconjg(rm) * rm + dconjg(re) * re).real() * fl;
	sums += rvalue;
	sum21 += ((rm + re) * fl);
      } // l10 loop
      sum21 *= -1.0;
      double scasec = cccs * sums;
      double extsec = -cccs * sum21.real();
      double abssec = extsec - scasec;
      c1->sscs[i14 - 1] = scasec;
      c1->sexs[i14 - 1] = extsec;
      c1->sabs[i14 - 1] = abssec;
      double gcss = c1->gcsv[i14 - 1];
      c1->sqscs[i14 - 1] = scasec / gcss;
      c1->sqexs[i14 - 1] = extsec / gcss;
      c1->sqabs[i14 - 1] = abssec / gcss;
      c1->fsas[i14 - 1] = sum21 * csam;
    }
    // label 12
    c3->scs += c1->sscs[iogi - 1];
    c3->ecs += c1->sexs[iogi - 1];
    c3->acs += c1->sabs[iogi - 1];
    c3->tfsas += c1->fsas[iogi - 1];
  } // i14 loop
}

void scr2(
	  double vk, double vkarg, double exri, double *duk, C1 *c1, C1_AddOns *c1ao,
	  C3 *c3, C4 *c4
) {
  const complex<double> cc0(0.0, 0.0);
  const complex<double> uim(0.0, 1.0);
  complex<double> s11, s21, s12, s22, rm, re, csam, cph, phas, cc;
  double ccs = 1.0 / (vk * vk);
  csam = -(ccs / (exri * vk)) * complex<double>(0.0, 0.5);
  const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
  double cfsq = 4.0 / (pi4sq * ccs * ccs);
  cph = uim * exri * vkarg;
  int ls = (c4->li < c4->le) ? c4->li : c4->le;
  c3->tsas[0][0] = cc0;
  c3->tsas[1][0] = cc0;
  c3->tsas[0][1] = cc0;
  c3->tsas[1][1] = cc0;
  for (int i14 = 1; i14 <= c4->nsph; i14++) {
    int i = i14 - 1;
    int iogi = c1->iog[i14 - 1];
    if (iogi >= i14) {
      int k = 0;
      s11 = cc0;
      s21 = cc0;
      s12 = cc0;
      s22 = cc0;
      for (int l10 = 1; l10 <= ls; l10++) {
	int l = l10 - 1;
	rm = 1.0 / c1->rmi[l][i];
	re = 1.0 / c1->rei[l][i];
	int ltpo = l10 + l10 + 1;
	for (int im10 = 1; im10 <= ltpo; im10++) {
	  k++;
	  int ke = k + c4->nlem;
	  s11 -= (c1->w[k - 1][2] * c1->w[k - 1][0] * rm + c1->w[ke - 1][2] * c1->w[ke - 1][0] * re);
	  s21 -= (c1->w[k - 1][3] * c1->w[k - 1][0] * rm + c1->w[ke - 1][3] * c1->w[ke - 1][0] * re);
	  s12 -= (c1->w[k - 1][2] * c1->w[k - 1][1] * rm + c1->w[ke - 1][2] * c1->w[ke - 1][1] * re);
	  s22 -= (c1->w[k - 1][3] * c1->w[k - 1][1] * rm + c1->w[ke - 1][3] * c1->w[ke - 1][1] * re);
	} // im10 loop
      } // l10 loop
      c1->sas[i][0][0] = s11 * csam;
      c1->sas[i][1][0] = s21 * csam;
      c1->sas[i][0][1] = s12 * csam;
      c1->sas[i][1][1] = s22 * csam;
    }
    // label 12
    phas = exp(cph * (duk[0] * c1->rxx[i] + duk[1] * c1->ryy[i] + duk[2] * c1->rzz[i]));
    c3->tsas[0][0] += (c1->sas[iogi - 1][0][0] * phas);
    c3->tsas[1][0] += (c1->sas[iogi - 1][1][0] * phas);
    c3->tsas[0][1] += (c1->sas[iogi - 1][0][1] * phas);
    c3->tsas[1][1] += (c1->sas[iogi - 1][1][1] * phas);
  } // i14 loop
  for (int i24 = 1; i24 <= c4->nsph; i24++) {
    int iogi = c1->iog[i24 - 1];
    if (iogi >= i24) {
      int j = 0;
      for (int ipo1 = 1; ipo1 <=2; ipo1++) {
	for (int jpo1 = 1; jpo1 <= 2; jpo1++) {
	  cc = dconjg(c1->sas[i24 - 1][jpo1 - 1][ipo1 - 1]);
	  for (int ipo2 = 1; ipo2 <= 2; ipo2++) {
	    for (int jpo2 = 1; jpo2 <= 2; jpo2++) {
	      j++;
	      c1ao->vints[i24 - 1][j - 1] = c1->sas[i24 - 1][jpo2 - 1][ipo2 - 1] * cc * cfsq;
	    } // jpo2 loop
	  } // ipo2 loop
	} // jpo1 loop
      } // ipo1 loop
    }
  } // i24 loop
  int j = 0;
  for (int ipo1 = 1; ipo1 <=2; ipo1++) {
    for (int jpo1 = 1; jpo1 <= 2; jpo1++) {
      cc = dconjg(c3->tsas[jpo1 - 1][ipo1 - 1]);
      for (int ipo2 = 1; ipo2 <= 2; ipo2++) {
	for (int jpo2 = 1; jpo2 <= 2; jpo2++) {
	  j++;
	  c1ao->vintt[j - 1] = c3->tsas[jpo2 - 1][ipo2 - 1] * cc * cfsq;
	} // jpo2 loop
      } // ipo2 loop
    } // jpo1 loop
  } // ipo1 loop
}

void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
  complex<double> *ylm;
  const double pi = acos(-1.0);
  c3->gcs = 0.0;
  double gcss = 0.0;
  for (int i18 = 1; i18 <= c4->nsph; i18++) {
    int iogi = c1->iog[i18 - 1];
    if (iogi >= i18) {
      gcss = pi * c1->ros[i18 - 1] * c1->ros[i18 - 1];
      c1->gcsv[i18 - 1] = gcss;
      int nsh = c1->nshl[i18 - 1];
      for (int j16 = 1; j16 <= nsh; j16++) {
	c1->rc[i18 - 1][j16 - 1] = rcf[i18 - 1][j16 - 1] * c1->ros[i18 - 1];
      } // j16 loop
    }
    c3->gcs += gcss;
  } // i18 loop
  int ylm_size = (c4->litpos > c4->lmtpos) ? c4->litpos : c4->lmtpos;
  ylm = new complex<double>[ylm_size]();
  int i = 0;
  for (int l1po28 = 1; l1po28 <= c4->lmpo; l1po28++) {
    int l1 = l1po28 - 1;
    for (int l2 = 1; l2 <= c4->lm; l2++) {
      r3j000(l1, l2, c6);
      c1ao->ind3j[l1po28 - 1][l2 - 1] = i;
      int lmnpo = (l2 > l1) ? l2 - l1 + 1 : l1 - l2 + 1;
      int lmxpo = l2 + l1 + 1;
      int lpo28 = lmnpo;
      int il = 0;
      while (lpo28 <= lmxpo) {
	i++;
	il++;
	c1ao->v3j0[i - 1] = c6->rac3j[il - 1];
	lpo28 += 2;
      }
    } // l2 loop
  } // l1po28 loop
  int nsphmo = c4->nsph - 1;
  int lit = c4->li + c4->li;
  int ivy = 0;
  for (int nf40 = 1; nf40 <= nsphmo; nf40++) { // GPU portable?
    int nfpo = nf40 + 1;
    for (int ns40 = nfpo; ns40 <= c4->nsph; ns40++) {
      double rx = c1->rxx[nf40 - 1] - c1->rxx[ns40 - 1];
      double ry = c1->ryy[nf40 - 1] - c1->ryy[ns40 - 1];
      double rz = c1->rzz[nf40 - 1] - c1->rzz[ns40 - 1];
      double rr = 0.0;
      double crth = 0.0, srth = 0.0, crph = 0.0, srph = 0.0;
      polar(rx, ry, rz, rr, crth, srth, crph, srph);
      sphar(crth, srth, crph, srph, lit, ylm);
      for (int iv38 = 1; iv38 <= c4->litpos; iv38++) {
	c1ao->vyhj[iv38 + ivy - 1] = dconjg(ylm[iv38 - 1]);
      } // iv38 loop
      ivy += c4->litpos;
    } // ns40 loop
  } // nf40 loop
  int lmt = c4->li + c4->le;
  ivy = 0;
  for (int nf50 = 1; nf50 <= c4->nsph; nf50++) {
    double rx = c1->rxx[nf50 - 1];
    double ry = c1->ryy[nf50 - 1];
    double rz = c1->rzz[nf50 - 1];
    if (rx != 0.0 || ry != 0.0 || rz != 0.0) {
      double rr = 0.0;
      double crth = 0.0, srth = 0.0, crph = 0.0, srph = 0.0;
      polar(rx, ry, rz, rr, crth, srth, crph, srph);
      sphar(crth, srth, crph, srph, lmt, ylm);
      for (int iv48 = 1; iv48 <= c4->lmtpos; iv48++) {
	c1ao->vyj0[iv48 + ivy - 1] = dconjg(ylm[iv48 - 1]);
      } // iv48 loop
    }
    ivy += c4->lmtpos;
  } // nf50 loop
  delete[] ylm;
}

void tqr(
	 double *u, double *up, double *un, double *tqev, double *tqsv, double &tep,
	 double &ten, double &tek, double &tsp, double &tsn, double &tsk
) {
  tep = up[0] * tqev[0] + up[1] * tqev[1] + up[2] * tqev[2];
  ten = un[0] * tqev[0] + un[1] * tqev[1] + un[2] * tqev[2];
  tek = u[0] * tqev[0] + u[1] * tqev[1] + u[2] * tqev[2];
  tsp = up[0] * tqsv[0] + up[1] * tqsv[1] + up[2] * tqsv[2];
  tsn = un[0] * tqsv[0] + un[1] * tqsv[1] + un[2] * tqsv[2];
  tsk = u[0] * tqsv[0] + u[1] * tqsv[1] + u[2] * tqsv[2];
}

void ztm(complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9) {
  complex<double> gie, gle, a1, a2, a3, a4, sum1, sum2, sum3, sum4;
  const complex<double> cc0(0.0, 0.0);
  int ndi = c4->nsph * c4->nlim;
  int i2 = 0;
  for (int n2 = 1; n2 <= c4->nsph; n2++) { // GPU portable?
    for (int l2 = 1; l2 <= c4->li; l2++) {
      int l2tpo = l2 + l2 + 1;
      int m2 = -l2 - 1;
      for (int im2 = 1; im2 <= l2tpo; im2++) {
	m2++;
	i2++;
	int i3 = 0;
	for (int l3 = 1; l3 <= c4->le; l3++) {
	  int l3tpo = l3 + l3 + 1;
	  int m3 = -l3 - 1;
	  for (int im3 = 1; im3 <= l3tpo; im3++) {
	    m3++;
	    i3++;
	    c9->gis[i2 - 1][i3 - 1] = ghit(2, 0, n2, l2, m2, l3, m3, c1, c1ao, c4, c6);
	    c9->gls[i2 - 1][i3 - 1] = ghit(2, 1, n2, l2, m2, l3, m3, c1, c1ao, c4, c6);
	  } // im3 loop
	} // l3 loop
      } // im2 loop
    } // l2 loop
  } // n2 loop
  for (int i1 = 1; i1 <= ndi; i1++) { // GPU portable?
    int i1e = i1 + ndi;
    for (int i3 = 1; i3 <= c4->nlem; i3++) {
      int i3e = i3 + c4->nlem;
      sum1 = cc0;
      sum2 = cc0;
      sum3 = cc0;
      sum4 = cc0;
      for (int i2 = 1; i2 <= ndi; i2++) {
	int i2e = i2 + ndi;
	gie = c9->gis[i2 - 1][i3 - 1];
	gle = c9->gls[i2 - 1][i3 - 1];
	a1 = am[i1 - 1][i2 - 1];
	a2 = am[i1 - 1][i2e - 1];
	a3 = am[i1e - 1][i2 - 1];
	a4 = am[i1e - 1][i2e - 1];
	sum1 += (a1 * gie + a2 * gle);
	sum2 += (a1 * gle + a2 * gie);
	sum3 += (a3 * gie + a4 * gle);
	sum4 += (a3 * gle + a4 * gie);
      } // i2 loop
      c9->sam[i1 - 1][i3 - 1] = sum1;
      c9->sam[i1 - 1][i3e - 1] = sum2;
      c9->sam[i1e - 1][i3 - 1] = sum3;
      c9->sam[i1e - 1][i3e - 1] = sum4;
    } // i3 loop
  } // i1 loop
  for (int i1 = 1; i1 <= ndi; i1++) {
    for (int i0 = 1; i0 <= c4->nlem; i0++) {
      c9->gis[i1 - 1][i0 - 1] = dconjg(c9->gis[i1 - 1][i0 - 1]);
      c9->gls[i1 - 1][i0 - 1] = dconjg(c9->gls[i1 - 1][i0 - 1]);
    } // i0 loop
  } // i1 loop
  int nlemt = c4->nlem + c4->nlem;
  for (int i0 = 1; i0 <= c4->nlem; i0++) {
    int i0e = i0 + c4->nlem;
    for (int i3 = 1; i3 <= nlemt; i3++) {
      sum1 = cc0;
      sum2 = cc0;
      for (int i1 = 1; i1 <= ndi; i1 ++) {
	int i1e = i1 + ndi;
	a1 = c9->sam[i1 - 1][i3 - 1];
	a2 = c9->sam[i1e - 1][i3 - 1];
	gie = c9->gis[i1 - 1][i0 - 1];
	gle = c9->gls[i1 - 1][i0 - 1];
	sum1 += (a1 * gie + a2 * gle);
	sum2 += (a1 * gle + a2 * gie);
      } // i1 loop
      c1ao->am0m[i0 - 1][i3 - 1] = -sum1;
      c1ao->am0m[i0e - 1][i3 - 1] = -sum2;
    } // i3 loop
  } // i0 loop
}
