From 8d84efcd2c2c94052d6acc5190aac815d1b8eb77 Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Sun, 3 Dec 2023 19:13:02 +0100
Subject: [PATCH] Improve format matching with legacy output
---
src/cluster/cluster.cpp | 184 +++++++++++++++++++++-
src/include/Commons.h | 2 +
src/include/clu_subs.h | 329 +++++++++++++++++++++++++++++++++++++---
src/libnptm/Commons.cpp | 6 +-
4 files changed, 494 insertions(+), 27 deletions(-)
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 7b7e27c5..cb662ce6 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -136,6 +136,18 @@ void cluster() {
double *unsmp = new double[3]();
double *upmp = new double[3]();
double *upsmp = new double[3]();
+ double *argi = new double[1]();
+ double *args = new double[1]();
+ double *duk = new double[3]();
+ double **cextlr, **cext;
+ cextlr = new double*[4];
+ cext = new double*[4];
+ for (int ci = 0; ci < 4; ci++) {
+ cextlr[ci] = new double[4]();
+ cext[ci] = new double[4]();
+ }
+ int isq, ibf;
+ double scan, cfmp, sfmp, cfsp, sfsp;
// End of global variables for CLU
fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");
fprintf(output, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",
@@ -155,7 +167,7 @@ void cluster() {
output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
ph, phstp, phlst, phs, phsstp, phslst
);
- fprintf(output, " READ(IR,*)JWTM\n");
+ fprintf(output, " READ(IR,*)JWTM\n");
fprintf(output, " %5d\n", gconf->jwtm);
fprintf(output, " READ(ITIN)NSPHT\n");
fprintf(output, " READ(ITIN)(IOG(I),I=1,NSPH)\n");
@@ -223,7 +235,7 @@ void cluster() {
if (sconf->idfc < 0) {
vk = sconf->xip * wn;
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
- fprintf(output, " \n");
+ fprintf(output, " \n");
}
for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {
int jaw = 1;
@@ -338,7 +350,7 @@ void cluster() {
if (c1->nshl[i] != 1) {
fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i]);
} else { // label 162
- fprintf(output, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE,%15.7lE\n", c2->vsz[i], c2->vkt[i].real(), c2->vkt[i].imag());
+ fprintf(output, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", c2->vsz[i], c2->vkt[i].real(), c2->vkt[i].imag());
}
// label 164
fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
@@ -373,14 +385,11 @@ void cluster() {
for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP portable?
ph = ph1;
double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
- // argi[NSPEF], with NSPEF=1 if IDOT=0, else NSPEF=NSPH
- double *argi;
for (int jph484 = 1; jph484 <= nph; jph484++) {
int jw = 0;
if (nk != 1 || jxi488 <= 1) {
upvmp(th, ph, 0, cost, sint, cosp, sinp, u, upmp, unmp);
if (isam >= 0) {
- argi = new double[1];
wmamp(
0, cost, sint, cosp, sinp, inpol, c4->le, 0,
nsph, argi, u, upmp, unmp, c1
@@ -400,6 +409,160 @@ void cluster() {
}
// label 184
double thsl = ths1;
+ double phsph = 0.0;
+ for (int jths = 1; jths <= nths; jths++) {
+ ths = thsl;
+ int icspnv = 0;
+ if (isam > 1) ths += thsca;
+ if (isam >= 1) {
+ phsph = 0.0;
+ if (ths < 0.0 || ths > 180.0) phsph = 180.0;
+ if (ths < 0.0) ths *= -1.0;
+ if (ths > 180.0) ths = 360.0 - ths;
+ if (phsph != 0.0) icspnv = 1;
+ }
+ // label 186
+ phs = phs1;
+ for (int jphs = 1; jphs <= nphs; jphs++) {
+ double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0;
+ if (isam >= 1) {
+ phs = ph + phsph;
+ if (phs > 360.0) phs -= 360.0;
+ }
+ // label 188
+ if (!((nks == 1 && jxi488 == 1) || jth486 > 1 || jph484 > 1)) {
+ upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
+ if (isam >= 0)
+ wmamp(
+ 2, costs, sints, cosps, sinps, inpol, c4->le,
+ 0, nsph, args, us, upsmp, unsmp, c1
+ );
+ }
+ // label 190
+ if (nkks != 1 || jxi488 <= 1) {
+ upvsp(
+ u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns,
+ duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp
+ );
+ if (isam < 0) {
+ wmasp(
+ cost, sint, cosp, sinp, costs, sints, cosps, sinps,
+ u, up, un, us, ups, uns, isq, ibf, inpol, c4->le,
+ 0, nsph, argi, args, c1
+ );
+ } else { // label 192
+ for (int i193 = 0; i193 < 3; i193++) {
+ up[i193] = upmp[i193];
+ un[i193] = unmp[i193];
+ ups[i193] = upsmp[i193];
+ uns[i193] = unsmp[i193];
+ }
+ }
+ }
+ // label 194
+ if (iavm == 1) crsm1(vk, exri, c1, c1ao, c4, c6);
+ if (isam < 0) {
+ apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
+ raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
+ jw = 1;
+ }
+ // label 196
+ tppoan.write(reinterpret_cast<char *>(&th), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&ph), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&ths), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&phs), sizeof(double));
+ tppoan.write(reinterpret_cast<char *>(&scan), sizeof(double));
+ if (jaw != 0) {
+ jaw = 0;
+ mextc(vk, exri, c1ao->fsacm, cextlr, cext);
+ // We now have 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 < 3; i++) {
+ double value = c1ao->scscm[i];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->scscpm[i].real();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->scscpm[i].imag();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->ecscm[i];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->ecscpm[i].real();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = c1ao->ecscpm[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 = gapm[i][j];
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = gappm[i][j].real();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ value = gappm[i][j].imag();
+ tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
+ }
+ }
+ fprintf(output, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n", iavm);
+ int jlr = 2;
+ for (int ilr210 = 1; ilr210 <= 2; ilr210++) {
+ int ipol = (ilr210 % 2 == 0) ? 1 : -1;
+ if (ilr210 == 2) jlr = 1;
+ double extsm = c1ao->ecscm[ilr210 - 1];
+ double qextm = extsm * sqsfi / c3->gcs;
+ double extrm = extsm / c3->ecs;
+ double scasm = c1ao->scscm[ilr210 - 1];
+ double albdm = scasm / extsm;
+ double qscam = scasm *sqsfi / c3->gcs;
+ double scarm = scasm / c3->scs;
+ double abssm = extsm - scasm;
+ double qabsm = abssm * sqsfi / c3->gcs;
+ double absrm = abssm / c3->acs;
+ double acsecs = c3->acs / c3->ecs;
+ if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;
+ complex<double> s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
+ double qschum = s0m.imag() * csch;
+ double pschum = s0m.real() * csch;
+ double s0magm = sqrt((s0m.real() + s0m.imag()) * (s0m.real() - s0m.imag())) * cs0;
+ double rfinrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].real() / c3->tfsas.real();
+ double extcrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag() / c3->tfsas.imag();
+ if (inpol == 0) {
+ fprintf(output, " LIN %2d\n", ipol);
+ } else { // label 206
+ fprintf(output, " CIRC %2d\n", ipol);
+ }
+ // label 208
+ fprintf(output, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n");
+ fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", scasm, abssm, extsm, albdm);
+ fprintf(output, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n");
+ fprintf(output, " %14.7lE%15.7lE%15.7lE\n", qscam, qabsm, qextm);
+ fprintf(output, " ---- SCCRT --- ABCRT --- EXCRT ----\n");
+ fprintf(output, " %14.7lE%15.7lE%15.7lE\n", scarm, absrm, extrm);
+ fprintf(
+ output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",
+ ilr210, ilr210, c1ao->fsacm[ilr210 - 1][ilr210 - 1].real(),
+ c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag(), jlr, ilr210,
+ c1ao->fsacm[jlr - 1][ilr210 - 1].real(), c1ao->fsacm[jlr - 1][ilr210 - 1].imag()
+ );
+ fprintf(
+ output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
+ ilr210, ilr210, rfinrm, ilr210, ilr210, extcrm
+ );
+ fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschum, pschum, s0magm);
+ double rapr = c1ao->ecscm[ilr210 - 1] - gapm[2][ilr210 - 1];
+ double cosav = gapm[2][ilr210 - 1] / c1ao->scscm[ilr210 - 1];
+ double fz = rapr;
+ fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
+ fprintf(output, " Fk=%15.7lE\n", fz);
+ } // ilr210 loop
+ // from RMBRIF
+ }
+ // label 212
+ } // jphs loop
+ } // jths loop
} // jph484 loop
} // jth486 loop
printf("INFO: done jxi488 iteration.\n");
@@ -466,6 +629,15 @@ void cluster() {
delete[] unsmp;
delete[] upmp;
delete[] upsmp;
+ delete[] argi;
+ delete[] args;
+ delete[] duk;
+ for (int ci = 3; ci > -1; ci--) {
+ delete[] cextlr[ci];
+ delete[] cext[ci];
+ }
+ delete[] cextlr;
+ delete[] cext;
} 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 81f0cf1b..812f8a74 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -228,6 +228,8 @@ public:
//! \brief QUESTION: definition?
std::complex<double> *vintm;
//! \brief QUESTION: definition?
+ std::complex<double> **vints;
+ //! \brief QUESTION: definition?
std::complex<double> *vintt;
//! \brief QUESTION: definition?
std::complex<double> **fsac;
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index f49166b6..38411b1a 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -40,11 +40,22 @@ extern void upvmp(
double thd, double phd, int icspnv, double &cost, double &sint,
double &cosp, double &sinp, double *u, double *up, double *un
);
+extern void upvsp(
+ double *u, double *upmp, double *unmp, double *us, double *upsmp, double *unsmp,
+ double *up, double *un, double *ups, double *uns, double *duk, int &isq,
+ int &ibf, double &scand, double &cfmp, double &sfmp, double &cfsp, double &sfsp
+);
extern void wmamp(
int iis, double cost, double sint, double cosp, double sinp, int inpol,
int lm, int idot, int nsph, double *arg, double *u, double *up,
double *un, C1 *c1
);
+extern void wmasp(
+ double cost, double sint, double cosp, double sinp, double costs, double sints,
+ double cosps, double sinps, double *u, double *up, double *un, double *us,
+ double *ups, double *uns, int isq, int ibf, int inpol, int lm, int idot,
+ int nsph, double *argi, double *args, C1 *c1
+);
// >>> END OF SPH_SUBS DECLARATION <<<
void logmat(std::complex<double> **mat, int rows, int cols);
@@ -245,6 +256,109 @@ void r3jjr(int j2, int j3, int m2, int m3, C6 *c6) {
}
}
+/*! \brief C++ porting of R3JMR
+ *
+ * \param j1: `int`
+ * \param j2: `int`
+ * \param j3: `int`
+ * \param m1: `int`
+ * \param c6: `C6 *`
+ */
+void r3jmr(int j1, int j2, int j3, int m1, C6 *c6) {
+ int mmx = (j2 < j3 - m1) ? j2 : j3 - m1;
+ int mmn = (-j2 > -(j3 + m1)) ? -j2 : -(j3 + m1);
+ int nmmo = mmx - mmn;
+ int j1po = j1 + 1;
+ int j1tpo = j1po + j1;
+ int isn = 1;
+ if ((j2 - j3 - m1) % 2 != 0) isn = -1;
+ if (nmmo <= 0) {
+ double sj = 1.0 * j1tpo;
+ double cnr = (1.0 / sqrt(sj)) * isn;
+ c6->rac3j[0] = cnr;
+ // returns
+ } else { // label 15
+ int j1s = j1 * j1po;
+ int j2po = j2 + 1;
+ int j2s = j2 * j2po;
+ int j3po = j3 + 1;
+ int j3s = j3 * j3po;
+ int id = j1s - j2s - j3s;
+ int m2 = mmx;
+ int m3 = m1 + m2;
+ double cm = sqrt(1.0 * (j2po - m2) * (j2 + m2) * (j3po - m3) * (j3 + m3));
+ double dm = 1.0 * (id + m2 * m3 * 2);
+ if (nmmo <= 1) {
+ c6->rac3j[0] = dm / cm;
+ double sj = (1.0 + c6->rac3j[0] * c6->rac3j[0]) * j1tpo;
+ double cnr = 1.0 / sqrt(sj) * isn;
+ c6->rac3j[1] = cnr;
+ c6->rac3j[0] *= cnr;
+ // returns
+ } else { // label 20
+ int nm = nmmo + 1;
+ int nmat = (nm + 1) / 2;
+ c6->rac3j[nm - 1] = 1.0;
+ c6->rac3j[nmmo - 1] = dm / cm;
+ double sjt = 1.0;
+ double sjr = 1.0;
+ if (nmat != nmmo) {
+ int nbr = nmmo - nmat;
+ for (int ibr45 = 1; ibr45 <= nbr; ibr45++) {
+ int irr = nm - ibr45;
+ m2--;
+ m3 = m1 + m2;
+ double cmp = cm;
+ cm = sqrt(1.0 * (j2po - m2) * (j2 + m2) * (j3po - m3) * (j3 + m3));
+ sjt = c6->rac3j[irr - 1] * c6->rac3j[irr - 1];
+ dm = 1.0 * (id + m2 * m3 * 2);
+ c6->rac3j[irr - 1] *= ((dm - c6->rac3j[irr] * cmp) / cm);
+ sjr += sjt;
+ } // ibr45 loop
+ }
+ // label 50
+ double osjt = sjt;
+ sjt = c6->rac3j[nmat - 1] * c6->rac3j[nmat - 1];
+ if (sjt >= osjt) {
+ sjr += sjt;
+ } else { // label 55
+ nmat++;
+ }
+ // label 60
+ double racmat = c6->rac3j[nmat - 1];
+ c6->rac3j[0] = 1.0;
+ m2 = mmn;
+ m3 = m1 + m2;
+ double cmp = sqrt(1.0 * (j2 - m2) * (j2po + m2) * (j3 - m3) * (j3po + m3));
+ dm = 1.0 * (id + m2 * m3 * 2);
+ c6->rac3j[1] = dm / cmp;
+ double sjl = 1.0;
+ int nmatmo = nmat - 1;
+ if (nmatmo > 1) {
+ for (int irl70 = 2; irl70 <= nmatmo; irl70++) {
+ m2++;
+ m3 = m1 + m2;
+ cm = cmp;
+ cmp = sqrt(1.0 * (j2 - m2) * (j2po + m2) * (j3 - m3) * (j3po + m3));
+ sjt = c6->rac3j[irl70 - 1] * c6->rac3j[irl70 - 1];
+ dm = 1.0 * (id + m2 * m3 * 2);
+ c6->rac3j[irl70] = (c6->rac3j[irl70 - 1] * dm - c6->rac3j[irl70 - 2] * cm) / cmp;
+ sjl += sjt;
+ }
+ }// label 75
+ double ratrac = racmat / c6->rac3j[nmat - 1];
+ double rats = ratrac * ratrac;
+ double sj = (sjr + sjl * rats) * j1tpo;
+ c6->rac3j[nmat - 1] = racmat;
+ double cnr = 1.0 / sqrt(sj) * isn;
+ for (int irr80 = nmat; irr80 <= nm; irr80++) c6->rac3j[irr80 - 1] *= cnr;
+ double cnl = cnr * ratrac;
+ for (int irl85 = 1; irl85 <= nmatmo; irl85++) c6->rac3j[irl85 - 1] *= cnl;
+ // returns
+ }
+ }
+}
+
/*! \brief C++ porting of GHIT
*
* \param ihi: `int`
@@ -663,7 +777,7 @@ void apcra(
// label 26
} // ilmp loop
} // l28 loop
- for (int l30 = 1; l30 <= le; l30++) { // This can be omitted, since svs was initialized at 0
+ 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++) {
@@ -757,9 +871,9 @@ void apcra(
sum1 = cc0;
sum2 = cc0;
for (int l68 = 1; l68 <= le; l68++) {
- int lpo = l68 + 1;
- int ltpo = l68 + lpo;
- int imm = l68 * lpo;
+ //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) {
@@ -921,6 +1035,135 @@ void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6) {
} // n1 loop
}
+/*! \brief C++ porting of CRSM1
+ *
+ * \param vk: `double`
+ * \param exri: `double`
+ * \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;
+
+ // 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`
@@ -1114,6 +1357,64 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier) {
delete[] v;
}
+/*! \brief C++ porting of MEXTC
+ *
+ * \param vk: `double`
+ * \param exri: `double`
+ * \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;
+ }
+ }
+}
+
/*! \brief C++ porting of PCRSM0
*
* \param vk: `double`
@@ -1515,7 +1816,7 @@ void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
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;
+ //double scs = 0.0, ecs = 0.0, acs = 0.0;
c3->tfsas = cc0;
for (int i14 = 1; i14 <= c4->nsph; i14++) {
int iogi = c1->iog[i14 - 1];
@@ -1544,9 +1845,9 @@ void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4) {
c1->fsas[i14 - 1] = sum21 * csam;
}
// label 12
- scs += c1->sscs[iogi - 1];
- ecs += c1->sexs[iogi - 1];
- acs += c1->sabs[iogi - 1];
+ 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
}
@@ -1670,18 +1971,6 @@ void ztm(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6, C9
} // im2 loop
} // l2 loop
} // n2 loop
- /* // DEBUG CODE
- std::complex<double> dbgtst(0.0, 0.0);
- for (int di = 0; di < ndi; di++) {
- for (int dj = 0; dj < c4->nlem; dj++) dbgtst += c9->gis[di][dj];
- }
- printf("DEBUG: in ZTM init sgis = (%lE, %lE)\n", dbgtst.real(), dbgtst.imag());
- dbgtst = std::complex<double>(0.0, 0.0);
- for (int di = 0; di < ndi; di++) {
- for (int dj = 0; dj < c4->nlem; dj++) dbgtst += c9->gls[di][dj];
- }
- printf("DEBUG: in ZTM init sgls = (%lE, %lE)\n", dbgtst.real(), dbgtst.imag());
- // END DEBUG */
for (int i1 = 1; i1 <= ndi; i1++) { // GPU portable?
int i1e = i1 + ndi;
for (int i3 = 1; i3 <= c4->nlem; i3++) {
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 62191875..71b60736 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -106,9 +106,11 @@ C1_AddOns::C1_AddOns(C4 *c4) {
for (int ai = 0; ai < nlemt; ai++) {
am0m[ai] = new complex<double>[nlemt]();
}
- vint = new complex<double>[16]();
+ vint = new complex<double>[16](); // QUESTION: is dimension 16 generally fixed?
vintm = new complex<double>[16]();
vintt = new complex<double>[16]();
+ vints = new complex<double>*[nsph];
+ for (int vi = 0; vi < nsph; vi++) vints[vi] = new complex<double>[16]();
fsac = new complex<double>*[2];
sac = new complex<double>*[2];
fsacm = new complex<double>*[2];
@@ -143,6 +145,8 @@ C1_AddOns::~C1_AddOns() {
delete[] vint;
delete[] vintm;
delete[] vintt;
+ for (int vi = nsph - 1; vi > -1; vi--) delete[] vints[vi];
+ delete[] vints;
for (int fi = 1; fi > -1; fi--) {
delete[] fsac[fi];
delete[] sac[fi];
--
GitLab