#include <cstdio>
#include <fstream>
#include <string>
#include <complex>
#ifndef INCLUDE_CONFIGURATION_H_
#include "../include/Configuration.h"
#endif
#ifndef INCLUDE_CLU_SUBS_H_
#include "../include/clu_subs.h"
#endif
using namespace std;
/*
* >>> WARNING: works only for IDFC >= 0, as the original code <<<
*
*/
//! \brief C++ implementation of CLU
void cluster() {
printf("INFO: making legacy configuration ...\n");
ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/cluster/DEDFB");
conf->write_formatted("c_OEDFB_clu");
conf->write_binary("c_TEDF_clu");
delete conf;
printf("INFO: reading binary configuration ...\n");
ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF_clu");
GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/cluster/DCLU");
if (sconf->number_of_spheres == gconf->number_of_spheres) {
// Shortcuts to variables stored in configuration objects
int nsph = gconf->number_of_spheres;
int mxndm = gconf->mxndm;
int inpol = gconf->in_pol;
int npnt = gconf->npnt;
int npntts = gconf->npntts;
int iavm = gconf->iavm;
int isam = gconf->meridional_type;
int nxi = sconf->number_of_scales;
double th = gconf->in_theta_start;
double thstp = gconf->in_theta_step;
double thlst = gconf->in_theta_end;
double ths = gconf->sc_theta_start;
double thsstp = gconf->sc_theta_step;
double thslst = gconf->sc_theta_end;
double ph = gconf->in_phi_start;
double phstp = gconf->in_phi_step;
double phlst = gconf->in_phi_end;
double phs = gconf->sc_phi_start;
double phsstp = gconf->sc_phi_step;
double phslst = gconf->sc_phi_end;
double wp = sconf->wp;
// Global variables for CLU
int lm = gconf->l_max;
if (gconf->li > lm) lm = gconf->li;
if (gconf->le > lm) lm = gconf->le;
C1 *c1 = new C1(nsph, lm, sconf->nshl_vec, sconf->iog_vec);
C3 *c3 = new C3();
for (int c1i = 0; c1i < nsph; c1i++) {
c1->rxx[c1i] = gconf->sph_x[c1i];
c1->ryy[c1i] = gconf->sph_y[c1i];
c1->rzz[c1i] = gconf->sph_z[c1i];
c1->ros[c1i] = sconf->radii_of_spheres[c1i];
}
C4 *c4 = new C4;
c4->li = gconf->li;
c4->le = gconf->le;
c4->lm = (c4->li > c4-> le) ? c4->li : c4->le;
c4->nv3j = (c4->lm * (c4->lm + 1) * (2 * c4->lm + 7)) / 6;
c4->nsph = nsph;
// The following is needed to initialize C1_AddOns
c4->litpo = c4->li + c4->li + 1;
c4->litpos = c4->litpo * c4->litpo;
c4->lmtpo = c4->li + c4->le + 1;
c4->lmtpos = c4->lmtpo * c4->lmtpo;
c4->nlim = c4->li * (c4->li + 2);
c4->nlem = c4->le * (c4->le + 2);
c4->lmpo = c4->lm + 1;
C1_AddOns *c1ao = new C1_AddOns(c4);
// End of add-ons initialization
C6 *c6 = new C6(c4->lmtpo);
FILE *output = fopen("c_OCLU", "w");
int jer = 0, lcalc = 0;
complex<double> arg(0.0, 0.0), ccsam(0.0, 0.0);
int max_ici = 0;
for (int insh = 0; insh < nsph; insh++) {
int nsh = sconf->nshl_vec[insh];
int ici = (nsh + 1) / 2;
if (ici > max_ici) max_ici = ici;
}
C2 *c2 = new C2(nsph, max_ici, npnt, npntts);
complex<double> **am = new complex<double>*[mxndm];
for (int ai = 0; ai < mxndm; ai++) am[ai] = new complex<double>[mxndm];
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];
tqspe = new complex<double>*[2];
tqss = new double*[2];
tqsps = new complex<double>*[2];
tqce = new double*[2];
tqcpe = new complex<double>*[2];
tqcs = new double*[2];
tqcps = new complex<double>*[2];
for (int ti = 0; ti < 2; ti++) {
tqse[ti] = new double[nsph]();
tqspe[ti] = new complex<double>[nsph]();
tqss[ti] = new double[nsph]();
tqsps[ti] = new complex<double>[nsph]();
tqce[ti] = new double[3]();
tqcpe[ti] = new complex<double>[3]();
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];
gappm = new complex<double>*[3];
gap = new double*[3];
gapm = new double*[3];
for (int gi = 0; gi < 3; gi++) {
gapp[gi] = new complex<double>[2]();
gappm[gi] = new complex<double>[2]();
gap[gi] = new double[2]();
gapm[gi] = new double[2]();
}
double *u = new double[3]();
double *us = new double[3]();
double *un = new double[3]();
double *uns = new double[3]();
double *up = new double[3]();
double *ups = new double[3]();
double *unmp = new double[3]();
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;
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;
// 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",
nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam
);
fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
for (int ri = 0; ri < nsph; ri++) fprintf(output, "%17.8lE%17.8lE%17.8lE\n",
gconf->sph_x[ri], gconf->sph_y[ri], gconf->sph_z[ri]
);
fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
fprintf(
output, " %10.3lE%10.3lE%10.3lE%10.3lE%10.3lE%10.3lE\n",
th, thstp, thlst, ths, thsstp, thslst
);
fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
fprintf(
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, " %5d\n", gconf->jwtm);
fprintf(output, " READ(ITIN)NSPHT\n");
fprintf(output, " READ(ITIN)(IOG(I),I=1,NSPH)\n");
fprintf(output, " READ(ITIN)EXDC,WP,XIP,IDFC,NXI\n");
fprintf(output, " READ(ITIN)(XIV(I),I=1,NXI)\n");
fprintf(output, " READ(ITIN)NSHL(I),ROS(I)\n");
fprintf(output, " READ(ITIN)(RCF(I,NS),NS=1,NSH)\n");
fprintf(output, " \n");
double small = 1.0e-3;
int nth = 0, nph = 0;
if (thstp != 0.0) nth = int((thlst - th) / thstp + small);
nth++;
if (phstp != 0.0) nph = int((phlst - ph) / phstp + small);
nph++;
int nths = 0, nphs = 0;
double thsca = 0.0;
if (isam > 1) {
nths = 1;
thsca = ths - th;
} else { // ISAM <= 1
if (thsstp == 0.0) nths = 0;
else nths = int ((thslst - ths) / thsstp + small);
nths++;
}
if (isam >= 1) {
nphs = 1;
} else {
if (phsstp == 0.0) nphs = 0;
else nphs = int((phslst - phs) / phsstp + small);
nphs++;
}
int nk = nth * nph;
int nks = nths * nphs;
int nkks = nk * nks;
double th1 = th, ph1 = ph, ths1 = ths, phs1 = phs;
str(sconf->rcf, c1, c1ao, c3, c4, c6);
double ****zpv = new double***[c4->lm]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
for (int zi = 0; zi < c4->lm; zi++) {
zpv[zi] = new double**[3];
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]();
}
}
}
thdps(c4->lm, zpv);
double exri = sqrt(sconf->exdc);
double vk = 0.0; // NOTE: Not really sure it should be initialized at 0
fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
fstream tppoan;
tppoan.open("c_TPPOAN", ios::out | ios::binary);
if (tppoan.is_open()) {
tppoan.write(reinterpret_cast<char *>(&iavm), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&isam), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&inpol), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nxi), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nth), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nph), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nths), sizeof(int));
tppoan.write(reinterpret_cast<char *>(&nphs), sizeof(int));
double wn = wp / 3.0e8;
double sqsfi = 1.0;
if (sconf->idfc < 0) {
vk = sconf->xip * wn;
fprintf(output, " VK=%15.7lE, XI IS SCALE FACTOR FOR LENGTHS\n", vk);
fprintf(output, " \n");
}
for (int jxi488 = 1; jxi488 <= nxi; jxi488++) {
int jaw = 1;
fprintf(output, "========== JXI =%3d ====================\n", jxi488);
double xi = sconf->scale_vec[jxi488 - 1];
double vkarg = 0.0;
if (sconf->idfc >= 0) {
vk = xi * wn;
vkarg = vk;
fprintf(output, " VK=%15.7lE, XI=%15.7lE\n", vk, xi);
} else {
vkarg = xi * vk;
sqsfi = 1.0 / (xi * xi);
fprintf(output, " XI=%15.7lE\n", xi);
}
hjv(exri, vkarg, jer, lcalc, arg, c1, c1ao, c4);
if (jer != 0) {
fprintf(output, " STOP IN HJV\n");
break; // jxi488 loop: goes to memory cleaning and return
}
for (int i132 = 1; i132 <= nsph; i132++) {
int iogi = c1->iog[i132 - 1];
if (iogi != i132) {
for (int l123 = 1; l123 <= gconf->li; 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];
} // l123 loop
} else {
int nsh = sconf->nshl_vec[i132 - 1];
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][jxi488 - 1];
} else {
if (jxi488 == 1) {
for (int ic = 0; ic < ici; ic++)
c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][0];
}
}
if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
dme(
c4->li, i132, npnt, npntts, vkarg, sconf->exdc, exri,
c1, c2, jer, lcalc, arg
);
if (jer != 0) {
fprintf(output, " STOP IN DME\n");
break;
}
}
if (jer != 0) break;
} // i132 loop
cms(am, c1, c1ao, c4, c6);
ccsam = summat(am, 960, 960);
printf("DEBUG: after CMS CCSAM = (%lE,%lE)\n", ccsam.real(), ccsam.imag());
int ndit = 2 * nsph * c4->nlim;
lucin(am, mxndm, ndit, jer);
ccsam = summat(am, 960, 960);
printf("DEBUG: after LUCIN CCSAM = (%lE,%lE)\n", ccsam.real(), ccsam.imag());
if (jer != 0) break; // jxi488 loop: goes to memory clean
ztm(am, c1, c1ao, c4, c6, c9);
if (sconf->idfc >= 0) {
if (jxi488 == gconf->jwtm) {
int is = 1;
fstream ttms_file;
ttms_file.open("c_TTMS", ios::out | ios::binary);
if (ttms_file.is_open()) {
ttms_file.write(reinterpret_cast<char *>(&is), sizeof(int));
ttms_file.write(reinterpret_cast<char *>(&lm), sizeof(int));
ttms_file.write(reinterpret_cast<char *>(&vk), sizeof(double));
ttms_file.write(reinterpret_cast<char *>(&exri), sizeof(double));
int nlemt = 2 * c4->nlem;
for (int ami = 0; ami < nlemt; ami++) {
for (int amj = 0; amj < nlemt; amj++) {
complex<double> value = c1ao->am0m[ami][amj];
double rval = value.real();
double ival = value.imag();
ttms_file.write(reinterpret_cast<char *>(&rval), sizeof(double));
ttms_file.write(reinterpret_cast<char *>(&ival), sizeof(double));
}
}
ttms_file.close();
} else { // Could not open TM file. Should never occur.
printf("ERROR: failed to open TTMS file.\n");
break;
}
}
}
// label 156: continue from here
if (inpol == 0) {
fprintf(output, " LIN\n");
} else { // label 158
fprintf(output, " CIRC\n");
}
// label 160
double cs0 = 0.25 * vk * vk * vk / acos(0.0);
double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
std::complex<double> s0(0.0, 0.0);
scr0(vk, exri, c1, c1ao, c3, c4);
printf("DEBUG: after SCR0 TFSAS = (%lE, %lE)\n", c3->tfsas.real(), c3->tfsas.imag());
double sqk = vk * vk * sconf->exdc;
aps(zpv, c4->li, nsph, c1, sqk, gaps);
rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps);
if (c4->li != c4->le) fprintf(output, " SPHERES; LMX=LI\n");
for (int i170 = 1; i170 <= nsph; i170++) {
if (c1->iog[i170 - 1] >= i170) {
int i = i170 - 1;
double albeds = c1->sscs[i] / c1->sexs[i];
c1->sqscs[i] *= sqsfi;
c1->sqabs[i] *= sqsfi;
c1->sqexs[i] *= sqsfi;
fprintf(output, " SPHERE %2d\n", i170);
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());
}
// label 164
fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE%15.7lE\n", c1->sscs[i], c1->sabs[i], c1->sexs[i], albeds);
fprintf(output, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n");
fprintf(output, " %14.7lE%15.7lE%15.7lE\n", c1->sqscs[i], c1->sqabs[i], c1->sqexs[i]);
fprintf(output, " FSAS=%15.7lE%15.7lE\n", c1->fsas[i].real(), c1->fsas[i].imag());
csch = 2.0 * vk * sqsfi / c1->gcsv[i];
s0 = c1->fsas[i] * exri;
qschu = s0.imag() * csch;
pschu = s0.real() * csch;
s0mag = abs(s0) * cs0;
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
double rapr = c1->sexs[i] - gaps[i];
double cosav = gaps[i] / c1->sscs[i];
fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
fprintf(output, " IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[0][i], tqss[0][i]);
fprintf(output, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[1][i], tqss[1][i]);
}
} // i170 loop
fprintf(output, " FSAT=%15.7lE%15.7lE\n", c3->tfsas.real(), c3->tfsas.imag());
csch = 2.0 * vk * sqsfi / c3->gcs;
s0 = c3->tfsas * exri;
qschu = s0.imag() * csch;
pschu = s0.real() * csch;
s0mag = abs(s0) * cs0;
fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
pcrsm0(vk, exri, inpol, c1, c1ao, c4);
apcra(zpv, c4->le, c1ao->am0m, inpol, sqk, gapm, gappm);
th = th1;
for (int jth486 = 1; jth486 <= nth; jth486++) { // OpenMP portable?
ph = ph1;
double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0;
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) {
wmamp(
0, cost, sint, cosp, sinp, inpol, c4->le, 0,
nsph, argi, u, upmp, unmp, c1
);
// label 182
apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
jw = 1;
}
} else { // label 180, NK == 1 AND JXI488 == 1
if (isam >= 0) {
// label 182
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 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
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
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 = abs(s0) * 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
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
printf("ERROR: failed to open TPPOAN file.\n");
}
fclose(output);
// Clean memory
delete c1;
delete c1ao;
delete c3;
delete c4;
delete c6;
delete c9;
for (int zi = c4->lm - 1; zi > -1; zi--) {
for (int zj = 2; zj > -1; zj--) {
delete[] zpv[zi][zj][1];
delete[] zpv[zi][zj][0];
delete[] zpv[zi][zj];
}
delete[] zpv[zi];
}
delete[] zpv;
for (int ai = mxndm - 1; ai > -1; ai--) delete[] am[ai];
delete[] am;
delete[] gaps;
for (int ti = 1; ti > -1; ti--) {
delete[] tqse[ti];
delete[] tqss[ti];
delete[] tqspe[ti];
delete[] tqsps[ti];
delete[] tqce[ti];
delete[] tqcpe[ti];
delete[] tqcs[ti];
delete[] tqcps[ti];
}
delete[] tqse;
delete[] tqss;
delete[] tqspe;
delete[] tqsps;
delete[] tqce;
delete[] tqcpe;
delete[] tqcs;
delete[] tqcps;
delete[] tqev;
delete[] tqsv;
for (int gi = 2; gi > -1; gi--) {
delete[] gapp[gi];
delete[] gappm[gi];
delete[] gap[gi];
delete[] gapm[gi];
}
delete[] gapp;
delete[] gappm;
delete[] gap;
delete[] gapm;
delete[] gapv;
delete[] u;
delete[] us;
delete[] un;
delete[] uns;
delete[] up;
delete[] ups;
delete[] unmp;
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[] 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."
);
}
delete sconf;
delete gconf;
printf("Done.\n");
}