diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index d7394ecdb0ff70cb255758d70f69569f3c445d8e..dd3470da68ed49ebc8f8fe31d9a4b9ef22f05072 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -110,46 +110,48 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
std::complex<double> summ, sume, suem, suee, sum;
double half_pi = acos(0.0);
double cofs = half_pi * 2.0 / sqk;
- for (int i40 = 1; i40 <= nsph; i40++) {
- int iogi = c1->iog[i40 - 1];
- if (iogi >= i40) {
+ for (int i40 = 0; i40 < nsph; i40++) {
+ int i = i40 + 1;
+ int iogi = c1->iog[i40];
+ if (iogi >= i) {
sum = cc0;
- for (int l30 = 1; l30 <= li; l30++) {
- int ltpo = l30 + l30 + 1;
- for (int ilmp = 1; ilmp <= 3; ilmp++) {
- bool goto30 = (l30 == 1 && ilmp == 1) || (l30 == li && ilmp == 3);
+ for (int l30 = 0; l30 < li; l30++) {
+ int l = l30 + 1;
+ int ltpo = l + l + 1;
+ for (int ilmp = 1; ilmp < 4; ilmp++) {
+ int ilmp30 = ilmp - 1;
+ bool goto30 = (l == 1 && ilmp == 1) || (l == li && ilmp == 3);
if (!goto30) {
int lmpml = ilmp - 2;
- int lmp = l30 + lmpml;
+ int lmp = l + lmpml;
double cofl = sqrt(1.0 * (ltpo * (lmp + lmp + 1)));
- summ = zpv[l30 - 1][ilmp - 1][0][0] /
+ summ = zpv[l30][ilmp30][0][0] /
(
- dconjg(c1->rmi[l30 - 1][i40 - 1]) *
- c1->rmi[lmp - 1][i40 - 1]
+ dconjg(c1->rmi[l30][i40]) *
+ c1->rmi[lmp - 1][i40]
);
- sume = zpv[l30 - 1][ilmp - 1][0][1] /
+ sume = zpv[l30][ilmp30][0][1] /
(
- dconjg(c1->rmi[l30 - 1][i40 - 1]) *
- c1->rei[lmp - 1][i40 - 1]
+ dconjg(c1->rmi[l30][i40]) *
+ c1->rei[lmp - 1][i40]
);
- suem = zpv[l30 - 1][ilmp - 1][1][0] /
+ suem = zpv[l30][ilmp30][1][0] /
(
- dconjg(c1->rei[l30 - 1][i40 - 1]) *
- c1->rmi[lmp - 1][i40 - 1]
+ dconjg(c1->rei[l30][i40]) *
+ c1->rmi[lmp - 1][i40]
);
- suee = zpv[l30 - 1][ilmp - 1][1][1] /
+ suee = zpv[l30][ilmp30][1][1] /
(
- dconjg(c1->rei[l30 - 1][i40 - 1]) *
- c1->rei[lmp - 1][i40 - 1]
+ dconjg(c1->rei[l30][i40]) *
+ c1->rei[lmp - 1][i40]
);
- sum += (cg1(lmpml, 0, l30, -1) * (summ - sume - suem + suee) +
- cg1(lmpml, 0, l30, 1) * (summ + sume + suem + suee)) * cofl;
+ sum += (cg1(lmpml, 0, l, -1) * (summ - sume - suem + suee) +
+ cg1(lmpml, 0, l, 1) * (summ + sume + suem + suee)) * cofl;
}
}
}
}
- gaps[i40 - 1] = sum.real() * cofs;
- //printf("DEBUG: gaps[%d] = %lE\n", i40, gaps[i40 - 1]);
+ gaps[i40] = sum.real() * cofs;
}
}
@@ -164,16 +166,19 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
* \param c2: `C2 *`
*/
void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) {
- double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1];
- double half_step = 0.5 * dif / npntmo;
+ const double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1];
+ const double half_step = 0.5 * dif / npntmo;
double rr = c1->rc[i - 1][ns - 1];
- std::complex<double> delta = c2->dc0[ic] - c2->dc0[ic - 1];
- int kpnt = npntmo + npntmo;
+ const std::complex<double> delta = c2->dc0[ic] - c2->dc0[ic - 1];
+ const int kpnt = npntmo + npntmo;
c2->ris[kpnt] = c2->dc0[ic];
c2->dlri[kpnt] = std::complex<double>(0.0, 0.0);
+ const int i90 = i - 1;
+ const int ns90 = ns - 1;
+ const int ic90 = ic - 1;
for (int np90 = 0; np90 < kpnt; np90++) {
- double ff = (rr - c1->rc[i - 1][ns - 1]) / dif;
- c2->ris[np90] = delta * ff * ff * (-2.0 * ff + 3.0) + c2->dc0[ic - 1];
+ double ff = (rr - c1->rc[i90][ns90]) / dif;
+ c2->ris[np90] = delta * ff * ff * (-2.0 * ff + 3.0) + c2->dc0[ic90];
c2->dlri[np90] = 3.0 * delta * ff * (1.0 - ff) / (dif * vk * c2->ris[np90]);
rr += half_step;
}
@@ -324,8 +329,8 @@ void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
cf0 = cf1;
cf1 = cf;
}
- double abs_csa = sqrt(csa.real() * csa.real() + csa.imag() * csa.imag());
- double abs_csb = sqrt(csb.real() * csb.real() + csb.imag() * csb.imag());
+ double abs_csa = abs(csa);
+ double abs_csb = abs(csb);
if (abs_csa > abs_csb) cs = csa / cf;
else cs = csb / cf0;
for (int k = 0; k <= nm; k++) {
@@ -430,33 +435,31 @@ void pwma(
double *up, double *un, std::complex<double> *ylm, int inpol, int lw,
int isq, C1 *c1
) {
- //std::complex<double> cp1, cm1, cc1, cp2, c2, cc2, uim, cl;
- double four_pi = 8.0 * acos(0.0);
+ const double four_pi = 8.0 * acos(0.0);
int is = isq;
if (isq == -1) is = 0;
int ispo = is + 1;
int ispt = is + 2;
int nlwm = lw * (lw + 2);
int nlwmt = nlwm + nlwm;
- double sqrtwi = 1.0 / sqrt(2.0);
- std::complex<double> uim(0.0, 1.0);
+ const double sqrtwi = 1.0 / sqrt(2.0);
+ const std::complex<double> uim(0.0, 1.0);
std::complex<double> cm1 = 0.5 * std::complex<double>(up[0], up[1]);
std::complex<double> cp1 = 0.5 * std::complex<double>(up[0], -up[1]);
double cz1 = up[2];
std::complex<double> cm2 = 0.5 * std::complex<double>(un[0], un[1]);
std::complex<double> cp2 = 0.5 * std::complex<double>(un[0], -un[1]);
double cz2 = un[2];
- //printf("DEBUG: in PWMA YLM(2) = (%lE, %lE)\n", ylm[1].real(), ylm[1].imag());
- for (int l20 = 1; l20 <= lw; l20++) {
- //if (ispo == 1) printf("DEBUG: l20 = %d\n", l20);
- int lf = l20 + 1;
- int lftl = lf * l20;
+ for (int l20 = 0; l20 < lw; l20++) {
+ int l = l20 + 1;
+ int lf = l + 1;
+ int lftl = lf * l;
double x = 1.0 * lftl;
std::complex<double> cl = std::complex<double>(four_pi / sqrt(x), 0.0);
- for (int powi = 1; powi <= l20; powi++) cl *= uim;
- int mv = l20 + lf;
+ for (int powi = 1; powi <= l; powi++) cl *= uim;
+ int mv = l + lf;
int m = -lf;
- for (int mf20 = 1; mf20 <= mv; mf20++) {
+ for (int mf20 = 0; mf20 < mv; mf20++) {
m += 1;
int k = lftl + m;
x = 1.0 * (lftl - m * (m + 1));
@@ -469,35 +472,27 @@ void pwma(
cm1 * cm * ylm[k - 1] +
cz1 * cz * ylm[k]
) * cl;
- //if (ispo == 1) printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k, ispo, c1->w[k - 1][ispo - 1].real(), c1->w[k - 1][ispo - 1].imag());
c1->w[k - 1][ispt - 1] = dconjg(
cp2 * cp * ylm[k + 1] +
cm2 * cm * ylm[k - 1] +
cz2 * cz * ylm[k]
) * cl;
- //printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k, ispt, c1->w[k - 1][ispt - 1].real(), c1->w[k - 1][ispt - 1].imag());
}
}
- for (int k30 = 1; k30 <= nlwm; k30++) {
+ for (int k30 = 0; k30 < nlwm; k30++) {
int i = k30 + nlwm;
- c1->w[i - 1][ispo - 1] = uim * c1->w[k30 - 1][ispt - 1];
- //printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispo, c1->w[i - 1][ispo - 1].real(), c1->w[i - 1][ispo - 1].imag());
- c1->w[i - 1][ispt - 1] = -uim * c1->w[k30 - 1][ispo - 1];
- //printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispt, c1->w[i - 1][ispt - 1].real(), c1->w[i - 1][ispt - 1].imag());
+ c1->w[i][ispo - 1] = uim * c1->w[k30][ispt - 1];
+ c1->w[i][ispt - 1] = -uim * c1->w[k30][ispo - 1];
}
if (inpol != 0) {
- for (int k40 = 1; k40 <= nlwm; k40++) {
+ for (int k40 = 0; k40 < nlwm; k40++) {
int i = k40 + nlwm;
- std::complex<double> cc1 = sqrtwi * (c1->w[k40 - 1][ispo - 1] + uim * c1->w[k40 - 1][ispt - 1]);
- std::complex<double> cc2 = sqrtwi * (c1->w[k40 - 1][ispo - 1] - uim * c1->w[k40 - 1][ispt - 1]);
- c1->w[k40 - 1][ispo - 1] = cc2;
- //printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k40, ispo, c1->w[k40 - 1][ispo - 1].real(), c1->w[k40 - 1][ispo - 1].imag());
- c1->w[i - 1][ispo - 1] = -cc2;
- //printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispo, c1->w[i - 1][ispo - 1].real(), c1->w[i - 1][ispo - 1].imag());
- c1->w[k40 - 1][ispt - 1] = cc1;
- //printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k40, ispt, c1->w[k40 - 1][ispt - 1].real(), c1->w[k40 - 1][ispt - 1].imag());
- c1->w[i - 1][ispt - 1] = cc1;
- //printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", i, ispt, c1->w[i - 1][ispt - 1].real(), c1->w[i - 1][ispt - 1].imag());
+ std::complex<double> cc1 = sqrtwi * (c1->w[k40][ispo - 1] + uim * c1->w[k40][ispt - 1]);
+ std::complex<double> cc2 = sqrtwi * (c1->w[k40][ispo - 1] - uim * c1->w[k40][ispt - 1]);
+ c1->w[k40][ispo - 1] = cc2;
+ c1->w[i][ispo - 1] = -cc2;
+ c1->w[k40][ispt - 1] = cc1;
+ c1->w[i][ispt - 1] = cc1;
}
} else {
if (isq == 0) {
@@ -505,16 +500,14 @@ void pwma(
}
}
if (isq != 0) {
- for (int i50 = 1; i50 <= 2; i50++) {
+ for (int i50 = 0; i50 < 2; i50++) {
int ipt = i50 + 2;
int ipis = i50 + is;
- for (int k50 = 1; k50 <= nlwmt; k50++) {
- c1->w[k50 - 1][ipt - 1] = dconjg(c1->w[k50 - 1][ipis - 1]);
- //printf("DEBUG: W( %d , %d) = (%lE,%lE)\n", k50, ipt, c1->w[k50 - 1][ipt - 1].real(), c1->w[k50 - 1][ipt - 1].imag());
+ for (int k50 = 0; k50 < nlwmt; k50++) {
+ c1->w[k50][ipt] = dconjg(c1->w[k50][ipis]);
}
}
}
- //printf("DEBUG: out PWMA W(1,1) = (%lE, %lE)\n", c1->w[0][0].real(), c1->w[0][0].imag());
}
/*! \brief C++ porting of RABAS
@@ -535,56 +528,50 @@ void rabas(
std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
std::complex<double> uim = std::complex<double>(0.0, 1.0);
double two_pi = 4.0 * acos(0.0);
- for (int i80 = 1; i80 <= nsph; i80++) {
- if(c1->iog[i80 - 1] >= i80) {
- tqse[0][i80 - 1] = 0.0;
- tqse[1][i80 - 1] = 0.0;
- tqspe[0][i80 - 1] = cc0;
- tqspe[1][i80 - 1] = cc0;
- tqss[0][i80 - 1] = 0.0;
- tqss[1][i80 - 1] = 0.0;
- tqsps[0][i80 - 1] = cc0;
- tqsps[1][i80 - 1] = cc0;
- for (int l70 = 1; l70 <= li; l70++) {
- double fl = 1.0 * l70 + l70 + 1;
- std::complex<double> rm = 1.0 / c1->rmi[l70 - 1][i80 - 1];
+ for (int i80 = 0; i80 < nsph; i80++) {
+ int i = i80 + 1;
+ if(c1->iog[i80] >= i) {
+ tqse[0][i80] = 0.0;
+ tqse[1][i80] = 0.0;
+ tqspe[0][i80] = cc0;
+ tqspe[1][i80] = cc0;
+ tqss[0][i80] = 0.0;
+ tqss[1][i80] = 0.0;
+ tqsps[0][i80] = cc0;
+ tqsps[1][i80] = cc0;
+ for (int l70 = 0; l70 < li; l70++) {
+ int l = l70 + 1;
+ double fl = 1.0 * (l + l + 1);
+ std::complex<double> rm = 1.0 / c1->rmi[l70][i80];
double rmm = (rm * dconjg(rm)).real();
- std::complex<double> re = 1.0 / c1->rei[l70 - 1][i80 - 1];
+ std::complex<double> re = 1.0 / c1->rei[l70][i80];
double rem = (re * dconjg(re)).real();
if (inpol == 0) {
std::complex<double> pce = ((rm + re) * uim) * fl;
std::complex<double> pcs = ((rmm + rem) * fl) * uim;
- tqspe[0][i80 - 1] -= pce;
- tqspe[1][i80 - 1] += pce;
- tqsps[0][i80 - 1] -= pcs;
- tqsps[1][i80 - 1] += pcs;
+ tqspe[0][i80] -= pce;
+ tqspe[1][i80] += pce;
+ tqsps[0][i80] -= pcs;
+ tqsps[1][i80] += pcs;
} else {
double ce = (rm + re).real() * fl;
double cs = (rmm + rem) * fl;
- tqse[0][i80 - 1] -= ce;
- tqse[1][i80 - 1] += ce;
- tqss[0][i80 - 1] -= cs;
- tqss[1][i80 - 1] += cs;
+ tqse[0][i80] -= ce;
+ tqse[1][i80] += ce;
+ tqss[0][i80] -= cs;
+ tqss[1][i80] += cs;
}
}
if (inpol == 0) {
- tqspe[0][i80 - 1] *= two_pi;
- tqspe[1][i80 - 1] *= two_pi;
- tqsps[0][i80 - 1] *= two_pi;
- tqsps[1][i80 - 1] *= two_pi;
- //printf("DEBUG: TQSPE(1, %d ) = ( %lE, %lE)\n", i80, tqspe[0][i80 - 1].real(), tqspe[0][i80 - 1].imag());
- //printf("DEBUG: TQSPE(2, %d ) = ( %lE, %lE)\n", i80, tqspe[1][i80 - 1].real(), tqspe[1][i80 - 1].imag());
- //printf("DEBUG: TQSPS(1, %d ) = ( %lE, %lE)\n", i80, tqsps[0][i80 - 1].real(), tqsps[0][i80 - 1].imag());
- //printf("DEBUG: TQSPS(2, %d ) = ( %lE, %lE)\n", i80, tqsps[1][i80 - 1].real(), tqsps[1][i80 - 1].imag());
+ tqspe[0][i80] *= two_pi;
+ tqspe[1][i80] *= two_pi;
+ tqsps[0][i80] *= two_pi;
+ tqsps[1][i80] *= two_pi;
} else {
- tqse[0][i80 - 1] *= two_pi;
- tqse[1][i80 - 1] *= two_pi;
- tqss[0][i80 - 1] *= two_pi;
- tqss[1][i80 - 1] *= two_pi;
- //printf("DEBUG: TQSE(1, %d ) = %lE)\n", i80, tqse[0][i80 - 1]);
- //printf("DEBUG: TQSE(2, %d ) = %lE)\n", i80, tqse[1][i80 - 1]);
- //printf("DEBUG: TQSS(1, %d ) = %lE)\n", i80, tqss[0][i80 - 1]);
- //printf("DEBUG: TQSS(2, %d ) = %lE)\n", i80, tqss[1][i80 - 1]);
+ tqse[0][i80] *= two_pi;
+ tqse[1][i80] *= two_pi;
+ tqss[0][i80] *= two_pi;
+ tqss[1][i80] *= two_pi;
}
}
}
@@ -618,8 +605,8 @@ void rbf(int n, double x, int &nm, double sj[]) {
if (a0 < 0.0) a0 *= -1.0;
nm = n;
if (a0 < 1.0e-60) {
- for (int k = 2; k <= n + 1; k++)
- sj[k - 1] = 0.0;
+ for (int k = 1; k <= n; k++)
+ sj[k] = 0.0;
sj[0] = 1.0;
return;
}
@@ -646,10 +633,8 @@ void rbf(int n, double x, int &nm, double sj[]) {
f1 = f;
}
double cs;
- double abs_sa = sa;
- if (abs_sa < 0.0) abs_sa *= -1.0;
- double abs_sb = sb;
- if (abs_sb < 0.0) abs_sb *= -1.0;
+ double abs_sa = (sa < 0.0) ? -sa : sa;
+ double abs_sb = (sb < 0.0) ? -sb : sb;
if (abs_sa > abs_sb) cs = sa / f;
else cs = sb / f0;
for (int k = 0; k <= nm; k++) {
@@ -678,7 +663,7 @@ void rkc(
std::complex<double> cy4, yy, c14, c21, c22, c23, c24;
double half_step = 0.5 * step;
double cl = 1.0 * lpo * (lpo - 1);
- for (int ipnt60 = 1; ipnt60 <= npntmo; ipnt60++) {
+ for (int ipnt60 = 0; ipnt60 < npntmo; ipnt60++) {
cy1 = cl / (x * x) - dcc;
cdy1 = -2.0 / x;
c11 = (cy1 * y1 + cdy1 * dy1) * step;
@@ -698,7 +683,7 @@ void rkc(
c21 = (cy1 * y2 + cdy1 * dy2) * step;
yc2 = y2 + dy2 * half_step;
c22 = (cy23 * yc2 + cdy23 * (dy2 + 0.5 * c21)) * step;
- c23= (cy23 * (yc2 + 0.25 * c21 * step) + cdy23 * (dy2 + 0.5 * c22)) * step;
+ c23 = (cy23 * (yc2 + 0.25 * c21 * step) + cdy23 * (dy2 + 0.5 * c22)) * step;
yy = y2 + dy2 * step;
c24 = (cy4 * (yc2 + 0.5 * c22 * step) + cdy4 * (dy2 + c23)) * step;
y2 = yy + (c21 + c22 + c23) * step / 6.0;
@@ -728,36 +713,38 @@ void rkt(
std::complex<double> cy4, cdy4, yy, c14, c21, c22, c23, c24;
double half_step = 0.5 * step;
double cl = 1.0 * lpo * (lpo - 1);
- for (int ipnt60 = 1; ipnt60 <= npntmo; ipnt60++) {
- int jpnt = ipnt60 + ipnt60 - 1;
- cy1 = cl / (x * x) - c2->ris[jpnt - 1];
+ for (int ipnt60 = 0; ipnt60 < npntmo; ipnt60++) {
+ int ipnt = ipnt60 + 1;
+ int jpnt = ipnt + ipnt - 1;
+ int jpnt60 = jpnt - 1;
+ cy1 = cl / (x * x) - c2->ris[jpnt60];
cdy1 = -2.0 / x;
c11 = (cy1 * y1 + cdy1 * dy1) * step;
double xh = x + half_step;
int jpntpo = jpnt + 1;
- cy23 = cl / (xh * xh) - c2->ris[jpntpo - 1];
+ cy23 = cl / (xh * xh) - c2->ris[jpnt];
cdy23 = -2.0 / xh;
yc2 = y1 + dy1 * half_step;
c12 = (cy23 * yc2 + cdy23 * (dy1 + 0.5 * c11)) * step;
c13= (cy23 * (yc2 + 0.25 * c11 *step) + cdy23 * (dy1 + 0.5 * c12)) * step;
double xn = x + step;
- int jpntpt = jpnt + 2;
- cy4 = cl / (xn * xn) - c2->ris[jpntpt - 1];
+ //int jpntpt = jpnt + 2;
+ cy4 = cl / (xn * xn) - c2->ris[jpntpo];
cdy4 = -2.0 / xn;
yy = y1 + dy1 * step;
c14 = (cy4 * (yy + 0.5 * c12 * step) + cdy4 * (dy1 + c13)) * step;
y1= yy + (c11 + c12 + c13) * step / 6.0;
dy1 += (0.5 * c11 + c12 + c13 + 0.5 * c14) /3.0;
- cy1 -= cdy1 * c2->dlri[jpnt - 1];
- cdy1 += 2.0 * c2->dlri[jpnt - 1];
+ cy1 -= cdy1 * c2->dlri[jpnt60];
+ cdy1 += 2.0 * c2->dlri[jpnt60];
c21 = (cy1 * y2 + cdy1 * dy2) * step;
- cy23 -= cdy23 * c2->dlri[jpntpo - 1];
- cdy23 += 2.0 * c2->dlri[jpntpo - 1];
+ cy23 -= cdy23 * c2->dlri[jpnt];
+ cdy23 += 2.0 * c2->dlri[jpnt];
yc2 = y2 + dy2 * half_step;
c22 = (cy23 * yc2 + cdy23 * (dy2 + 0.5 * c21)) * step;
c23 = (cy23 * (yc2 + 0.25 * c21 * step) + cdy23 * (dy2 + 0.5 * c22)) * step;
- cy4 -= cdy4 * c2->dlri[jpntpt - 1];
- cdy4 += 2.0 * c2->dlri[jpntpt - 1];
+ cy4 -= cdy4 * c2->dlri[jpntpo];
+ cdy4 += 2.0 * c2->dlri[jpntpo];
yy = y2 + dy2 * step;
c24 = (cy4 * (yc2 + 0.5 * c22 * step) + cdy4 * (dy2 + c23)) * step;
y2 = yy + (c21 + c22 + c23) * step / 6.0;
@@ -830,23 +817,25 @@ void rnf(int n, double x, int &nm, double sy[]) {
*/
void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
std::complex<double> sum21, rm, re, csam;
- std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
- double exdc = exri * exri;
+ const std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
+ const 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);
tfsas = cc0;
- for (int i12 = 1; i12 <= nsph; i12++) {
- int iogi = c1->iog[i12 - 1];
- if (iogi >= i12) {
+ for (int i12 = 0; i12 < nsph; i12++) {
+ int i = i12 + 1;
+ int iogi = c1->iog[i12];
+ if (iogi >= i) {
double sums = 0.0;
std::complex<double> sum21 = cc0;
- for (int l10 = 1; l10 <= lm; l10++) {
- double fl = 1.0 + l10 + l10;
- rm = 1.0 / c1->rmi[l10 - 1][i12 - 1];
- re = 1.0 / c1->rei[l10 - 1][i12 - 1];
- std::complex<double> rm_cnjg = std::complex<double>(rm.real(), -rm.imag());
- std::complex<double> re_cnjg = std::complex<double>(re.real(), -re.imag());
+ for (int l10 = 0; l10 < lm; l10++) {
+ int l = l10 + 1;
+ double fl = 1.0 + l + l;
+ rm = 1.0 / c1->rmi[l10][i12];
+ re = 1.0 / c1->rei[l10][i12];
+ std::complex<double> rm_cnjg = dconjg(rm);
+ std::complex<double> re_cnjg = dconjg(re);
sums += (rm_cnjg * rm + re_cnjg * re).real() * fl;
sum21 += (rm + re) * fl;
}
@@ -854,15 +843,14 @@ void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri
double scasec = cccs * sums;
double extsec = -cccs * sum21.real();
double abssec = extsec - scasec;
- c1->sscs[i12 - 1] = scasec;
- c1->sexs[i12 - 1] = extsec;
- c1->sabs[i12 - 1] = abssec;
- double gcss = c1-> gcsv[i12 - 1];
- c1->sqscs[i12 - 1] = scasec / gcss;
- c1->sqexs[i12 - 1] = extsec / gcss;
- c1->sqabs[i12 - 1] = abssec / gcss;
- c1->fsas[i12 - 1] = sum21 * csam;
- //printf("DEBUG: FSAS( %d ) = (%lE,%lE)\n", i12, c1->fsas[i12 - 1].real(), c1->fsas[i12 - 1].imag());
+ c1->sscs[i12] = scasec;
+ c1->sexs[i12] = extsec;
+ c1->sabs[i12] = abssec;
+ double gcss = c1->gcsv[i12];
+ c1->sqscs[i12] = scasec / gcss;
+ c1->sqexs[i12] = extsec / gcss;
+ c1->sqabs[i12] = abssec / gcss;
+ c1->fsas[i12] = sum21 * csam;
}
tfsas += c1->fsas[iogi - 1];
}
@@ -878,61 +866,52 @@ void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri
*/
void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
std::complex<double> s11, s21, s12, s22, rm, re, csam, cc;
- std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
+ const std::complex<double> cc0(0.0, 0.0);
double ccs = 1.0 / (vk * vk);
csam = -(ccs / (exri * vk)) * std::complex<double>(0.0, 0.5);
- double pigfsq = 64.0 * acos(0.0) * acos(0.0);
+ const double pigfsq = 64.0 * acos(0.0) * acos(0.0);
double cfsq = 4.0 / (pigfsq * ccs * ccs);
int nlmm = lm * (lm + 2);
- //printf("DEBUG: in SSCR2 W(1,1) = (%lE,%lE)\n", c1->w[0][0].real(), c1->w[0][0].imag());
- for (int i14 = 1; i14 <= nsph; i14++) {
- int iogi = c1->iog[i14 - 1];
- if (iogi >= i14) {
+ for (int i14 = 0; i14 < nsph; i14++) {
+ int i = i14 + 1;
+ int iogi = c1->iog[i14];
+ if (iogi >= i) {
int k = 0;
s11 = cc0;
s21 = cc0;
s12 = cc0;
s22 = cc0;
- for (int l10 = 1; l10 <= lm; l10++) {
- rm = 1.0 / c1->rmi[l10 - 1][i14 - 1];
- re = 1.0 / c1->rei[l10 - 1][i14 - 1];
- //printf("DEBUG: rm = (%lE,%lE)\n", rm.real(), rm.imag());
- //printf("DEBUG: re = (%lE,%lE)\n", re.real(), re.imag());
- int ltpo = l10 + l10 + 1;
- for (int im10 = 1; im10 <= ltpo; im10++) {
+ for (int l10 = 0; l10 < lm; l10++) {
+ int l = l10 + 1;
+ rm = 1.0 / c1->rmi[l10][i14];
+ re = 1.0 / c1->rei[l10][i14];
+ int ltpo = l + l + 1;
+ for (int im10 = 0; im10 < ltpo; im10++) {
k += 1;
int ke = k + nlmm;
- //printf("DEBUG: W( %d, 3) = (%lE,%lE)\n", k, c1->w[k - 1][2].real(), c1->w[k - 1][2].imag());
- //printf("DEBUG: W( %d, 1) = (%lE,%lE)\n", k, c1->w[k - 1][0].real(), c1->w[k - 1][0].imag());
- //printf("DEBUG: W( %d, 3) = (%lE,%lE)\n", ke, c1->w[ke - 1][2].real(), c1->w[ke - 1][2].imag());
- //printf("DEBUG: W( %d, 1) = (%lE,%lE)\n", ke, c1->w[ke - 1][0].real(), c1->w[ke - 1][0].imag());
s11 = s11 - c1->w[k - 1][2] * c1->w[k - 1][0] * rm - c1->w[ke - 1][2] * c1->w[ke - 1][0] * re;
- //printf("DEBUG: W( %d, 4) = (%lE,%lE)\n", k, c1->w[k - 1][3].real(), c1->w[k - 1][3].imag());
- //printf("DEBUG: W( %d, 4) = (%lE,%lE)\n", ke, c1->w[ke - 1][3].real(), c1->w[ke - 1][3].imag());
s21 = s21 - c1->w[k - 1][3] * c1->w[k - 1][0] * rm - c1->w[ke - 1][3] * c1->w[ke - 1][0] * re;
- //printf("DEBUG: W( %d, 2) = (%lE,%lE)\n", k, c1->w[k - 1][1].real(), c1->w[k - 1][1].imag());
- //printf("DEBUG: W( %d, 2) = (%lE,%lE)\n", ke, c1->w[ke - 1][1].real(), c1->w[ke - 1][1].imag());
s12 = s12 - c1->w[k - 1][2] * c1->w[k - 1][1] * rm - c1->w[ke - 1][2] * c1->w[ke - 1][1] * re;
s22 = s22 - c1->w[k - 1][3] * c1->w[k - 1][1] * rm - c1->w[ke - 1][3] * c1->w[ke - 1][1] * re;
}
}
- c1->sas[i14 - 1][0][0] = s11 * csam;
- c1->sas[i14 - 1][1][0] = s21 * csam;
- c1->sas[i14 - 1][0][1] = s12 * csam;
- c1->sas[i14 - 1][1][1] = s22 * csam;
+ c1->sas[i14][0][0] = s11 * csam;
+ c1->sas[i14][1][0] = s21 * csam;
+ c1->sas[i14][0][1] = s12 * csam;
+ c1->sas[i14][1][1] = s22 * csam;
}
} // loop i14
- for (int i24 = 1; i24 <= nsph; i24++) {
- int iogi = c1->iog[i24 - 1];
- if (iogi >= i24) {
+ for (int i24 = 0; i24 < nsph; i24++) {
+ int i = i24 + 1;
+ int iogi = c1->iog[i24];
+ if (iogi >= i) {
int j = 0;
for (int ipo1 = 0; ipo1 < 2; ipo1++) {
for (int jpo1 = 0; jpo1 < 2; jpo1++) {
- std::complex<double> cc = dconjg(c1->sas[i24 - 1][jpo1][ipo1]);
+ std::complex<double> cc = dconjg(c1->sas[i24][jpo1][ipo1]);
for (int ipo2 = 0; ipo2 < 2; ipo2++) {
for (int jpo2 = 0; jpo2 < 2; jpo2++) {
- j += 1;
- c1->vints[i24 - 1][j - 1] = c1->sas[i24 - 1][jpo2][ipo2] * cc * cfsq;
+ c1->vints[i24][j++] = c1->sas[i24][jpo2][ipo2] * cc * cfsq;
}
}
}
@@ -954,25 +933,20 @@ void sphar(
double cosrth, double sinrth, double cosrph, double sinrph,
int ll, std::complex<double> *ylm
) {
- int rmp_size = ll;
- int plegn_size = (ll + 1) * ll / 2 + ll + 1;
+ const int rmp_size = ll;
+ const int plegn_size = (ll + 1) * ll / 2 + ll + 1;
double sinrmp[rmp_size], cosrmp[rmp_size], plegn[plegn_size];
double four_pi = 8.0 * acos(0.0);
double pi4irs = 1.0 / sqrt(four_pi);
double x = cosrth;
double y = sinrth;
- //printf("DEBUG: X = %lE\n", x);
if (y < 0.0) y *= -1.0;
- //printf("DEBUG: Y = %lE\n", y);
double cllmo = 3.0;
double cll = 1.5;
double ytol = y;
plegn[0] = 1.0;
- //printf("DEBUG: PLEGN( %d ) = %lE\n", 1, plegn[0]);
plegn[1] = x * sqrt(cllmo);
- //printf("DEBUG: PLEGN( %d ) = %lE\n", 2, plegn[1]);
plegn[2] = ytol * sqrt(cll);
- //printf("DEBUG: PLEGN( %d ) = %lE\n", 3, plegn[2]);
sinrmp[0] = sinrph;
cosrmp[0] = cosrph;
if (ll >= 2) {
@@ -992,18 +966,15 @@ void sphar(
double cdm = 1.0 * ls * (ltmo - 2);
plegn[mpopk - 1] = plegn[mpopk - l20 - 1] * x * sqrt(cn / cd) -
plegn[mpopk - ltmo - 1] * sqrt(cnm / cdm);
- //printf("DEBUG: PLEGN( %d ) = %lE\n", mpopk, plegn[mpopk - 1]);
}
int lpk = l20 + k;
double cltpo = 1.0 * ltpo;
plegn[lpk - 1] = plegn[k - 1] * x * sqrt(cltpo);
- //printf("DEBUG: PLEGN( %d ) = %lE\n", lpk, plegn[lpk - 1]);
k = lpk + 1;
double clt = 1.0 * (ltpo - 1);
cll *= (cltpo / clt);
ytol *= y;
plegn[k - 1] = ytol * sqrt(cll);
- //printf("DEBUG: PLEGN( %d ) = %lE\n", k, plegn[k - 1]);
sinrmp[l20 - 1] = sinrph * cosrmp[lmo - 1] + cosrph * sinrmp[lmo - 1];
cosrmp[l20 - 1] = cosrph * cosrmp[lmo - 1] - sinrph * sinrmp[lmo - 1];
} // end l20 loop
@@ -1018,7 +989,6 @@ label40:
l0y = k + 1;
l0p = k / 2 + 1;
ylm[l0y - 1] = pi4irs * plegn[l0p - 1];
- //printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", l0y, ylm[l0y - 1].real(), ylm[l0y - 1].imag());
goto label45;
label44:
lmp = l0p + m;
@@ -1026,10 +996,8 @@ label44:
lmy = l0y + m;
ylm[lmy - 1] = save * std::complex<double>(cosrmp[m - 1], sinrmp[m - 1]);
if (m % 2 != 0) ylm[lmy - 1] *= -1.0;
- //printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", lmy, ylm[lmy - 1].real(), ylm[lmy - 1].imag());
lmy = l0y - m;
ylm[lmy - 1] = save * std::complex<double>(cosrmp[m - 1], -sinrmp[m - 1]);
- //printf("DEBUG: YLM( %d ) = (%lE,%lE)\n", lmy, ylm[lmy - 1].real(), ylm[lmy - 1].imag());
label45:
if (m >= l) goto label47;
m += 1;
@@ -1046,42 +1014,40 @@ label47:
* \param zpv: `double ****`
*/
void thdps(int lm, double ****zpv) {
- //double szpv = 0.0;
- for (int l10 = 0; l10 < lm; l10++) {
- for (int ilmp = 0; ilmp < 3; ilmp++) {
- zpv[l10][ilmp][0][0] = 0.0;
- zpv[l10][ilmp][0][1] = 0.0;
- zpv[l10][ilmp][1][0] = 0.0;
- zpv[l10][ilmp][1][1] = 0.0;
- }
- }
- for (int l15 = 1; l15 <= lm; l15++) {
- double xd = 1.0 * l15 * (l15 + 1);
+ //for (int l10 = 0; l10 < lm; l10++) { // 0-init, can be omitted
+ // for (int ilmp = 0; ilmp < 3; ilmp++) {
+ // zpv[l10][ilmp][0][0] = 0.0;
+ // zpv[l10][ilmp][0][1] = 0.0;
+ // zpv[l10][ilmp][1][0] = 0.0;
+ // zpv[l10][ilmp][1][1] = 0.0;
+ // }
+ //}
+ for (int l15 = 0; l15 < lm; l15++) {
+ int l = l15 + 1;
+ double xd = 1.0 * l * (l + 1);
double zp = -1.0 / sqrt(xd);
- zpv[l15 - 1][1][0][1] = zp;
- zpv[l15 - 1][1][1][0] = zp;
- //szpv += 2.0 * zp;
+ zpv[l15][1][0][1] = zp;
+ zpv[l15][1][1][0] = zp;
}
if (lm != 1) {
- for (int l20 = 2; l20 <= lm; l20++) {
- double xn = 1.0 * (l20 - 1) * (l20 + 1);
- double xd = 1.0 * l20 * (l20 + l20 + 1);
+ for (int l20 = 1; l20 < lm; l20++) {
+ int l = l20 + 1;
+ double xn = 1.0 * (l - 1) * (l + 1);
+ double xd = 1.0 * l * (l + l + 1);
double zp = sqrt(xn / xd);
- zpv[l20 - 1][0][0][0] = zp;
- zpv[l20 - 1][0][1][1] = zp;
- //szpv += 2.0 * zp;
+ zpv[l20][0][0][0] = zp;
+ zpv[l20][0][1][1] = zp;
}
int lmmo = lm - 1;
- for (int l25 = 1; l25 <= lmmo; l25++) {
- double xn = 1.0 * l25 * (l25 + 2);
- double xd = (l25 + 1) * (l25 + l25 + 1);
+ for (int l25 = 0; l25 < lmmo; l25++) {
+ int l = l25 + 1;
+ double xn = 1.0 * l * (l + 2);
+ double xd = (l + 1) * (l + l + 1);
double zp = -1.0 * sqrt(xn / xd);
- zpv[l25 - 1][2][0][0] = zp;
- zpv[l25 - 1][2][1][1] = zp;
- //szpv += 2.0 * zp;
+ zpv[l25][2][0][0] = zp;
+ zpv[l25][2][1][1] = zp;
}
}
- //printf("DEBUG: szpv = %lE\n", szpv);
}
/*! \brief C++ porting of UPVMP
@@ -1109,10 +1075,6 @@ void upvmp(
sint = sin(th);
cosp = cos(ph);
sinp = sin(ph);
- //printf("DEBUG: cost = %lE\n", cost);
- //printf("DEBUG: sint = %lE\n", sint);
- //printf("DEBUG: cosp = %lE\n", cosp);
- //printf("DEBUG: sinp = %lE\n", sinp);
u[0] = cosp * sint;
u[1] = sinp * sint;
u[2] = cost;
@@ -1161,8 +1123,7 @@ void upvsp(
double small = 1.0e-6;
isq = 0;
scand = u[0] * us[0] + u[1] * us[1] + u[2] * us[2];
- double abs_scand = scand - 1.0;
- if (abs_scand < 0.0) abs_scand *= -1.0;
+ double abs_scand = (scand >= 1.0) ? scand - 1.0 : 1.0 - scand;
if (abs_scand >= small) {
abs_scand = scand + 1.0;
if (abs_scand < 0.0) abs_scand *= -1.0;
@@ -1172,8 +1133,7 @@ void upvsp(
duk[1] = u[1] - us[1];
duk[2] = u[2] - us[2];
ibf = 0;
- } else {
- // label 15
+ } else { // label 15
scand = 180.0;
duk[0] = 2.0 * u[0];
duk[1] = 2.0 * u[1];
@@ -1186,8 +1146,7 @@ void upvsp(
uns[1] = -unsmp[1];
uns[2] = -unsmp[2];
}
- } else {
- // label 10
+ } else { // label 10
scand = 0.0;
duk[0] = 0.0;
duk[1] = 0.0;
@@ -1245,9 +1204,9 @@ void wmamp(
int lm, int idot, int nsph, double *arg, double *u, double *up,
double *un, C1 *c1
) {
- int ylm_size = (lm + 1) * (lm + 1) + 1;
+ const int ylm_size = (lm + 1) * (lm + 1) + 1;
std::complex<double> *ylm = new std::complex<double>[ylm_size];
- int nlmp = lm * (lm + 2) + 2;
+ const int nlmp = lm * (lm + 2) + 2;
ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
if (idot != 0) {
if (idot != 1) {
@@ -1301,27 +1260,27 @@ void wmasp(
double *ups, double *uns, int isq, int ibf, int inpol, int lm, int idot,
int nsph, double *argi, double *args, C1 *c1
) {
- int ylm_size = (lm + 1) * (lm + 1) + 1;
+ const int ylm_size = (lm + 1) * (lm + 1) + 1;
std::complex<double> *ylm = new std::complex<double>[ylm_size];
- int nlmp = lm * (lm + 2) + 2;
+ const int nlmp = lm * (lm + 2) + 2;
ylm[nlmp - 1] = std::complex<double>(0.0, 0.0);
if (idot != 0) {
if (idot != 1) {
- for (int n40 = 1; n40 <= nsph; n40++) {
- argi[n40 - 1] = u[0] * c1->rxx[n40 - 1] + u[1] * c1->ryy[n40 - 1] + u[2] * c1->rzz[n40 - 1];
+ for (int n40 = 0; n40 < nsph; n40++) {
+ argi[n40] = u[0] * c1->rxx[n40] + u[1] * c1->ryy[n40] + u[2] * c1->rzz[n40];
if (ibf != 0) {
- args[n40 - 1] = argi[n40 - 1] * ibf;
+ args[n40] = argi[n40] * ibf;
} else {
- args[n40 - 1] = -1.0 * (us[0] * c1->rxx[n40 - 1] + us[1] * c1->ryy[n40 - 1] + us[2] * c1->rzz[n40 - 1]);
+ args[n40] = -1.0 * (us[0] * c1->rxx[n40] + us[1] * c1->ryy[n40] + us[2] * c1->rzz[n40]);
}
}
} else { // label 50
- for (int n60 = 1; n60 <= nsph; n60++) {
- argi[n60 - 1] = cost * c1->rzz[n60 - 1];
+ for (int n60 = 0; n60 < nsph; n60++) {
+ argi[n60] = cost * c1->rzz[n60];
if (ibf != 0) {
- args[n60 - 1] = argi[n60 - 1] * ibf;
+ args[n60] = argi[n60] * ibf;
} else {
- args[n60 - 1] = -costs * c1->rzz[n60 - 1];
+ args[n60] = -costs * c1->rzz[n60];
}
}
}
@@ -1337,7 +1296,6 @@ void wmasp(
delete[] ylm;
}
-
/*! \brief C++ porting of DME
*
* \param li: `int`
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 85005091273ac3fa34f4a0ff96eb11919e3216a8..7565ed594e2f9caca942a4d04d350cf620b4bab8 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -152,17 +152,18 @@ void sphere() {
const double half_pi = acos(0.0);
const double pi = 2.0 * half_pi;
double gcs = 0.0;
- for (int i116 = 1; i116 <= nsph; i116++) {
- int iogi = c1->iog[i116 - 1];
- if (iogi >= i116) {
- double gcss = pi * c1->ros[i116 - 1] * c1->ros[i116 - 1];
- c1->gcsv[i116 - 1] = gcss;
- int nsh = c1->nshl[i116 - 1];
- for (int j115 = 1; j115 <= nsh; j115++) {
- c1->rc[i116 - 1][j115 - 1] = sconf->rcf[i116 - 1][j115 - 1] * c1->ros[i116 - 1];
+ for (int i116 = 0; i116 < nsph; i116++) {
+ int i = i116 + 1;
+ int iogi = c1->iog[i116];
+ if (iogi >= i) {
+ double gcss = pi * c1->ros[i116] * c1->ros[i116];
+ c1->gcsv[i116] = gcss;
+ int nsh = c1->nshl[i116];
+ for (int j115 = 0; j115 < nsh; j115++) {
+ c1->rc[i116][j115] = sconf->rcf[i116][j115] * c1->ros[i116];
}
}
- gcs += c1->gcsv[iogi - 1];
+ gcs += c1->gcsv[iogi];
}
double ****zpv = new double***[gconf->l_max]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
for (int zi = 0; zi < gconf->l_max; zi++) {
@@ -170,7 +171,7 @@ void sphere() {
for (int zj = 0; zj < 3; zj++) {
zpv[zi][zj] = new double*[2];
for (int zk = 0; zk < 2; zk++) {
- zpv[zi][zj][zk] = new double[2];
+ zpv[zi][zj][zk] = new double[2]();
}
}
}
@@ -203,9 +204,10 @@ void sphere() {
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
fprintf(output, " \n");
}
- for (int jxi488 = 1; jxi488 <= sconf->number_of_scales; jxi488++) {
- fprintf(output, "========== JXI =%3d ====================\n", jxi488);
- double xi = sconf->scale_vec[jxi488 - 1];
+ for (int jxi488 = 0; jxi488 < sconf->number_of_scales; jxi488++) {
+ int jxi = jxi488 + 1;
+ fprintf(output, "========== JXI =%3d ====================\n", jxi);
+ double xi = sconf->scale_vec[jxi488];
if (sconf->idfc >= 0) {
vk = xi * wn;
vkarg = vk;
@@ -216,24 +218,25 @@ void sphere() {
fprintf(output, " XI=%15.7lE\n", xi);
}
tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
- for (int i132 = 1; i132 <= nsph; i132++) {
- int iogi = sconf->iog_vec[i132 - 1];
- if (iogi != i132) {
- for (int l123 = 1; l123 <= gconf->l_max; l123++) {
- c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];
- c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];
+ for (int i132 = 0; i132 < nsph; i132++) {
+ int i = i132 + 1;
+ int iogi = sconf->iog_vec[i132];
+ if (iogi != i) {
+ for (int l123 = 0; l123 < gconf->l_max; l123++) {
+ c1->rmi[l123][i132] = c1->rmi[l123][iogi - 1];
+ c1->rei[l123][i132] = c1->rei[l123][iogi - 1];
}
continue; // i132
}
- int nsh = sconf->nshl_vec[i132 - 1];
+ int nsh = sconf->nshl_vec[i132];
int ici = (nsh + 1) / 2;
if (sconf->idfc == 0) {
for (int ic = 0; ic < ici; ic++)
- c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][0]; // WARNING: IDFC=0 is not tested!
+ c2->dc0[ic] = sconf->dc0_matrix[ic][i132][0]; // WARNING: IDFC=0 is not tested!
} else { // IDFC != 0
- if (jxi488 == 1) {
+ if (jxi == 1) {
for (int ic = 0; ic < ici; ic++) {
- c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][jxi488 - 1];
+ c2->dc0[ic] = sconf->dc0_matrix[ic][i132][jxi488];
}
}
}
@@ -241,7 +244,7 @@ void sphere() {
int jer = 0;
int lcalc = 0;
dme(
- gconf->l_max, i132, gconf->npnt, gconf->npntts, vkarg,
+ gconf->l_max, i, gconf->npnt, gconf->npntts, vkarg,
sconf->exdc, exri, c1, c2, jer, lcalc, arg
);
if (jer != 0) {
@@ -252,7 +255,7 @@ void sphere() {
return;
}
} // i132
- if (sconf->idfc >= 0 and nsph == 1 and jxi488 == gconf->jwtm) {
+ if (sconf->idfc >= 0 and nsph == 1 and jxi == gconf->jwtm) {
// This is the condition that writes the transition matrix to output.
int is = 1111;
fstream ttms;
@@ -273,50 +276,51 @@ void sphere() {
} else { // Failed to open output file. Should never happen.
printf("ERROR: could not open TTMS file.\n");
tppoan.close();
- fclose(output);
}
}
double cs0 = 0.25 * vk * vk * vk / half_pi;
+ //printf("DEBUG: cs0 = %lE\n", cs0);
sscr0(tfsas, nsph, gconf->l_max, vk, exri, c1);
printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag());
double sqk = vk * vk * sconf->exdc;
aps(zpv, gconf->l_max, nsph, c1, sqk, gaps);
rabas(gconf->in_pol, gconf->l_max, nsph, c1, tqse, tqspe, tqss, tqsps);
- for (int i170 = 1; i170 <= nsph; i170++) {
- if (c1->iog[i170 - 1] >= i170) {
- double albeds = c1->sscs[i170 - 1] / c1->sexs[i170 - 1];
- c1->sqscs[i170 - 1] *= sqsfi;
- c1->sqabs[i170 - 1] *= sqsfi;
- c1->sqexs[i170 - 1] *= sqsfi;
- fprintf(output, " SPHERE %2d\n", i170);
- if (c1->nshl[i170 - 1] != 1) {
- fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i170 - 1]);
+ for (int i170 = 0; i170 < nsph; i170++) {
+ int i = i170 + 1;
+ if (c1->iog[i170] >= i) {
+ double albeds = c1->sscs[i170] / c1->sexs[i170];
+ c1->sqscs[i170] *= sqsfi;
+ c1->sqabs[i170] *= sqsfi;
+ c1->sqexs[i170] *= sqsfi;
+ fprintf(output, " SPHERE %2d\n", i);
+ if (c1->nshl[i170] != 1) {
+ fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i170]);
} else {
fprintf(
output,
" SIZE=%15.7lE, REFRACTIVE INDEX=(%15.7lE,%15.7lE)\n",
- c2->vsz[i170 -1],
- c2->vkt[i170 - 1].real(),
- c2->vkt[i170 - 1].imag()
+ c2->vsz[i170],
+ c2->vkt[i170].real(),
+ c2->vkt[i170].imag()
);
}
fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
fprintf(
output,
" %14.7lE%15.7lE%15.7lE%15.7lE\n",
- c1->sscs[i170 - 1], c1->sabs[i170 - 1],
- c1->sexs[i170 - 1], albeds
+ c1->sscs[i170], c1->sabs[i170],
+ c1->sexs[i170], albeds
);
fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
fprintf(
output,
" %14.7lE%15.7lE%15.7lE\n",
- c1->sqscs[i170 - 1], c1->sqabs[i170 - 1],
- c1->sqexs[i170 - 1]
+ c1->sqscs[i170], c1->sqabs[i170],
+ c1->sqexs[i170]
);
- fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i170 - 1].real(), c1->fsas[i170 - 1].imag());
- double csch = 2.0 * vk * sqsfi / c1->gcsv[i170 -1];
- s0 = c1->fsas[i170 - 1] * exri;
+ fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i170].real(), c1->fsas[i170].imag());
+ double csch = 2.0 * vk * sqsfi / c1->gcsv[i170];
+ s0 = c1->fsas[i170] * exri;
double qschu = csch * s0.imag();
double pschu = csch * s0.real();
double s0mag = cs0 * abs(s0);
@@ -325,32 +329,32 @@ void sphere() {
" QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
qschu, pschu, s0mag
);
- double rapr = c1->sexs[i170 - 1] - gaps[i170 - 1];
- double cosav = gaps[i170 - 1] / c1->sscs[i170 - 1];
+ double rapr = c1->sexs[i170] - gaps[i170];
+ double cosav = gaps[i170] / c1->sscs[i170];
fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
- fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 1, tqse[0][i170 - 1], tqss[0][i170 - 1]);
- fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170 - 1], tqss[1][i170 - 1]);
- tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170 - 1])), sizeof(double));
- tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170 - 1])), sizeof(double));
- double val = tqspe[0][i170 - 1].real();
+ fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 1, tqse[0][i170], tqss[0][i170]);
+ fprintf(output, " IPO=%2d, TQEk=%15.7lE, TQSk=%15.7lE\n", 2, tqse[1][i170], tqss[1][i170]);
+ tppoan.write(reinterpret_cast<char *>(&(tqse[0][i170])), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&(tqss[0][i170])), sizeof(double));
+ double val = tqspe[0][i170].real();
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqspe[0][i170 - 1].imag();
+ val = tqspe[0][i170].imag();
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[0][i170 - 1].real();
+ val = tqsps[0][i170].real();
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[0][i170 - 1].imag();
+ val = tqsps[0][i170].imag();
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- tppoan.write(reinterpret_cast<char *>(&(tqse[1][i170 - 1])), sizeof(double));
- tppoan.write(reinterpret_cast<char *>(&(tqss[1][i170 - 1])), sizeof(double));
- val = tqspe[1][i170 - 1].real();
+ tppoan.write(reinterpret_cast<char *>(&(tqse[1][i170])), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&(tqss[1][i170])), sizeof(double));
+ val = tqspe[1][i170].real();
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqspe[1][i170 - 1].imag();
+ val = tqspe[1][i170].imag();
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[1][i170 - 1].real();
+ val = tqsps[1][i170].real();
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- val = tqsps[1][i170 - 1].imag();
+ val = tqsps[1][i170].imag();
tppoan.write(reinterpret_cast<char *>(&val), sizeof(double));
- } // End if iog[i170 - 1] >= i170
+ } // End if iog[i170] >= i
} // i170 loop
if (nsph != 1) {
fprintf(output, " FSAT=(%15.7lE,%15.7lE)\n", tfsas.real(), tfsas.imag());
@@ -366,10 +370,12 @@ void sphere() {
);
}
th = th1;
- for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP parallelizable section
+ for (int jth486 = 0; jth486 < nth; jth486++) { // OpenMP parallelizable section
+ int jth = jth486 + 1;
ph = ph1;
- for (int jph484 = 1; jph484 <= nph; jph484++) {
- bool goto182 = (nk == 1) && (jxi488 > 1);
+ for (int jph484 = 0; jph484 < nph; jph484++) {
+ int jph = jph484 + 1;
+ bool goto182 = (nk == 1) && (jxi > 1);
if (!goto182) {
upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
}
@@ -385,7 +391,8 @@ void sphere() {
}
double thsl = ths1;
double phsph = 0.0;
- for (int jths482 = 1; jths482 <= nths; jths482++) {
+ for (int jths482 = 0; jths482 < nths; jths482++) {
+ int jths = jths482 + 1;
double ths = thsl;
int icspnv = 0;
if (gconf->meridional_type > 1) ths = th + thsca;
@@ -397,19 +404,20 @@ void sphere() {
if (phsph != 0.0) icspnv = 1;
}
double phs = phs1;
- for (int jphs480 = 1; jphs480 <= nphs; jphs480++) {
+ for (int jphs480 = 0; jphs480 < nphs; jphs480++) {
+ int jphs = jphs480 + 1;
if (gconf->meridional_type >= 1) {
phs = ph + phsph;
if (phs >= 360.0) phs -= 360.0;
}
- bool goto190 = (nks == 1) && ((jxi488 > 1) || (jth486 > 1) || (jph484 > 1));
+ bool goto190 = (nks == 1) && ((jxi > 1) || (jth > 1) || (jph > 1));
if (!goto190) {
upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
if (gconf->meridional_type >= 0) {
wmamp(2, costs, sints, cosps, sinps, gconf->in_pol, gconf->l_max, 0, nsph, args, us, upsmp, unsmp, c1);
}
}
- if (nkks != 0 || jxi488 == 1) {
+ if (nkks != 0 || jxi == 1) {
upvsp(u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns, duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp);
if (gconf->meridional_type < 0) {
wmasp(
@@ -418,9 +426,9 @@ void sphere() {
gconf->l_max, 0, nsph, argi, args, c1
);
}
- for (int i193 = 1; i193 <= 3; i193++) {
- un[i193 - 1] = unmp[i193 - 1];
- uns[i193 - 1] = unsmp[i193 - 1];
+ for (int i193 = 0; i193 < 3; i193++) {
+ un[i193] = unmp[i193];
+ uns[i193] = unsmp[i193];
}
}
if (gconf->meridional_type < 0) jw = 1;
@@ -438,7 +446,7 @@ void sphere() {
fprintf(
output,
"********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n",
- jth486, jph484, jths482, jphs480
+ jth, jph, jths, jphs
);
fprintf(
output,
@@ -455,25 +463,25 @@ void sphere() {
fprintf(output, " UN=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
}
sscr2(nsph, gconf->l_max, vk, exri, c1);
- for (int ns226 = 1; ns226 <= nsph; ns226++) {
- fprintf(output, " SPHERE %2d\n", ns226);
+ for (int ns226 = 0; ns226 < nsph; ns226++) {
+ int ns = ns226 + 1;
+ fprintf(output, " SPHERE %2d\n", ns);
fprintf(
output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
- c1->sas[ns226 - 1][0][0].real(), c1->sas[ns226 - 1][0][0].imag(),
- c1->sas[ns226 - 1][1][0].real(), c1->sas[ns226 - 1][1][0].imag()
+ c1->sas[ns226][0][0].real(), c1->sas[ns226][0][0].imag(),
+ c1->sas[ns226][1][0].real(), c1->sas[ns226][1][0].imag()
);
fprintf(
output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
- c1->sas[ns226 - 1][0][1].real(), c1->sas[ns226 - 1][0][1].imag(),
- c1->sas[ns226 - 1][1][1].real(), c1->sas[ns226 - 1][1][1].imag()
+ c1->sas[ns226][0][1].real(), c1->sas[ns226][0][1].imag(),
+ c1->sas[ns226][1][1].real(), c1->sas[ns226][1][1].imag()
);
- if (jths482 == 1 && jphs480 == 1)
+ if (jths == 1 && jphs == 1)
fprintf(
output, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n",
frx, fry, frz
);
- for (int i225 = 0; i225 < 16; i225++)
- vint[i225] = c1->vints[ns226 - 1][i225];
+ for (int i225 = 0; i225 < 16; i225++) vint[i225] = c1->vints[ns226][i225];
mmulc(vint, cmullr, cmul);
fprintf(output, " MULS\n ");
for (int imul = 0; imul < 4; imul++) {
@@ -519,6 +527,16 @@ void sphere() {
fclose(output);
delete c1;
delete c2;
+ for (int zi = gconf->l_max - 1; zi > -1; zi--) {
+ for (int zj = 0; zj < 3; zj++) {
+ for (int zk = 0; zk < 2; zk++) {
+ delete[] zpv[zi][zj][zk];
+ }
+ delete[] zpv[zi][zj];
+ }
+ delete[] zpv[zi];
+ }
+ delete[] zpv;
delete[] duk;
delete[] u;
delete[] us;
@@ -534,7 +552,7 @@ void sphere() {
delete[] argi;
delete[] args;
delete[] gaps;
- for (int i = 0; i < 4; i++) {
+ for (int i = 3; i > -1; i--) {
delete[] cmul[i];
delete[] cmullr[i];
}