diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index cb662ce6b1a37afb3680d3ba943afa822a6ad07c..99a648cf9aabf74953b687a4256a832679299b80 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -94,6 +94,8 @@ void cluster() {
const int ndi = c4->nsph * c4->nlim;
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
double *gaps = new double[nsph]();
+ double *tqev = new double[3]();
+ double *tqsv = new double[3]();
double **tqse, **tqss, **tqce, **tqcs;
complex<double> **tqspe, **tqsps, **tqcpe, **tqcps;
tqse = new double*[2];
@@ -114,6 +116,7 @@ void cluster() {
tqcs[ti] = new double[3]();
tqcps[ti] = new complex<double>[3]();
}
+ double *gapv = new double[3]();
complex<double> **gapp, **gappm;
double **gap, **gapm;
gapp = new complex<double>*[3];
@@ -140,11 +143,16 @@ void cluster() {
double *args = new double[1]();
double *duk = new double[3]();
double **cextlr, **cext;
+ double **cmullr, **cmul;
cextlr = new double*[4];
cext = new double*[4];
+ cmullr = new double*[4];;
+ cmul = new double*[4];
for (int ci = 0; ci < 4; ci++) {
cextlr[ci] = new double[4]();
cext[ci] = new double[4]();
+ cmullr[ci] = new double[4]();
+ cmul[ci] = new double[4]();
}
int isq, ibf;
double scan, cfmp, sfmp, cfsp, sfsp;
@@ -558,12 +566,311 @@ void cluster() {
fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
fprintf(output, " Fk=%15.7lE\n", fz);
} // ilr210 loop
- // from RMBRIF
+ double rmbrif = (c1ao->fsacm[0][0].real() - c1ao->fsacm[1][1].real()) / c1ao->fsacm[0][0].real();
+ double rmdchr = (c1ao->fsacm[0][0].imag() - c1ao->fsacm[1][1].imag()) / c1ao->fsacm[0][0].imag();
+ fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rmbrif);
+ fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rmdchr);
}
// label 212
- } // jphs loop
- } // jths loop
+ fprintf(output, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n", jth486, jph484, jths, jphs);
+ fprintf(output, " TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n", th, ph, ths, phs);
+ fprintf(output, " SCAND=%10.3lE\n", scan);
+ fprintf(output, " CFMP=%15.7lE, SFMP=%15.7lE\n", cfmp, sfmp);
+ fprintf(output, " CFSP=%15.7lE, SFSP=%15.7lE\n", cfsp, sfsp);
+ if (isam >= 0) {
+ fprintf(output, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
+ fprintf(output, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n", uns[0], uns[1], uns[2]);
+ } else { // label 214
+ fprintf(output, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", un[0], un[1], un[2]);
+ }
+ // label 220
+ if (inpol == 0) {
+ fprintf(output, " LIN\n");
+ } else { // label 222
+ fprintf(output, " CIRC\n");
+ }
+ // label 224
+ scr2(vk, vkarg, exri, duk, c1, c1ao, c3, c4);
+ if (c4->li != c4->le) fprintf(output, " SPHERES; LMX=MIN0(LI,LE)\n");
+ for (int i226 = 1; i226 <= nsph; i226++) {
+ if (c1->iog[i226 - 1] >= i226) {
+ fprintf(output, " SPHERE %2d\n", i226);
+ fprintf(
+ output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
+ c1->sas[i226 - 1][0][0].real(), c1->sas[i226 - 1][0][0].imag(),
+ c1->sas[i226 - 1][1][0].real(), c1->sas[i226 - 1][1][0].imag()
+ );
+ fprintf(
+ output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",
+ c1->sas[i226 - 1][0][1].real(), c1->sas[i226 - 1][0][1].imag(),
+ c1->sas[i226 - 1][1][1].real(), c1->sas[i226 - 1][1][1].imag()
+ );
+ for (int j225 = 0; j225 < 16; j225++) { // QUESTION: check that 16 is a fixed dimension
+ c1ao->vint[j225] = c1ao->vints[i226 - 1][j225];
+ } // j225 loop
+ mmulc(c1ao->vint, cmullr, cmul);
+ fprintf(output, " MULS\n");
+ for (int i1 = 0; i1 < 4; i1++) {
+ fprintf(
+ output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ cmul[i1][0], cmul[i1][1], cmul[i1][2], cmul[i1][3]
+ );
+ } // i1 loop
+ fprintf(output, " MULSLR\n");
+ for (int i1 = 0; i1 < 4; i1++) {
+ fprintf(
+ output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ cmullr[i1][0], cmullr[i1][1], cmullr[i1][2], cmullr[i1][3]
+ );
+ } // i1 loop
+ }
+ } // i226 loop
+ fprintf(
+ output, " SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n",
+ c3->tsas[0][0].real(), c3->tsas[0][0].imag(),
+ c3->tsas[1][0].real(), c3->tsas[1][0].imag()
+ );
+ fprintf(
+ output, " SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n",
+ c3->tsas[0][1].real(), c3->tsas[0][1].imag(),
+ c3->tsas[1][1].real(), c3->tsas[1][1].imag()
+ );
+ fprintf(output, " CLUSTER\n");
+ pcros(vk, exri, c1, c1ao, c4);
+ mextc(vk, exri, c1ao->fsac, cextlr, cext);
+ mmulc(c1ao->vint, cmullr, cmul);
+ if (jw != 0) {
+ jw = 0;
+ // Some implicit loops writing to binary.
+ for (int i = 0; i < 4; i++) {
+ for (int j = 0; j < 4; j++) {
+ double value = cext[i][j];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ }
+ for (int i = 0; i < 2; i++) {
+ double value = c1ao->scsc[i];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->scscp[i].real();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->scscp[i].imag();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->ecsc[i];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->ecscp[i].real();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->ecscp[i].imag();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 2; j++) {
+ double value = gap[i][j];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = gapp[i][j].real();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = gapp[i][j].imag();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ }
+ for (int i = 0; i < 2; i++) {
+ for (int j = 0; j < 3; j++) {
+ double value = tqce[i][j];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = tqcpe[i][j].real();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = tqcpe[i][j].imag();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ }
+ for (int i = 0; i < 2; i++) {
+ for (int j = 0; j < 3; j++) {
+ double value = tqcs[i][j];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = tqcps[i][j].real();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = tqcps[i][j].imag();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ }
+ for (int i = 0; i < 3; i++) {
+ double value = u[i];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = up[i];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = un[i];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ }
+ // label 254
+ for (int i = 0; i < 16; i++) {
+ double value = c1ao->vint[i].real();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->vint[i].imag();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ for (int i = 0; i < 4; i++) {
+ for (int j = 0; j < 4; j++) {
+ double value = cmul[i][j];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ }
+ int jlr = 2;
+ for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
+ int ipol = (ilr290 % 2 == 0) ? 1 : -1;
+ if (ilr290 == 2) jlr = 1;
+ double extsec = c1ao->ecsc[ilr290 - 1];
+ double qext = extsec * sqsfi / c3->gcs;
+ double extrat = extsec / c3->ecs;
+ double scasec = c1ao->scsc[ilr290 - 1];
+ double albedc = scasec / extsec;
+ double qsca = scasec * sqsfi / c3->gcs;
+ double scarat = scasec / c3->scs;
+ double abssec = extsec - scasec;
+ double qabs = abssec * sqsfi / c3->gcs;
+ double absrat = 1.0;
+ double ratio = c3->acs / c3->ecs;
+ if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / c3->acs;
+ s0 = c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
+ double qschu = s0.imag() * csch;
+ double pschu = s0.real() * csch;
+ s0mag = sqrt((s0.real() + s0.imag()) * (s0.real() - s0.imag())) * cs0;
+ double refinr = c1ao->fsac[ilr290 - 1][ilr290 - 1].real() / c3->tfsas.real();
+ double extcor = c1ao->fsac[ilr290 - 1][ilr290 - 1].imag() / c3->tfsas.imag();
+ if (inpol == 0) {
+ fprintf(output, " LIN %2d\n", ipol);
+ } else { // label 273
+ fprintf(output, " CIRC %2d\n", ipol);
+ }
+ // label 275
+ fprintf(output, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");
+ fprintf(
+ output, " %14.7lE%15.7lE%15.7lE%15.7lE\n",
+ scasec, abssec, extsec, albedc
+ );
+ fprintf(output, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");
+ fprintf(
+ output, " %14.7lE%15.7lE%15.7lE\n",
+ qsca, qabs, qext
+ );
+ fprintf(output, " ---- SCCRT --- ABCRT --- EXCRT ----\n");
+ fprintf(
+ output, " %14.7lE%15.7lE%15.7lE\n",
+ scarat, absrat, extrat
+ );
+ fprintf(
+ output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
+ ilr290, ilr290, c1ao->fsac[ilr290 - 1][ilr290 - 1].real(), c1ao->fsac[ilr290 - 1][ilr290 - 1].imag(),
+ jlr, ilr290, c1ao->fsac[jlr - 1][ilr290 - 1].real(), c1ao->fsac[jlr - 1][ilr290 - 1].imag()
+ );
+ fprintf(
+ output, " SAC(%1d,%1d)=%15.7lE%15.7lE SAC(%1d,%1d)=%15.7lE%15.7lE\n",
+ ilr290, ilr290, c1ao->sac[ilr290 - 1][ilr290 - 1].real(), c1ao->sac[ilr290 - 1][ilr290 - 1].imag(),
+ jlr, ilr290, c1ao->sac[jlr - 1][ilr290 - 1].real(), c1ao->sac[jlr - 1][ilr290 - 1].imag()
+ );
+ fprintf(
+ output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
+ ilr290, ilr290, refinr, ilr290, ilr290, extcor
+ );
+ fprintf(
+ output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n",
+ qschu, pschu, s0mag
+ );
+ bool goto190 = isam >= 0 && (jths > 1 || jphs > 1);
+ if (!goto190) {
+ gapv[0] = gap[0][ilr290 - 1];
+ gapv[1] = gap[1][ilr290 - 1];
+ gapv[2] = gap[2][ilr290 - 1];
+ double extins = c1ao->ecsc[ilr290 - 1];
+ double scatts = c1ao->scsc[ilr290 - 1];
+ double rapr, cosav, fp, fn, fk, fx, fy, fz;
+ rftr(u, up, un, gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz);
+ fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
+ fprintf(output, " Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n", fp, fn, fk);
+ fprintf(output, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n", fx, fy, fz);
+ tqev[0] = tqce[ilr290 - 1][0];
+ tqev[1] = tqce[ilr290 - 1][1];
+ tqev[2] = tqce[ilr290 - 1][2];
+ tqsv[0] = tqcs[ilr290 - 1][0];
+ tqsv[1] = tqcs[ilr290 - 1][1];
+ tqsv[2] = tqcs[ilr290 - 1][2];
+ double tep, ten, tek, tsp, tsn, tsk;
+ tqr(u, up, un, tqev, tqsv, tep, ten, tek, tsp, tsn, tsk);
+ fprintf(output, " TQEl=%15.7lE, TQEr=%15.7lE, TQEk=%15.7lE\n", tep, ten, tek);
+ fprintf(output, " TQSl=%15.7lE, TQSr=%15.7lE, TQSk=%15.7lE\n", tsp, tsn, tsk);
+ fprintf(
+ output, " TQEx=%15.7lE, TQEy=%15.7lE, TQEz=%15.7lE\n",
+ tqce[ilr290 - 1][0], tqce[ilr290 - 1][1], tqce[ilr290 - 1][2]
+ );
+ fprintf(
+ output, " TQSx=%15.7lE, TQSy=%15.7lE, TQSz=%15.7lE\n",
+ tqcs[ilr290 - 1][0], tqcs[ilr290 - 1][1], tqcs[ilr290 - 1][2]
+ );
+ }
+ } //ilr290 loop
+ double rbirif = (c1ao->fsac[0][0].real() - c1ao->fsac[1][1].real()) / c1ao->fsac[0][0].real();
+ double rdichr = (c1ao->fsac[0][0].imag() - c1ao->fsac[1][1].imag()) / c1ao->fsac[0][0].imag();
+ fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rbirif);
+ fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rdichr);
+ fprintf(output, " MULC\n");
+ for (int i = 0; i < 4; i++) {
+ fprintf(
+ output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ cmul[i][0], cmul[i][1], cmul[i][2], cmul[i][3]
+ );
+ }
+ fprintf(output, " MULCLR\n");
+ for (int i = 0; i < 4; i++) {
+ fprintf(
+ output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ cmullr[i][0], cmullr[i][1], cmullr[i][2], cmullr[i][3]
+ );
+ }
+ if (iavm != 0) {
+ mmulc(c1ao->vintm, cmullr, cmul);
+ // Some implicit loops writing to binary.
+ for (int i = 0; i < 16; i++) {
+ double value;
+ value = c1ao->vintm[i].real();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->vintm[i].imag();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ for (int i = 0; i < 4; i++) {
+ for (int j = 0; j < 4; j++) {
+ double value = cmul[i][j];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ }
+ fprintf(output, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
+ if (inpol == 0) {
+ fprintf(output, " LIN\n");
+ } else { // label 316
+ fprintf(output, " CIRC\n");
+ }
+ // label 318
+ fprintf(output, " MULC\n");
+ for (int i = 0; i < 4; i++) {
+ fprintf(
+ output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ cmul[i][0], cmul[i][1], cmul[i][2], cmul[i][3]
+ );
+ }
+ fprintf(output, " MULCLR\n");
+ for (int i = 0; i < 4; i++) {
+ fprintf(
+ output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
+ cmullr[i][0], cmullr[i][1], cmullr[i][2], cmullr[i][3]
+ );
+ }
+ }
+ // label 420, continues jphs loop
+ if (isam < 1) phs += phsstp;
+ } // jphs loop, labeled 480
+ if (isam <= 1) thsl += thsstp;
+ } // jths loop, labeled 482
+ ph += phstp;
} // jph484 loop
+ th += thstp;
} // jth486 loop
printf("INFO: done jxi488 iteration.\n");
} // jxi488 loop
@@ -609,6 +916,8 @@ void cluster() {
delete[] tqcpe;
delete[] tqcs;
delete[] tqcps;
+ delete[] tqev;
+ delete[] tqsv;
for (int gi = 2; gi > -1; gi--) {
delete[] gapp[gi];
delete[] gappm[gi];
@@ -619,6 +928,7 @@ void cluster() {
delete[] gappm;
delete[] gap;
delete[] gapm;
+ delete[] gapv;
delete[] u;
delete[] us;
delete[] un;
@@ -635,9 +945,13 @@ void cluster() {
for (int ci = 3; ci > -1; ci--) {
delete[] cextlr[ci];
delete[] cext[ci];
+ delete[] cmullr[ci];
+ delete[] cmul[ci];
}
delete[] cextlr;
delete[] cext;
+ delete[] cmullr;
+ delete[] cmul;
} else { // NSPH mismatch between geometry and scatterer configurations.
throw UnrecognizedConfigurationException(
"Inconsistent geometry and scatterer configurations."
diff --git a/src/include/Commons.h b/src/include/Commons.h
index 812f8a74d34cb916f7f13a8572f80bc4d16b0030..7a45deafdbc46ab08d11b59443d241d4bd19f22b 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -238,8 +238,12 @@ public:
//! \brief QUESTION: definition?
std::complex<double> **fsacm;
//! \brief QUESTION: definition?
+ double *scsc;
+ //! \brief QUESTION: definition?
std::complex<double> *scscp;
//! \brief QUESTION: definition?
+ double *ecsc;
+ //! \brief QUESTION: definition?
double *ecscm;
//! \brief QUESTION: definition?
double *scscm;
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 38411b1a9f7707d9fbb10f6d64cca44624cc58fb..3e6c9a393b909734d094bd8ec5e3d718221606b3 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -34,6 +34,7 @@ extern void rabas(
);
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(
@@ -57,7 +58,6 @@ extern void wmasp(
int nsph, double *argi, double *args, C1 *c1
);
// >>> END OF SPH_SUBS DECLARATION <<<
-void logmat(std::complex<double> **mat, int rows, int cols);
/*! \brief C++ porting of CDTP
*
@@ -1415,6 +1415,74 @@ void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr,
}
}
+/*! \brief C++ porting of PCROS
+ *
+ * This function is intended to evaluate the particle cross-section. QUESTIUON: correct?
+ * \param vk: `double`
+ * \param exri: `double`
+ * \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) * amp);
+ 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`
@@ -1800,6 +1868,39 @@ void raba(
delete[] ctqcs;
}
+/*! \brief C++ porting of RFTR
+ *
+ * \param u: `double *`
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param gapv: `double *`
+ * \param extins: `double`
+ * \param scatts: `double`
+ * \param rapr: `double &`
+ * \param cosav: `double &`
+ * \param fp: `double &`
+ * \param fn: `double &`
+ * \param fk: `double &`
+ * \param fx: `double &`
+ * \param fy: `double &`
+ * \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];
+}
+
/*! \brief C++ porting of SCR0
*
* \param vk: `double` QUESTION: definition?
@@ -1852,6 +1953,99 @@ void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
} // i14 loop
}
+/*! \brief C++ porting of SCR2
+ *
+ * \param vk: `double` QUESTION: definition?
+ * \param vkarg: `double` QUESTION: definition?
+ * \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.
+ * \param c3: `C3 *` Pointer to a C3 instance.
+ * \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
+}
+
/*! \brief C++ porting of STR
*
* \param rcf: `double **` Matrix of double.
@@ -1936,6 +2130,32 @@ void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6) {
delete[] ylm;
}
+/*! \brief C++ porting of TQR
+ *
+ * \param u: `double *`
+ * \param up: `double *`
+ * \param un: `double *`
+ * \param tqev: `double *`
+ * \param tqsv: `double *`
+ * \param tep: `double &`
+ * \param ten: `double &`
+ * \param tek: `double &`
+ * \param tsp: `double &`
+ * \param tsn: `double &`
+ * \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];
+}
+
/*! \brief C++ porting of ZTM
*
* \param am: Matrix of complex.
@@ -2025,23 +2245,6 @@ void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9
} // i0 loop
}
-/*! \brief Write a matrix to a log file (debug function).
- *
- * \param mat: Matrix of complex.
- * \param rows: `int`
- * \param cols: `int`
- */
-void logmat(std::complex<double> **mat, int rows, int cols) {
- FILE* logfile = fopen("c_matrix.log", "w");
- for (int i = 0; i < rows; i++) {
- for (int j = 0; j < cols; j++) {
- fprintf(logfile, "R:%3d,%3d,%15.7lE\n", i + 1, j + 1, mat[i][j].real());
- fprintf(logfile, "I:%3d,%3d,%15.7lE\n", i + 1, j + 1, mat[i][j].imag());
- }
- }
- fclose(logfile);
-}
-
/*! \brief Sum all the elements of a matrix (debug function).
*
* \param mat: Matrix of complex.
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 71b6073654a6f283669ae0f087a6452ab4af5d30..6717bac3a03d13b933ff57bcbe46cfc9ae85fd3a 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -127,6 +127,8 @@ C1_AddOns::C1_AddOns(C4 *c4) {
sscs = new double[nsph]();
ecscm = new double[2]();
scscm = new double[2]();
+ scsc = new double[2]();
+ ecsc = new double[2]();
}
C1_AddOns::~C1_AddOns() {
@@ -161,6 +163,8 @@ C1_AddOns::~C1_AddOns() {
delete[] ecscpm;
delete[] ecscm;
delete[] scscm;
+ delete[] scsc;
+ delete[] ecsc;
}
void C1_AddOns::allocate_vectors(C4 *c4) {