diff --git a/src/cluster/Makefile b/src/cluster/Makefile
index 71a106ed2ecb24c6467f168f19506133457bdd18..181647521ebfdeadc696de00957a8d471153e6c3 100644
--- a/src/cluster/Makefile
+++ b/src/cluster/Makefile
@@ -15,8 +15,8 @@ clu: clu.o
 edfb: edfb.o
 	$(FC) $(FCFLAGS) -o $(BUILDDIR)/edfb $(BUILDDIR)/edfb.o $(LFLAGS)
 
-np_cluster: $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/cluster.o
-	$(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_cluster $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sphere.o $(BUILDDIR)/cluster.o
+np_cluster: $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/clu_subs.o $(BUILDDIR)/cluster.o
+	$(CXX) $(CXXFLAGS) $(CXXLFLAGS) -o $(BUILDDIR)/np_cluster $(BUILDDIR)/np_cluster.o $(BUILDDIR)/Commons.o $(BUILDDIR)/Configuration.o $(BUILDDIR)/Parsers.o $(BUILDDIR)/sph_subs.o $(BUILDDIR)/clu_subs.o $(BUILDDIR)/cluster.o
 
 $(BUILDDIR)/np_cluster.o:
 	$(CXX) $(CXXFLAGS) -c np_cluster.cpp -o $(BUILDDIR)/np_cluster.o
@@ -33,8 +33,11 @@ $(BUILDDIR)/Parsers.o:
 $(BUILDDIR)/cluster.o:
 	$(CXX) $(CXXFLAGS) -c cluster.cpp -o $(BUILDDIR)/cluster.o
 
+$(BUILDDIR)/clu_subs.o:
+	$(CXX) $(CXXFLAGS) -c ../libnptm/clu_subs.cpp -o $(BUILDDIR)/clu_subs.o
+
 $(BUILDDIR)/sph_subs.o:
-	$(CXX) $(CXXFLAGS) -c ../libnptm/sph_subs.cpp -o $(BUILDDIR)/sphere.o
+	$(CXX) $(CXXFLAGS) -c ../libnptm/sph_subs.cpp -o $(BUILDDIR)/sph_subs.o
 
 clean:
 	rm -f $(BUILDDIR)/*.o
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 9e13b5267940c34b75059b578a497f93372b9408..af6ece2314a68490e3db244ac90221867f30c98d 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -2,13 +2,15 @@
  *
  * \brief C++ porting of CLU functions and subroutines.
  *
- * Remember that FORTRAN passes arguments by reference, so, every time we use
- * a subroutine call, we need to add a referencing layer to the C++ variable.
- * All the functions defined below need to be properly documented and ported
- * to C++.
  *
- * Currently, only basic documenting information about functions and parameter
- * types are given, to avoid doxygen warning messages.
+ * 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`
+ * 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_COMMONS_H_
@@ -27,11 +29,11 @@ extern double cg1(int lmpml, int mu, int l, int m);
 extern void dme(
 		int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
 		C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg
-		);
+);
 extern void rabas(
 		  int inpol, int li, int nsph, C1 *c1, double **tqse, std::complex<double> **tqspe,
 		  double **tqss, std::complex<double> **tqsps
-		  );
+);
 extern void rbf(int n, double x, int &nm, double sj[]);
 extern void rnf(int n, double x, int &nm, double sy[]);
 extern void mmulc(std::complex<double> *vint, double **cmullr, double **cmul);
@@ -40,552 +42,25 @@ extern void thdps(int lm, double ****zpv);
 extern void upvmp(
 		  double thd, double phd, int icspnv, double &cost, double &sint,
 		  double &cosp, double &sinp, double *u, double *up, double *un
-		  );
+);
 extern void upvsp(
 		  double *u, double *upmp, double *unmp, double *us, double *upsmp, double *unsmp,
 		  double *up, double *un, double *ups, double *uns, double *duk, int &isq,
 		  int &ibf, double &scand, double &cfmp, double &sfmp, double &cfsp, double &sfsp
-		  );
+);
 extern void wmamp(
 		  int iis, double cost, double sint, double cosp, double sinp, int inpol,
 		  int lm, int idot, int nsph, double *arg, double *u, double *up,
 		  double *un, C1 *c1
-		  );
+);
 extern void wmasp(
 		  double cost, double sint, double cosp, double sinp, double costs, double sints,
 		  double cosps, double sinps, double *u, double *up, double *un, double *us,
 		  double *ups, double *uns, int isq, int ibf, int inpol, int lm, int idot,
 		  int nsph, double *argi, double *args, C1 *c1
-		  );
+);
 // >>> END OF SPH_SUBS DECLARATION <<<
 
-/*! \brief C++ porting of CDTP
- *
- * \param z: `complex<double>`
- * \param am: Matrix of complex.
- * \param i: `int`
- * \param jf: `int`
- * \param k: `int`
- * \param nj: `int`
- */
-std::complex<double> cdtp(
-			  std::complex<double> z, std::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.
-   */
-  std::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;
-}
-
-/*! \brief C++ porting of CGEV
- *
- * \param ipamo: `int`
- * \param mu: `int`
- * \param l: `int`
- * \param m: `int`
- * \return result: `double`
- */
-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;
-}
-
-/*! \brief C++ porting of R3JJR
- *
- * \param j2: `int`
- * \param j3: `int`
- * \param m2: `int`
- * \param m3: `int`
- * \param c6: `C6 *`
- */
-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);
-    // In old version, CJP was defined here. Did not work.
-    // double cjp = 0.0
-    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;
-    }
-  }
-}
-
-/*! \brief C++ porting of R3JMR
- *
- * \param j1: `int`
- * \param j2: `int`
- * \param j3: `int`
- * \param m1: `int`
- * \param c6: `C6 *`
- */
-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
-    }
-  }
-}
-
-/*! \brief C++ porting of GHIT
- *
- * \param ihi: `int`
- * \param ipamo: `int`
- * \param nbl: `int`
- * \param l1: `int`
- * \param m1: `int`
- * \param l2: `int`
- * \param m2: `int`
- * \param c1: `C1 *`
- * \param c1ao: `C1_AddOns *`
- * \param c4: `C4 *`
- * \param c6: `C6 *`
- */
-std::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 std::complex<double> cc0(0.0, 0.0);
-  const std::complex<double> uim(0.0, 1.0);
-  std::complex<double> csum(0.0, 0.0), cfun(0.0, 0.0);
-  std::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 = std::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;
-		//printf("DEBUG: VH( %d ) = (%lE, %lE)\n", (nbhj + lt20), c1ao->vh[nbhj + lt20 - 1].real(), c1ao->vh[nbhj + lt20 - 1].imag());
-		//printf("DEBUG: VYHJ( %d ) = (%lE, %lE)\n", (nby + ny), c1ao->vyhj[nby + ny - 1].real(), c1ao->vyhj[nby + ny - 1].imag());
-		cfun = (c1ao->vh[nbhj + lt20 - 1] * c1ao->vyhj[nby + ny - 1]) * f3j;
-		csum += cfun; // we were here
-	      }
-	      // 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;
-}
-
 /*! \brief C++ porting of APC
  *
  * \param zpv: `double ****`
@@ -599,122 +74,7 @@ std::complex<double> ghit(
 void apc(
 	 double ****zpv, int le, std::complex<double> **am0m, std::complex<double> **w,
 	 double sqk, double **gapr, std::complex<double> **gapp
-	 ) {
-  std::complex<double> **ac, **gap;
-  const std::complex<double> cc0(0.0, 0.0);
-  const std::complex<double> uim(0.0, 1.0);
-  std::complex<double> uimmp, summ, sume, suem, suee, summp, sumep;
-  std::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 std::complex<double>*[nlemt];
-  gap = new std::complex<double>*[3];
-  for (int ai = 0; ai < nlemt; ai++) ac[ai] = new std::complex<double>[2]();
-  for (int gi = 0; gi < 3; gi++) gap[gi] = new std::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;
-	      //printf("DEBUG: i=%d, ipo=%d, imp=%d\n", i, ipo, imp);
-	      //fflush(stdout);
-	      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;
-}
+);
 
 /*! \brief C++ porting of APCRA
  *
@@ -729,222 +89,31 @@ void apc(
 void apcra(
 	   double ****zpv, const int le, std::complex<double> **am0m, int inpol, double sqk,
 	   double **gaprm, std::complex<double> **gappm
-	   ) {
-  const std::complex<double> cc0(0.0, 0.0);
-  const std::complex<double> uim(0.0, 1.0);
-  std::complex<double> uimtl, uimtls, ca11, ca12, ca21, ca22;
-  std::complex<double> a11, a12, a21, a22, sum1, sum2, fc;
-  double ****svw = new double***[le];
-  std::complex<double> ****svs = new std::complex<double>***[le];
-  for (int i = 0; i < le; i++) {
-    svw[i] = new double**[3];
-    svs[i] = new std::complex<double>**[3];
-    for (int j = 0; j < 3; j++) {
-      svw[i][j] = new double*[2];
-      svs[i][j] = new std::complex<double>*[2];
-      for (int k = 0; k < 2; k++) {
-	svw[i][j][k] = new double[2]();
-	svs[i][j][k] = new std::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 = std::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 = std::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;
-}
+);
+
+/*! \brief C++ porting of CDTP
+ *
+ * \param z: `complex<double>`
+ * \param am: Matrix of complex.
+ * \param i: `int`
+ * \param jf: `int`
+ * \param k: `int`
+ * \param nj: `int`
+ */
+std::complex<double> cdtp(
+			  std::complex<double> z, std::complex<double> **am, int i, int jf,
+			  int k, int nj
+);
+
+/*! \brief C++ porting of CGEV
+ *
+ * \param ipamo: `int`
+ * \param mu: `int`
+ * \param l: `int`
+ * \param m: `int`
+ * \return result: `double`
+ */
+double cgev(int ipamo, int mu, int l, int m);
 
 /*! \brief C++ porting of CMS
  *
@@ -954,87 +123,7 @@ void apcra(
  * \param c4: `C4 *`
  * \param c6: `C6 *`
  */
-void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
-  std::complex<double> dm, de, cgh, cgk;
-  const std::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 cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
 
 /*! \brief C++ porting of CRSM1
  *
@@ -1045,125 +134,26 @@ void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
  * \param c4: `C4 *`
  * \param c6: `C6 *`
  */
-void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
-  std::complex<double> ***svf, ***svw, **svs;
-  const std::complex<double> cc0(0.0, 0.0);
-  std::complex<double> cam(0.0, 0.0);
-  const int le4po = 4 * c4->le + 1;
-  svf = new std::complex<double>**[le4po];
-  svw = new std::complex<double>**[le4po];
-  svs = new std::complex<double>*[le4po];
-  for (int si = 0; si < le4po; si++) {
-    svf[si] = new std::complex<double>*[le4po];
-    svw[si] = new std::complex<double>*[4];
-    svs[si] = new std::complex<double>[4]();
-    for (int sj = 0; sj < le4po; sj++) svf[si][sj] = new std::complex<double>[4]();
-    for (int sj = 0; sj < 4; sj++) svw[si][sj] = new std::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;
+void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
 
-  // 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;
-}
+/*! \brief C++ porting of GHIT
+ *
+ * \param ihi: `int`
+ * \param ipamo: `int`
+ * \param nbl: `int`
+ * \param l1: `int`
+ * \param m1: `int`
+ * \param l2: `int`
+ * \param m2: `int`
+ * \param c1: `C1 *`
+ * \param c1ao: `C1_AddOns *`
+ * \param c4: `C4 *`
+ * \param c6: `C6 *`
+ */
+std::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
+);
 
 /*! \brief C++ porting of HJV
  *
@@ -1179,72 +169,7 @@ void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
 void hjv(
 	 double exri, double vk, int &jer, int &lcalc, std::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 = std::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] = std::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;
-}
+);
 
 /*! \brief C++ porting of LUCIN
  *
@@ -1253,110 +178,7 @@ void hjv(
  * \param n: `int`
  * \param ier: `int &`
  */
-void lucin(std::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];
-  std::complex<double> ctemp, cfun;
-  std::complex<double> cc0 = std::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 lucin(std::complex<double> **am, const int nddmst, int n, int &ier);
 
 /*! \brief C++ porting of MEXTC
  *
@@ -1366,55 +188,7 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
  * \param cextlr: `double **`
  * \param cext: `double **`
  */
-void mextc(double vk, double exri, std::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 mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr, double **cext);
 
 /*! \brief C++ porting of PCROS
  *
@@ -1425,64 +199,7 @@ void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr,
  * \param c1ao: `C1_AddOns *`
  * \param c4: `C4 *`
  */
-void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
-  const std::complex<double> cc0(0.0, 0.0);
-  std::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)) * std::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 pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4);
 
 /*! \brief C++ porting of PCRSM0
  *
@@ -1493,71 +210,7 @@ void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
  * \param c1ao: `C1_AddOns *`
  * \param c4: `C4 *`
  */
-void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4) {
-  const std::complex<double> cc0(0.0, 0.0);
-  const std::complex<double> uim(0.0, 1.0);
-  std::complex<double> sum1, sum2, sum3, sum4, sumpd;
-  std::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)) * std::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 pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4);
 
 /*! \brief C++ porting of POLAR
  *
@@ -1573,62 +226,7 @@ void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4)
 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
-    }
-  }
-}
+);
 
 /*! \brief C++ porting of R3J000
  *
@@ -1636,97 +234,27 @@ void polar(
  * \param j3: `int`
  * \param c6: `C6 *` Pointer to a C6 instance.
  */
-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 r3j000(int j2, int j3, C6 *c6);
+
+/*! \brief C++ porting of R3JJR
+ *
+ * \param j2: `int`
+ * \param j3: `int`
+ * \param m2: `int`
+ * \param m3: `int`
+ * \param c6: `C6 *`
+ */
+void r3jjr(int j2, int j3, int m2, int m3, C6 *c6);
+
+/*! \brief C++ porting of R3JMR
+ *
+ * \param j1: `int`
+ * \param j2: `int`
+ * \param j3: `int`
+ * \param m1: `int`
+ * \param c6: `C6 *`
+ */
+void r3jmr(int j1, int j2, int j3, int m1, C6 *c6);
 
 /*! \brief C++ porting of RABA
  *
@@ -1741,133 +269,7 @@ void r3j000(int j2, int j3, C6 *c6) {
 void raba(
 	  int le, std::complex<double> **am0m, std::complex<double> **w, double **tqce,
 	  std::complex<double> **tqcpe, double **tqcs, std::complex<double> **tqcps
-	  ) {
-  std::complex<double> **a, **ctqce, **ctqcs;
-  std::complex<double> acw, acwp, aca, acap, c1, c2, c3;
-  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);
-  int nlem = le * (le + 2);
-  const int nlemt = nlem + nlem;
-  a = new std::complex<double>*[nlemt];
-  ctqce = new std::complex<double>*[2];
-  ctqcs = new std::complex<double>*[2];
-  for (int ai = 0; ai < nlemt; ai++) a[ai] = new std::complex<double>[2]();
-  for (int ci = 0; ci < 2; ci++) {
-    ctqce[ci] = new std::complex<double>[3]();
-    ctqcs[ci] = new std::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;
-}
+);
 
 /*! \brief C++ porting of RFTR
  *
@@ -1890,17 +292,7 @@ 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];
-}
+);
 
 /*! \brief C++ porting of SCR0
  *
@@ -1911,51 +303,7 @@ void rftr(
  * \param c3: `C3 *` Pointer to a C3 instance.
  * \param c4: `C4 *` Pointer to a C4 structure.
  */
-void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
-  const std::complex<double> cc0(0.0, 0.0);
-  double exdc = exri * exri;
-  double ccs = 4.0 * acos(0.0) / (vk * vk);
-  double cccs = ccs / exdc;
-  std::complex<double> sum21, rm, re, csam;
-  csam = -(ccs / (exri * vk)) * std::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 scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4);
 
 /*! \brief C++ porting of SCR2
  *
@@ -1970,85 +318,8 @@ void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
  */
 void scr2(
 	  double vk, double vkarg, double exri, double *duk, C1 *c1, C1_AddOns *c1ao,
-	  C3 *c3, C4 *c4) {
-  const std::complex<double> cc0(0.0, 0.0);
-  const std::complex<double> uim(0.0, 1.0);
-  std::complex<double> s11, s21, s12, s22, rm, re, csam, cph, phas, cc;
-  double ccs = 1.0 / (vk * vk);
-  csam = -(ccs / (exri * vk)) * std::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
-}
+	  C3 *c3, C4 *c4
+);
 
 /*! \brief C++ porting of STR
  *
@@ -2059,81 +330,7 @@ void scr2(
  * \param c4: `C4 *` Pointer to a C4 structure.
  * \param c6: `C6 *` Pointer to a C6 instance.
  */
-void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
-  std::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 std::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 str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6);
 
 /*! \brief C++ porting of TQR
  *
@@ -2152,14 +349,7 @@ void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
 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];
-}
+);
 
 /*! \brief C++ porting of ZTM
  *
@@ -2170,98 +360,6 @@ void tqr(
  * \param c6: `C6 *` Pointer to a C6 instance.
  * \param c9: `C9 *` Pointer to a C9 instance.
  */
-void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9) {
-  std::complex<double> gie, gle, a1, a2, a3, a4, sum1, sum2, sum3, sum4;
-  const std::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
-}
-
-/*! \brief Sum all the elements of a matrix (debug function).
- *
- *  \param mat: Matrix of complex.
- *  \param rows: `int`
- *  \param cols: `int`
- */
-std::complex<double> summat(std::complex<double> **mat, int rows, int cols) {
-  std::complex<double> result(0.0, 0.0);
-  for (int i = 0; i < rows; i++) {
-    for (int j = 0; j < cols; j++) result += mat[i][j];
-  }
-  return result;
-}
+void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9 * c9);
 
 #endif
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index 781f679514c4857753348b45fea0498d4a9aadfe..e698cfdb19c37408e90543f7bde9969e9461af15 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -2,7 +2,7 @@
  *
  * \brief C++ porting of SPH functions and subroutines.
  *
- * This library includes a collection of function that are used to solve the
+ * This library includes a collection of functions that are used to solve the
  * scattering problem in the case of a single sphere. Some of these functions
  * are also generalized to the case of clusters of spheres. In most cases, the
  * results of calculations do not fall back to fundamental data types. They are
diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b1fb2d3def86e8f7aff84b6dadd1bc3b3c7c6a89
--- /dev/null
+++ b/src/libnptm/clu_subs.cpp
@@ -0,0 +1,1971 @@
+/*! \file clu_subs.cpp
+ *
+ * \brief C++ implementation of CLUSTER subroutines.
+ */
+
+#include "../include/clu_subs.h"
+
+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
+}