Skip to content
Snippets Groups Projects
Commit 36c854d4 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Reduce indentation and verbosity in cluster.cpp

parent c33f52e8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/cluster/cluster.cpp
+ 934
− 934
View file @ 36c854d4
...@@ -23,950 +23,950 @@ using namespace std; ...@@ -23,950 +23,950 @@ using namespace std;
* \param output_path: `string` Directory to write the output files in. * \param output_path: `string` Directory to write the output files in.
*/ */
void cluster(string config_file, string data_file, string output_path) { void cluster(string config_file, string data_file, string output_path) {
printf("INFO: making legacy configuration..."); printf("INFO: making legacy configuration...");
ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(config_file); ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(config_file);
sconf->write_formatted(output_path + "/c_OEDFB_clu"); sconf->write_formatted(output_path + "/c_OEDFB_clu");
sconf->write_binary(output_path + "/c_TEDF_clu"); sconf->write_binary(output_path + "/c_TEDF_clu");
GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(data_file); GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(data_file);
printf(" done.\n"); printf(" done.\n");
if (sconf->number_of_spheres == gconf->number_of_spheres) { if (sconf->number_of_spheres == gconf->number_of_spheres) {
// Shortcuts to variables stored in configuration objects // Shortcuts to variables stored in configuration objects
int nsph = gconf->number_of_spheres; int nsph = gconf->number_of_spheres;
int mxndm = gconf->mxndm; int mxndm = gconf->mxndm;
int inpol = gconf->in_pol; int inpol = gconf->in_pol;
int npnt = gconf->npnt; int npnt = gconf->npnt;
int npntts = gconf->npntts; int npntts = gconf->npntts;
int iavm = gconf->iavm; int iavm = gconf->iavm;
int isam = gconf->meridional_type; int isam = gconf->meridional_type;
int nxi = sconf->number_of_scales; int nxi = sconf->number_of_scales;
double th = gconf->in_theta_start; double th = gconf->in_theta_start;
double thstp = gconf->in_theta_step; double thstp = gconf->in_theta_step;
double thlst = gconf->in_theta_end; double thlst = gconf->in_theta_end;
double ths = gconf->sc_theta_start; double ths = gconf->sc_theta_start;
double thsstp = gconf->sc_theta_step; double thsstp = gconf->sc_theta_step;
double thslst = gconf->sc_theta_end; double thslst = gconf->sc_theta_end;
double ph = gconf->in_phi_start; double ph = gconf->in_phi_start;
double phstp = gconf->in_phi_step; double phstp = gconf->in_phi_step;
double phlst = gconf->in_phi_end; double phlst = gconf->in_phi_end;
double phs = gconf->sc_phi_start; double phs = gconf->sc_phi_start;
double phsstp = gconf->sc_phi_step; double phsstp = gconf->sc_phi_step;
double phslst = gconf->sc_phi_end; double phslst = gconf->sc_phi_end;
double wp = sconf->wp; double wp = sconf->wp;
// Global variables for CLU // Global variables for CLU
int lm = gconf->l_max; int lm = gconf->l_max;
if (gconf->li > lm) lm = gconf->li; if (gconf->li > lm) lm = gconf->li;
if (gconf->le > lm) lm = gconf->le; if (gconf->le > lm) lm = gconf->le;
C1 *c1 = new C1(nsph, lm, sconf->nshl_vec, sconf->iog_vec); C1 *c1 = new C1(nsph, lm, sconf->nshl_vec, sconf->iog_vec);
C3 *c3 = new C3(); C3 *c3 = new C3();
for (int c1i = 0; c1i < nsph; c1i++) { for (int c1i = 0; c1i < nsph; c1i++) {
c1->rxx[c1i] = gconf->sph_x[c1i]; c1->rxx[c1i] = gconf->sph_x[c1i];
c1->ryy[c1i] = gconf->sph_y[c1i]; c1->ryy[c1i] = gconf->sph_y[c1i];
c1->rzz[c1i] = gconf->sph_z[c1i]; c1->rzz[c1i] = gconf->sph_z[c1i];
c1->ros[c1i] = sconf->radii_of_spheres[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((output_path + "/c_OCLU").c_str(), "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;
string tppoan_name = output_path + "/c_TPPOAN";
tppoan.open(tppoan_name.c_str(), 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++) {
printf("INFO: running scale iteration %d...", 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
//printf("INFO: initializing matrix...");
cms(am, c1, c1ao, c4, c6);
//printf(" done.\n");
//ccsam = summat(am, 960, 960);
//printf("DEBUG: after CMS CCSAM = (%lE,%lE)\n", ccsam.real(), ccsam.imag());
int ndit = 2 * nsph * c4->nlim;
//printf("INFO: inverting matrix...");
lucin(am, mxndm, ndit, jer);
//printf(" done.\n");
//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;
string ttms_name = output_path + "/c_TTMS";
ttms_file.open(ttms_name.c_str(), 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));
} }
C4 *c4 = new C4; }
c4->li = gconf->li; ttms_file.close();
c4->le = gconf->le; } else { // Could not open TM file. Should never occur.
c4->lm = (c4->li > c4-> le) ? c4->li : c4->le; printf("\nERROR: failed to open TTMS file.\n");
c4->nv3j = (c4->lm * (c4->lm + 1) * (2 * c4->lm + 7)) / 6; break;
c4->nsph = nsph; }
// The following is needed to initialize C1_AddOns }
c4->litpo = c4->li + c4->li + 1; }
c4->litpos = c4->litpo * c4->litpo; // label 156: continue from here
c4->lmtpo = c4->li + c4->le + 1; if (inpol == 0) {
c4->lmtpos = c4->lmtpo * c4->lmtpo; fprintf(output, " LIN\n");
c4->nlim = c4->li * (c4->li + 2); } else { // label 158
c4->nlem = c4->le * (c4->le + 2); fprintf(output, " CIRC\n");
c4->lmpo = c4->lm + 1; }
C1_AddOns *c1ao = new C1_AddOns(c4); // label 160
// End of add-ons initialization double cs0 = 0.25 * vk * vk * vk / acos(0.0);
C6 *c6 = new C6(c4->lmtpo); double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w"); std::complex<double> s0(0.0, 0.0);
int jer = 0, lcalc = 0; scr0(vk, exri, c1, c1ao, c3, c4);
complex<double> arg(0.0, 0.0), ccsam(0.0, 0.0); //printf("DEBUG: after SCR0 TFSAS = (%lE, %lE)\n", c3->tfsas.real(), c3->tfsas.imag());
int max_ici = 0; double sqk = vk * vk * sconf->exdc;
for (int insh = 0; insh < nsph; insh++) { aps(zpv, c4->li, nsph, c1, sqk, gaps);
int nsh = sconf->nshl_vec[insh]; rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps);
int ici = (nsh + 1) / 2; if (c4->li != c4->le) fprintf(output, " SPHERES; LMX=LI\n");
if (ici > max_ici) max_ici = ici; 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;
} }
C2 *c2 = new C2(nsph, max_ici, npnt, npntts); // label 188
complex<double> **am = new complex<double>*[mxndm]; bool goto190 = (nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
for (int ai = 0; ai < mxndm; ai++) am[ai] = new complex<double>[mxndm]; if (!goto190) {
const int ndi = c4->nsph * c4->nlim; upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, us, upsmp, unsmp);
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem); if (isam >= 0)
double *gaps = new double[nsph](); wmamp(
double *tqev = new double[3](); 2, costs, sints, cosps, sinps, inpol, c4->le,
double *tqsv = new double[3](); 0, nsph, args, us, upsmp, unsmp, c1
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](); // label 190
complex<double> **gapp, **gappm; if (nkks != 1 || jxi488 <= 1) {
double **gap, **gapm; upvsp(
gapp = new complex<double>*[3]; u, upmp, unmp, us, upsmp, unsmp, up, un, ups, uns,
gappm = new complex<double>*[3]; duk, isq, ibf, scan, cfmp, sfmp, cfsp, sfsp
gap = new double*[3]; );
gapm = new double*[3]; if (isam < 0) {
for (int gi = 0; gi < 3; gi++) { wmasp(
gapp[gi] = new complex<double>[2](); cost, sint, cosp, sinp, costs, sints, cosps, sinps,
gappm[gi] = new complex<double>[2](); u, up, un, us, ups, uns, isq, ibf, inpol, c4->le,
gap[gi] = new double[2](); 0, nsph, argi, args, c1
gapm[gi] = new double[2](); );
} 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];
}
}
} }
double *u = new double[3](); // label 194
double *us = new double[3](); if (iavm == 1) crsm1(vk, exri, c1, c1ao, c4, c6);
double *un = new double[3](); if (isam < 0) {
double *uns = new double[3](); apc(zpv, c4->le, c1ao->am0m, c1->w, sqk, gap, gapp);
double *up = new double[3](); raba(c4->le, c1ao->am0m, c1->w, tqce, tqcpe, tqcs, tqcps);
double *ups = new double[3](); jw = 1;
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) { // label 196
nphs = 1; tppoan.write(reinterpret_cast<char *>(&th), sizeof(double));
} else { tppoan.write(reinterpret_cast<char *>(&ph), sizeof(double));
if (phsstp == 0.0) nphs = 0; tppoan.write(reinterpret_cast<char *>(&ths), sizeof(double));
else nphs = int((phslst - phs) / phsstp + small); tppoan.write(reinterpret_cast<char *>(&phs), sizeof(double));
nphs++; 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 = abs(s0m) * 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);
} }
int nk = nth * nph; // label 212
int nks = nths * nphs; fprintf(output, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n", jth486, jph484, jths, jphs);
int nkks = nk * nks; fprintf(output, " TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n", th, ph, ths, phs);
double th1 = th, ph1 = ph, ths1 = ths, phs1 = phs; fprintf(output, " SCAND=%10.3lE\n", scan);
str(sconf->rcf, c1, c1ao, c3, c4, c6); fprintf(output, " CFMP=%15.7lE, SFMP=%15.7lE\n", cfmp, sfmp);
double ****zpv = new double***[c4->lm]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2] fprintf(output, " CFSP=%15.7lE, SFSP=%15.7lE\n", cfsp, sfsp);
for (int zi = 0; zi < c4->lm; zi++) { if (isam >= 0) {
zpv[zi] = new double**[3]; fprintf(output, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n", un[0], un[1], un[2]);
for (int zj = 0; zj < 3; zj++) { fprintf(output, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n", uns[0], uns[1], uns[2]);
zpv[zi][zj] = new double*[2]; } else { // label 214
for (int zk = 0; zk < 2; zk++) { fprintf(output, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n", un[0], un[1], un[2]);
zpv[zi][zj][zk] = new double[2]();
}
}
} }
thdps(c4->lm, zpv); // label 220
double exri = sqrt(sconf->exdc); if (inpol == 0) {
double vk = 0.0; // NOTE: Not really sure it should be initialized at 0 fprintf(output, " LIN\n");
fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri); } else { // label 222
fstream tppoan; fprintf(output, " CIRC\n");
string tppoan_name = output_path + "/c_TPPOAN";
tppoan.open(tppoan_name.c_str(), 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++) {
printf("INFO: running scale iteration...\n");
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
printf("INFO: initializing matrix...");
cms(am, c1, c1ao, c4, c6);
printf(" done.\n");
//ccsam = summat(am, 960, 960);
//printf("DEBUG: after CMS CCSAM = (%lE,%lE)\n", ccsam.real(), ccsam.imag());
int ndit = 2 * nsph * c4->nlim;
printf("INFO: inverting matrix...");
lucin(am, mxndm, ndit, jer);
printf(" done.\n");
//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;
string ttms_name = output_path + "/c_TTMS";
ttms_file.open(ttms_name.c_str(), 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
bool goto190 = (nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1));
if (!goto190) {
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 = abs(s0m) * 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 scale.\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); // label 224
// Clean memory scr2(vk, vkarg, exri, duk, c1, c1ao, c3, c4);
delete c1; if (c4->li != c4->le) fprintf(output, " SPHERES; LMX=MIN0(LI,LE)\n");
delete c1ao; for (int i226 = 1; i226 <= nsph; i226++) {
delete c3; if (c1->iog[i226 - 1] >= i226) {
delete c4; fprintf(output, " SPHERE %2d\n", i226);
delete c6; fprintf(
delete c9; output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",
for (int zi = c4->lm - 1; zi > -1; zi--) { c1->sas[i226 - 1][0][0].real(), c1->sas[i226 - 1][0][0].imag(),
for (int zj = 2; zj > -1; zj--) { c1->sas[i226 - 1][1][0].real(), c1->sas[i226 - 1][1][0].imag()
delete[] zpv[zi][zj][1]; );
delete[] zpv[zi][zj][0]; fprintf(
delete[] zpv[zi][zj]; 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(),
delete[] zpv[zi]; 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));
}
} }
delete[] zpv; // label 254
for (int ai = mxndm - 1; ai > -1; ai--) delete[] am[ai]; for (int i = 0; i < 16; i++) {
delete[] am; double value = c1ao->vint[i].real();
delete[] gaps; tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
for (int ti = 1; ti > -1; ti--) { value = c1ao->vint[i].imag();
delete[] tqse[ti]; tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
delete[] tqss[ti];
delete[] tqspe[ti];
delete[] tqsps[ti];
delete[] tqce[ti];
delete[] tqcpe[ti];
delete[] tqcs[ti];
delete[] tqcps[ti];
} }
delete[] tqse; for (int i = 0; i < 4; i++) {
delete[] tqss; for (int j = 0; j < 4; j++) {
delete[] tqspe; double value = cmul[i][j];
delete[] tqsps; tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
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; int jlr = 2;
delete[] gappm; for (int ilr290 = 1; ilr290 <= 2; ilr290++) {
delete[] gap; int ipol = (ilr290 % 2 == 0) ? 1 : -1;
delete[] gapm; if (ilr290 == 2) jlr = 1;
delete[] gapv; double extsec = c1ao->ecsc[ilr290 - 1];
delete[] u; double qext = extsec * sqsfi / c3->gcs;
delete[] us; double extrat = extsec / c3->ecs;
delete[] un; double scasec = c1ao->scsc[ilr290 - 1];
delete[] uns; double albedc = scasec / extsec;
delete[] up; double qsca = scasec * sqsfi / c3->gcs;
delete[] ups; double scarat = scasec / c3->scs;
delete[] unmp; double abssec = extsec - scasec;
delete[] unsmp; double qabs = abssec * sqsfi / c3->gcs;
delete[] upmp; double absrat = 1.0;
delete[] upsmp; double ratio = c3->acs / c3->ecs;
delete[] argi; if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / c3->acs;
delete[] args; s0 = c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
delete[] duk; double qschu = s0.imag() * csch;
for (int ci = 3; ci > -1; ci--) { double pschu = s0.real() * csch;
delete[] cextlr[ci]; s0mag = abs(s0) * cs0;
delete[] cext[ci]; double refinr = c1ao->fsac[ilr290 - 1][ilr290 - 1].real() / c3->tfsas.real();
delete[] cmullr[ci]; double extcor = c1ao->fsac[ilr290 - 1][ilr290 - 1].imag() / c3->tfsas.imag();
delete[] cmul[ci]; 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]
);
} }
delete[] cextlr; fprintf(output, " MULCLR\n");
delete[] cext; for (int i = 0; i < 4; i++) {
delete[] cmullr; fprintf(
delete[] cmul; output, " %15.7lE%15.7lE%15.7lE%15.7lE\n",
} else { // NSPH mismatch between geometry and scatterer configurations. cmullr[i][0], cmullr[i][1], cmullr[i][2], cmullr[i][3]
throw UnrecognizedConfigurationException( );
"Inconsistent geometry and scatterer configurations." }
); if (iavm != 0) {
} mmulc(c1ao->vintm, cmullr, cmul);
delete sconf; // Some implicit loops writing to binary.
delete gconf; for (int i = 0; i < 16; i++) {
printf("Finished: output written to %s.\n", (output_path + "/c_OCLU").c_str()); 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(" done.\n");
} // jxi488 loop
tppoan.close();
} else { // In case TPPOAN could not be opened. Should never happen.
printf("\nERROR: 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("Finished: output written to %s.\n", (output_path + "/c_OCLU").c_str());
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment