/*! \file clu_subs.cpp
*
* \brief C++ implementation of CLUSTER subroutines.
*/
#include <complex>
#ifndef INCLUDE_COMMONS_H_
#include "../include/Commons.h"
#endif
#ifndef INCLUDE_SPH_SUBS_H_
#include "../include/sph_subs.h"
#endif
#ifndef INCLUDE_CLU_SUBS_H_
#include "../include/clu_subs.h"
#endif
using namespace std;
void apc(
double ****zpv, int le, complex<double> **am0m, complex<double> **w,
double sqk, double **gapr, complex<double> **gapp
) {
complex<double> **ac, **gap;
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
complex<double> uimmp, summ, sume, suem, suee, summp, sumep;
complex<double> suemp, sueep;
double cof = 1.0 / sqk;
double cimu = cof / sqrt(2.0);
int nlem = le * (le + 2);
const int nlemt = nlem + nlem;
ac = new complex<double>*[nlemt];
gap = new complex<double>*[3];
for (int ai = 0; ai < nlemt; ai++) ac[ai] = new complex<double>[2]();
for (int gi = 0; gi < 3; gi++) gap[gi] = new complex<double>[2]();
for (int j45 = 1; j45 <= nlemt; j45++) {
int j = j45 - 1;
ac[j][0] = cc0;
ac[j][1] = cc0;
for (int i45 = 1; i45 <= nlemt; i45++) {
int i = i45 - 1;
ac[j][0] += (am0m[j][i] * w[i][0]);
ac[j][1] += (am0m[j][i] * w[i][1]);
} //i45 loop
} //j45 loop
for (int imu90 = 1; imu90 <=3; imu90++) {
int mu = imu90 - 2;
gap[imu90 - 1][0] = cc0;
gap[imu90 - 1][1] = cc0;
gapp[imu90 - 1][0] = cc0;
gapp[imu90 - 1][1] = cc0;
for (int l80 =1; l80 <= le; l80++) {
int lpo = l80 + 1;
int ltpo = lpo + l80;
int imm = l80 * lpo;
for (int ilmp = 1; ilmp <= 3; ilmp++) {
if ((l80 == 1 && ilmp == 1) || (l80 == le && ilmp == 3)) continue; // ilmp loop
int lmpml = ilmp - 2;
int lmp = l80 + lmpml;
uimmp = (-1.0 * lmpml) * uim;
int impmmmp = lmp * (lmp + 1);
for (int im70 = 1; im70 <= ltpo; im70++) {
int m = im70 - lpo;
int mmp = m - mu;
int abs_mmp = (mmp > 0) ? mmp : -mmp;
if (abs_mmp <= lmp) {
int i = imm + m;
int ie = i + nlem;
int imp = impmmmp + mmp;
int impe = imp + nlem;
double cgc = cg1(lmpml, mu, l80, m);
int jpo = 2;
for (int ipo = 1; ipo <= 2; ipo++) {
if (ipo == 2) jpo = 1;
summ = dconjg(ac[i - 1][ipo - 1]) * ac[imp - 1][ipo - 1];
sume = dconjg(ac[i - 1][ipo - 1]) * ac[impe - 1][ipo - 1];
suem = dconjg(ac[ie - 1][ipo - 1]) * ac[imp - 1][ipo - 1];
suee = dconjg(ac[ie - 1][ipo - 1]) * ac[impe - 1][ipo - 1];
summp = dconjg(ac[i - 1][jpo - 1]) * ac[imp - 1][ipo - 1];
sumep = dconjg(ac[i - 1][jpo - 1]) * ac[impe - 1][ipo - 1];
suemp = dconjg(ac[ie - 1][jpo - 1]) * ac[imp - 1][ipo - 1];
sueep = dconjg(ac[ie - 1][jpo - 1]) * ac[impe - 1][ipo - 1];
if (lmpml != 0) {
summ *= uimmp;
sume *= uimmp;
suem *= uimmp;
suee *= uimmp;
summp *= uimmp;
sumep *= uimmp;
suemp *= uimmp;
sueep *= uimmp;
}
// label 55
gap[imu90 - 1][ipo - 1] += (
(
summ * zpv[l80 - 1][ilmp - 1][0][0]
+ sume * zpv[l80 - 1][ilmp - 1][0][1]
+ suem * zpv[l80 - 1][ilmp - 1][1][0]
+ suee * zpv[l80 - 1][ilmp - 1][1][1]
) * cgc
);
gapp[imu90 - 1][ipo - 1] += (
(
summp * zpv[l80 - 1][ilmp - 1][0][0]
+ sumep * zpv[l80 - 1][ilmp - 1][0][1]
+ suemp * zpv[l80 - 1][ilmp - 1][1][0]
+ sueep * zpv[l80 - 1][ilmp - 1][1][1]
) * cgc
);
} // ipo loop
} // ends im70 loop
} // im70 loop
} // ilmp loop
} // l80 loop
} // imu90 loop
for (int ipo95 = 1; ipo95 <= 2; ipo95++) {
sume = gap[0][ipo95 - 1] * cimu;
suee = gap[1][ipo95 - 1] * cof;
suem = gap[2][ipo95 - 1] * cimu;
gapr[0][ipo95 - 1] = (sume - suem).real();
gapr[1][ipo95 - 1] = ((sume + suem) * uim).real();
gapr[2][ipo95 - 1] = suee.real();
sumep = gapp[0][ipo95 - 1] * cimu;
sueep = gapp[1][ipo95 - 1] * cof;
suemp = gapp[2][ipo95 - 1] * cimu;
gapp[0][ipo95 - 1] = sumep - suemp;
gapp[1][ipo95 - 1] = (sumep + suemp) * uim;
gapp[2][ipo95 - 1] = sueep;
} // ipo95 loop
// Clean memory
for (int ai = nlemt - 1; ai > -1; ai--) delete[] ac[ai];
for (int gi = 2; gi > -1; gi--) delete[] gap[gi];
delete[] ac;
delete[] gap;
}
void apcra(
double ****zpv, const int le, complex<double> **am0m, int inpol, double sqk,
double **gaprm, complex<double> **gappm
) {
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
complex<double> uimtl, uimtls, ca11, ca12, ca21, ca22;
complex<double> a11, a12, a21, a22, sum1, sum2, fc;
double ****svw = new double***[le];
complex<double> ****svs = new complex<double>***[le];
for (int i = 0; i < le; i++) {
svw[i] = new double**[3];
svs[i] = new complex<double>**[3];
for (int j = 0; j < 3; j++) {
svw[i][j] = new double*[2];
svs[i][j] = new complex<double>*[2];
for (int k = 0; k < 2; k++) {
svw[i][j][k] = new double[2]();
svs[i][j][k] = new complex<double>[2]();
}
}
}
int nlem = le * (le + 2);
for (int l28 = 1; l28 <= le; l28++) {
int lpo = l28 + 1;
int ltpo = lpo + l28;
double fl = sqrt(1.0 * ltpo);
for (int ilmp = 1; ilmp <= 3; ilmp++) {
if ((l28 == 1 && ilmp == 1) || (l28 == le && ilmp == 3)) continue; // ilmp loop
int lmpml = ilmp - 2;
int lmp = l28 + lmpml;
double flmp = sqrt(1.0 * (lmp + lmp + 1));
double fllmp = flmp / fl;
double cgmmo = fllmp * cg1(lmpml, 0, l28, 1);
double cgmpo = fllmp * cg1(lmpml, 0, l28, -1);
if (inpol == 0) {
double cgs = cgmpo + cgmmo;
double cgd = cgmpo - cgmmo;
svw[l28 - 1][ilmp - 1][0][0] = cgs;
svw[l28 - 1][ilmp - 1][0][1] = cgd;
svw[l28 - 1][ilmp - 1][1][0] = cgd;
svw[l28 - 1][ilmp - 1][1][1] = cgs;
} else { // label 22
svw[l28 - 1][ilmp - 1][0][0] = cgmpo;
svw[l28 - 1][ilmp - 1][1][0] = cgmpo;
svw[l28 - 1][ilmp - 1][0][1] = -cgmmo;
svw[l28 - 1][ilmp - 1][1][1] = cgmmo;
}
// label 26
} // ilmp loop
} // l28 loop
for (int l30 = 1; l30 <= le; l30++) { // 0-init: can be omitted
for (int ilmp = 1; ilmp <= 3; ilmp++) {
for (int ipa = 1; ipa <= 2; ipa++) {
for (int ipamp = 1; ipamp <= 2; ipamp++) {
svs[l30 - 1][ilmp - 1][ipa - 1][ipamp - 1] = cc0;
}
} // ipa loop
} // ilmp loop
} // l30 loop
for (int l58 = 1; l58 <= le; l58 ++) {
int lpo = l58 + 1;
int ltpo = l58 + lpo;
int imm = l58 * lpo;
for (int ilmp = 1; ilmp <= 3; ilmp++) {
if ((l58 == 1 && ilmp == 1) || (l58 == le && ilmp == 3)) continue; // ilmp loop
int lmpml = ilmp - 2;
int lmp = l58 + lmpml;
int impmm = lmp * (lmp + 1);
uimtl = uim * (1.0 * lmpml);
if (lmpml == 0) uimtl = complex<double>(1.0, 0.0);
for (int im54 = 1; im54 <= ltpo; im54++) {
int m = im54 - lpo;
int i = imm + m;
int ie = i + nlem;
for (int imu52 = 1; imu52 <= 3; imu52++) {
int mu = imu52 - 2;
int mmp = m - mu;
int abs_mmp = (mmp > 0) ? mmp : -mmp;
if (abs_mmp <= lmp) {
int imp = impmm + mmp;
int impe = imp + nlem;
double cgc = cg1(lmpml, -mu, l58, -m);
for (int ls = 1; ls <= le; ls++) {
int lspo = ls + 1;
int lstpo = ls + lspo;
int ismm = ls * lspo;
for (int ilsmp = 1; ilsmp <= 3; ilsmp++) {
if ((ls == 1 && ilsmp == 1) || (ls == le && ilsmp == 3)) continue; // ilsmp loop
int lsmpml = ilsmp - 2;
int lsmp = ls + lsmpml;
int ismpmm = lsmp * (lsmp + 1);
uimtls = -uim * (1.0 * lsmpml);
if (lsmpml == 0) uimtls = complex<double>(1.0, 0.0);
for (int ims = 1; ims <= lstpo; ims++) {
int ms = ims - lspo;
int msmp = ms - mu;
int abs_msmp = (msmp > 0) ? msmp : -msmp;
if (abs_msmp <= lsmp) {
int is = ismm + ms;
int ise = is + nlem;
int ismp = ismpmm + msmp;
int ismpe = ismp + nlem;
double cgcs = cg1(lsmpml, mu, ls, ms);
fc = (uimtl * uimtls) * (cgc * cgcs);
ca11 = dconjg(am0m[is - 1][i - 1]);
ca12 = dconjg(am0m[is - 1][ie - 1]);
ca21 = dconjg(am0m[ise - 1][i - 1]);
ca22 = dconjg(am0m[ise - 1][ie - 1]);
a11 = am0m[ismp - 1][imp - 1];
a12 = am0m[ismp - 1][impe - 1];
a21 = am0m[ismpe - 1][imp - 1];
a22 = am0m[ismpe - 1][impe - 1];
double z11 = zpv[ls - 1][ilsmp - 1][0][0];
double z12 = zpv[ls - 1][ilsmp - 1][0][1];
double z21 = zpv[ls - 1][ilsmp - 1][1][0];
double z22 = zpv[ls - 1][ilsmp - 1][1][1];
svs[l58 - 1][ilmp - 1][0][0] += ((ca11 * a11 * z11
+ ca11 * a21 * z12
+ ca21 * a11 * z21
+ ca21 * a21 * z22) * fc);
svs[l58 - 1][ilmp - 1][0][1] += ((ca11 * a12 * z11
+ ca11 * a22 * z12
+ ca21 * a12 * z21
+ ca21 * a22 * z22) * fc);
svs[l58 - 1][ilmp - 1][1][0] += ((ca12 * a11 * z11
+ ca12 * a21 * z12
+ ca22 * a11 * z21
+ ca22 * a21 * z22) * fc);
svs[l58 - 1][ilmp - 1][1][1] += ((ca12 * a12 * z11
+ ca12 * a22 * z12
+ ca22 * a12 * z21
+ ca22 * a22 * z22) * fc);
} // ends ims loop
} // ims loop
} // ilsmp loop
} // ls loop
} // ends imu52 loop
} // imu52 loop
} // im54 loop
} // ilmp loop
} // l58 loop
sum1 = cc0;
sum2 = cc0;
for (int l68 = 1; l68 <= le; l68++) {
//int lpo = l68 + 1;
//int ltpo = l68 + lpo;
//int imm = l68 * lpo;
for (int ilmp = 1; ilmp <= 3; ilmp++) {
if ((l68 == 1 && ilmp == 1) || (l68 == le && ilmp == 3)) continue; // ilmp loop
if (inpol == 0) {
sum1 += (
svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][0][0]
+ svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][0][1]
+ svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][1][0]
+ svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][1][1]
);
sum2 += (
svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][0][0]
+ svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][0][1]
+ svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][1][0]
+ svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][1][1]
);
} else { // label 62
sum1 += (
svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][0][0]
+ svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][0][1]
+ svw[l68 - 1][ilmp - 1][0][0] * svs[l68 - 1][ilmp - 1][1][0]
+ svw[l68 - 1][ilmp - 1][1][0] * svs[l68 - 1][ilmp - 1][1][1]
);
sum2 += (
svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][0][0]
+ svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][0][1]
+ svw[l68 - 1][ilmp - 1][0][1] * svs[l68 - 1][ilmp - 1][1][0]
+ svw[l68 - 1][ilmp - 1][1][1] * svs[l68 - 1][ilmp - 1][1][1]
);
} // label 66, ends ilmp loop
} // ilmp loop
} // l68 loop
const double half_pi = acos(0.0);
double cofs = half_pi * 2.0 / sqk;
gaprm[0][0] = 0.0;
gaprm[0][1] = 0.0;
gaprm[1][0] = 0.0;
gaprm[1][1] = 0.0;
gappm[0][0] = cc0;
gappm[0][1] = cc0;
gappm[1][0] = cc0;
gappm[1][1] = cc0;
if (inpol == 0) {
sum1 *= cofs;
sum2 *= cofs;
gaprm[2][0] = sum1.real();
gaprm[2][1] = sum1.real();
gappm[2][0] = sum2 * uim;
gappm[2][1] = -gappm[2][0];
} else { // label 72
cofs *= 2.0;
gaprm[2][0] = sum1.real() * cofs;
gaprm[2][1] = sum2.real() * cofs;
gappm[2][0] = cc0;
gappm[2][1] = cc0;
}
// Clean memory
for (int i = le - 1; i > -1; i--) {
for (int j = 2; j > -1; j--) {
for (int k = 1; k > -1; k--) {
delete[] svw[i][j][k];
delete[] svs[i][j][k];
}
delete[] svw[i][j];
delete[] svs[i][j];
}
delete[] svw[i];
delete[] svs[i];
}
delete[] svw;
delete[] svs;
}
complex<double> cdtp(
complex<double> z, complex<double> **am, int i, int jf,
int k, int nj
) {
/* NOTE: the original FORTRAN code treats the AM matrix as a
* vector. This is not directly allowed in C++ and it requires
* accounting for the different dimensions.
*/
complex<double> result = z;
if (nj > 0) {
int jl = jf + nj - 1;
for (int j = jf; j <= jl; j++) {
result += (am[i - 1][j - 1] * am[j - 1][k - 1]);
}
}
return result;
}
double cgev(int ipamo, int mu, int l, int m) {
double result = 0.0;
double xd = 0.0, xn = 0.0;
if (ipamo == 0) {
if (m != 0 || mu != 0) { // label 10
if (mu != 0) {
xd = 2.0 * l * (l + 1);
if (mu <= 0) {
xn = 1.0 * (l + m) * (l - m + 1);
result = sqrt(xn / xd);
} else { // label 15
xn = 1.0 * (l - m) * (l + m + 1);
result = -sqrt(xn / xd);
}
} else { // label 20
xd = 1.0 * (l + 1) * l;
xn = -1.0 * m;
result = xn / sqrt(xd);
}
}
} else { // label 30
xd = 2.0 * l * (l * 2 - 1);
if (mu < 0) { // label 35
xn = 1.0 * (l - 1 + m) * (l + m);
} else if (mu == 0) { // label 40
xn = 2.0 * (l - m) * (l + m);
} else { // mu > 0, label 45
xn = 1.0 * (l - 1 - m) * (l - m);
}
result = sqrt(xn / xd);
}
return result;
}
void cms(complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
complex<double> dm, de, cgh, cgk;
const complex<double> cc0(0.0, 0.0);
int ndi = c4->nsph * c4->nlim;
int nbl = 0;
int nsphmo = c4->nsph - 1;
for (int n1 = 1; n1 <= nsphmo; n1++) { // GPU portable?
int in1 = (n1 - 1) * c4->nlim;
int n1po = n1 + 1;
for (int n2 = n1po; n2 <= c4->nsph; n2++) {
int in2 = (n2 - 1) * c4->nlim;
nbl++;
for (int l1 = 1; l1 <= c4->li; l1++) {
int l1po = l1 + 1;
int il1 = l1po * l1;
int l1tpo = l1po + l1;
for (int im1 = 1; im1 <= l1tpo; im1++) {
int m1 = im1 - l1po;
int ilm1 = il1 + m1;
int ilm1e = ilm1 + ndi;
int i1 = in1 + ilm1;
int i1e = in1 + ilm1e;
int j1 = in2 + ilm1;
int j1e = in2 + ilm1e;
for (int l2 = 1; l2 <= c4->li; l2++) {
int l2po = l2 + 1;
int il2 = l2po * l2;
int l2tpo = l2po + l2;
int ish = ((l2 + l1) % 2 == 0) ? 1 : -1;
int isk = -ish;
for (int im2 = 1; im2 <= l2tpo; im2++) {
int m2 = im2 - l2po;
int ilm2 = il2 + m2;
int ilm2e = ilm2 + ndi;
int i2 = in2 + ilm2;
int i2e = in2 + ilm2e;
int j2 = in1 + ilm2;
int j2e = in1 + ilm2e;
cgh = ghit(0, 0, nbl, l1, m1, l2, m2, c1, c1ao, c4, c6);
cgk = ghit(0, 1, nbl, l1, m1, l2, m2, c1, c1ao, c4, c6);
am[i1 - 1][i2 - 1] = cgh;
am[i1 - 1][i2e - 1] = cgk;
am[i1e - 1][i2 - 1] = cgk;
am[i1e - 1][i2e - 1] = cgh;
am[j1 - 1][j2 - 1] = cgh * (1.0 * ish);
am[j1 - 1][j2e - 1] = cgk * (1.0 * isk);
am[j1e - 1][j2 - 1] = cgk * (1.0 * isk);
am[j1e - 1][j2e - 1] = cgh * (1.0 * ish);
}
}
} // im1 loop
} // l1 loop
} // n2 loop
} // n1 loop
for (int n1 = 1; n1 <= c4->nsph; n1++) { // GPU portable?
int in1 = (n1 - 1) * c4->nlim;
for (int l1 = 1; l1 <= c4->li; l1++) {
dm = c1->rmi[l1 - 1][n1 - 1];
de = c1->rei[l1 - 1][n1 - 1];
int l1po = l1 + 1;
int il1 = l1po * l1;
int l1tpo = l1po + l1;
for (int im1 = 1; im1 <= l1tpo; im1++) {
int m1 = im1 - l1po;
int ilm1 = il1 + m1;
int i1 = in1 + ilm1;
int i1e = i1 + ndi;
for (int ilm2 = 1; ilm2 <= c4->nlim; ilm2++) {
int i2 = in1 + ilm2;
int i2e = i2 + ndi;
am[i1 - 1][i2 - 1] = cc0;
am[i1 - 1][i2e - 1] = cc0;
am[i1e - 1][i2 - 1] = cc0;
am[i1e - 1][i2e - 1] = cc0;
}
am[i1 - 1][i1 - 1] = dm;
am[i1e - 1][i1e - 1] = de;
} // im1 loop
} // l1 loop
} // n1 loop
}
void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
complex<double> ***svf, ***svw, **svs;
const complex<double> cc0(0.0, 0.0);
complex<double> cam(0.0, 0.0);
const int le4po = 4 * c4->le + 1;
svf = new complex<double>**[le4po];
svw = new complex<double>**[le4po];
svs = new complex<double>*[le4po];
for (int si = 0; si < le4po; si++) {
svf[si] = new complex<double>*[le4po];
svw[si] = new complex<double>*[4];
svs[si] = new complex<double>[4]();
for (int sj = 0; sj < le4po; sj++) svf[si][sj] = new complex<double>[4]();
for (int sj = 0; sj < 4; sj++) svw[si][sj] = new complex<double>[4]();
}
double exdc = exri * exri;
double ccs = 1.0 / (vk * vk);
const double pi4sq = 64.0 * acos(0.0) * acos(0.0);
double cint = ccs / (pi4sq * exdc);
int letpo = c4->le + c4->le + 1;
for (int i20 = 0; i20 < 16; i20++) c1ao->vintm[i20] = cc0; // 0-init: can be omitted
for (int lpo40 = 1; lpo40 <= letpo; lpo40++) {
int l = lpo40 - 1;
int ltpo = lpo40 + l;
int immn = letpo - l;
int immx = letpo + l;
for (int imf = immn; imf <= immx; imf++) { // 0-init: can be omitted
for (int ims = immn; ims <= immx; ims++) {
for (int ipo = 1; ipo <= 4; ipo++) {
svf[imf - 1][ims - 1][ipo - 1] = cc0;
} // ipo loop
} // ims loop
} // imf loop
for (int l1 = 1; l1 <= c4->le; l1++) {
int il1 = l1 * (l1 + 1);
for (int l2 = 1; l2 <= c4->le; l2++) {
int abs_l2ml1 = (l2 > l1) ? l2 - l1 : l1 - l2;
if (l < abs_l2ml1 || l > l2 + l1) continue; // l2 loop
int il2 = l2 * (l2 + 1);
for (int im = immn; im >= immx; im++) { // 0-init: can be omitted
for (int ipa = 1; ipa <= 4; ipa++) {
svs[im - 1][ipa - 1] = cc0;
for (int ipo = 1; ipo <= 4; ipo++) {
svw[im - 1][ipa - 1][ipo - 1] = cc0;
} // ipo loop
} // ipa loop
} // im loop
for (int im = immn; im <= immx; im++) {
int m = im - letpo;
r3jmr(l, l1, l2, m, c6);
int m1mnmo = (-l1 > -l2 - m) ? -(l1 + 1) : -(l2 + m + 1);
int nm1 = (l1 < l2 - m) ? (l1 - m1mnmo) : (l2 - m - m1mnmo);
for (int im1 = 1; im1 <= nm1; im1++) {
int m1 = -im1 - m1mnmo;
int isn = 1;
if (m1 % 2 != 0) isn = -1;
double cg3j = c6->rac3j[im1 - 1] * isn;
int ilm1 = il1 + m1;
int ilm2 = il2 + m1 - m;
int ipa = 0;
for (int ipa1 = 1; ipa1 <= 2; ipa1++) {
int i1 = ilm1;
if (ipa1 == 2) i1 = ilm1 + c4->nlem;
for (int ipa2 = 1; ipa2 <= 2; ipa2++) {
int i2 = ilm2;
if (ipa2 == 2) i2 = ilm2 + c4->nlem;
ipa++;
svs[im - 1][ipa - 1] += (c1ao->am0m[i1 - 1][i2 - 1] * cg3j);
int ipo = 0;
for (int ipo2 = 1; ipo2 <= 2; ipo2++) {
for (int ipo1 = 3; ipo1 <= 4; ipo1++) {
ipo++;
svw[im - 1][ipa - 1][ipo - 1] += (c1->w[i1 - 1][ipo1 - 1] * c1->w[i2 - 1][ipo2 - 1] * cg3j);
} // ipo1 loop
} // ipo2 loop
} // ipa2 loop
} // ipa1 loop
} // im1 loop
// label 32 loops
for (int imf = immn; imf <= immx; imf++) {
for (int ims = immn; ims <= immx; ims++) {
for (int ipo = 1; ipo <= 4; ipo++) {
for (int ipa = 1; ipa <= 4; ipa++) {
svf[imf - 1][ims - 1][ipo - 1] += (svw[imf - 1][ipa - 1][ipo - 1] * svs[ims - 1][ipa - 1]);
} // ipa loop
} // ipo loop
} // ims loop
} // imf loop
// ends loop level 34, which are l2 loop and l1 loop
} // im loop
} // l2 loop
} // l1 loop
for (int imf = immn; imf <= immx; imf++) {
for (int ims = immn; ims <= immx; ims++) {
int i = 0;
for (int ipo1 = 1; ipo1 <= 4; ipo1++) {
cam = dconjg(svf[imf - 1][ims - 1][ipo1 - 1]);
for (int ipo2 = 1; ipo2 <= 4; ipo2++) {
i++;
c1ao->vintm[i - 1] += (svf[imf - 1][ims - 1][ipo2 - 1] * cam * (1.0 * ltpo));
}
} // ipo1 loop
} // ims loop
} // imf loop
} // lpo40 loop
for (int i42 = 0; i42 < 16; i42++) c1ao->vintm[i42] *= cint;
// Clean memory
for (int si = le4po - 1; si > -1; si--) {
for (int sj = le4po - 1; sj > -1; sj--) delete[] svf[si][sj];
for (int sj = 3; sj > -1; sj--) delete[] svw[si][sj];
delete[] svf[si];
delete[] svw[si];
delete[] svs[si];
}
delete[] svf;
delete[] svw;
delete[] svs;
}
complex<double> ghit(
int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
C1_AddOns *c1ao, C4 *c4, C6 *c6
) {
/* NBL identifies transfer vector going from N2 to N1;
* IHI=0 for Hankel, IHI=1 for Bessel, IHI=2 for Bessel from origin;
* depending on IHI, IPAM=0 gives H or I, IPAM= 1 gives K or L. */
const complex<double> cc0(0.0, 0.0);
const complex<double> uim(0.0, 1.0);
complex<double> csum(0.0, 0.0), cfun(0.0, 0.0);
complex<double> result = cc0;
if (ihi == 2) {
if (c1->rxx[nbl - 1] == 0.0 && c1->ryy[nbl - 1] == 0.0 && c1->rzz[nbl - 1] == 0.0) {
if (ipamo == 0) {
if (l1 == l2 && m1 == m2) result = complex<double>(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
}