diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index ab8fff17f926789b629f1bd82f06044e97532798..9e13b5267940c34b75059b578a497f93372b9408 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -27,36 +27,36 @@ 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
-);
+ 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);
extern void sphar(double cth, double sth, double cph, double sph, int lm, std::complex<double> *ylm);
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
-);
+ 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
-);
+ 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
-);
+ 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
-);
+ 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
@@ -69,21 +69,21 @@ extern void wmasp(
* \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;
+ 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
@@ -95,37 +95,37 @@ std::complex<double> cdtp(
* \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);
+ 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);
}
- return result;
+ } 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
@@ -137,123 +137,123 @@ double cgev(int ipamo, int mu, int l, int m) {
* \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;
- }
+ 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
@@ -265,98 +265,98 @@ void r3jjr(int j2, int j3, int m2, int m3, C6 *c6) {
* \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
- }
+ 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
@@ -374,216 +374,216 @@ void r3jmr(int j1, int j2, int j3, int m1, C6 *c6) {
* \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;
+ 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;
- }
+ 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;
}
- // 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;
+ // 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;
}
- // 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
+ // 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 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;
+ // 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
@@ -597,123 +597,123 @@ std::complex<double> ghit(
* \param gapp: Matrix of complex.
*/
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;
+ 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
@@ -727,223 +727,223 @@ void apc(
* \param gappm: Matrix of complex.
*/
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;
+ 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;
}
- // 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;
+ } // 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 CMS
@@ -955,295 +955,295 @@ void apcra(
* \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
+ 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
}
/*! \brief C++ porting of CRSM1
*
- * \param vk: `double`
- * \param exri: `double`
+ * \param vk: `double` Wave number.
+ * \param exri: `double` External medium refractive index.
* \param c1: `C1 *`
* \param c1ao: `C1_AddOns *`
* \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;
+ 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;
- // 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;
+ // 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 HJV
*
- * \param exri: `double`
- * \param vk: `double`
- * \param jer: `int &`
- * \param lcalc: `int &`
+ * \param exri: `double` External medium refractive index.
+ * \param vk: `double` Wave number.
+ * \param jer: `int &` Reference to error code flag.
+ * \param lcalc: `int &` Reference to the highest order accounted for in calculation.
* \param arg: `complex\<double\> &`
* \param c1: `C1 *`
* \param c1ao: `C1_AddOns *`
* \param c4: `C4 *`
*/
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
+ 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
@@ -1254,380 +1254,380 @@ void hjv(
* \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;
+ /* 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;
}
/*! \brief C++ porting of MEXTC
*
- * \param vk: `double`
- * \param exri: `double`
+ * \param vk: `double` Wave number.
+ * \param exri: `double` External medium refractive index.
* \param fsac: Matrix of complex
* \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;
- }
- }
+ 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;
+ }
+ }
}
/*! \brief C++ porting of PCROS
*
* This function is intended to evaluate the particle cross-section. QUESTIUON: correct?
- * \param vk: `double`
- * \param exri: `double`
+ * \param vk: `double` Wave number.
+ * \param exri: `double` External medium refractive index.
* \param c1: `C1 *`
* \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
+ 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
}
/*! \brief C++ porting of PCRSM0
*
- * \param vk: `double`
- * \param exri: `double`
- * \param inpol: `int`
+ * \param vk: `double` Wave number.
+ * \param exri: `double` External medium refractive index.
+ * \param inpol: `int` Incident field polarization type.
* \param c1: `C1 *`
* \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;
+ 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;
}
/*! \brief C++ porting of POLAR
*
- * \param x: `double`
- * \param y: `double`
- * \param z: `double`
- * \param r: `double &`
- * \param cth: `double &`
- * \param sth: `double &`
- * \param cph: `double &`
- * \param sph: `double &`
+ * \param x: `double` X-axis Cartesian coordinate.
+ * \param y: `double` Y-axis Cartesian coordinate.
+ * \param z: `double` Z-axis Cartesian coordinate.
+ * \param r: `double &` Reference to radial vector (output value).
+ * \param cth: `double &` Reference to the cosine of the azimuth coordinate (output value).
+ * \param sth: `double &` Reference to the sine of the azimuth coordinate (output value).
+ * \param cph: `double &` Reference to the cosine of the elevation coordinate (output value).
+ * \param sph: `double &` Reference to the sine of the elevation coordinate (output value).
*/
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
- }
- }
+ 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
@@ -1637,95 +1637,95 @@ void polar(
* \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;
- }
+ 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;
+ }
}
/*! \brief C++ porting of RABA
@@ -1739,134 +1739,134 @@ void r3j000(int j2, int j3, C6 *c6) {
* \param tqcps: Matrix of complex.
*/
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]();
+ 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);
}
- 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;
+ // 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
@@ -1887,81 +1887,81 @@ void raba(
* \param fz: `double &`
*/
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];
+ 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
*
- * \param vk: `double` QUESTION: definition?
- * \param exri: `double` External medium refractive index. QUESTION: correct?
+ * \param vk: `double` Wave number
+ * \param exri: `double` External medium refractive index.
* \param c1: `C1 *` Pointer to a C1 instance.
* \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.
* \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
+ 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
}
/*! \brief C++ porting of SCR2
*
- * \param vk: `double` QUESTION: definition?
+ * \param vk: `double` Wave number.
* \param vkarg: `double` QUESTION: definition?
- * \param exri: `double` External medium refractive index
+ * \param exri: `double` External medium refractive index.
* \param duk: `double *` QUESTION: definition?
* \param c1: `C1 *` Pointer to a C1 instance.
* \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.
@@ -1969,85 +1969,85 @@ void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
* \param c4: `C4 *` Pointer to a C4 structure.
*/
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
+ 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
}
/*! \brief C++ porting of STR
@@ -2060,79 +2060,79 @@ void scr2(
* \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;
+ 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;
}
/*! \brief C++ porting of TQR
@@ -2150,15 +2150,15 @@ void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
* \param tsk: `double &`
*/
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];
+ 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
@@ -2171,83 +2171,83 @@ void tqr(
* \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
+ 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).
@@ -2257,11 +2257,11 @@ void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9
* \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;
+ 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;
}
#endif