diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index b685453d6aaf15996a52419508ee5525cfc7a115..4c08ba9e21ff7cdb82e15175d06a9eaf97fbabf6 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -23,950 +23,950 @@ using namespace std;
* \param output_path: `string` Directory to write the output files in.
*/
void cluster(string config_file, string data_file, string output_path) {
- printf("INFO: making legacy configuration...");
- ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(config_file);
- sconf->write_formatted(output_path + "/c_OEDFB_clu");
- sconf->write_binary(output_path + "/c_TEDF_clu");
- GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(data_file);
- printf(" done.\n");
- 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];
+ printf("INFO: making legacy configuration...");
+ ScattererConfiguration *sconf = ScattererConfiguration::from_dedfb(config_file);
+ sconf->write_formatted(output_path + "/c_OEDFB_clu");
+ sconf->write_binary(output_path + "/c_TEDF_clu");
+ GeometryConfiguration *gconf = GeometryConfiguration::from_legacy(data_file);
+ printf(" done.\n");
+ 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((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;
- 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;
+ }
+ ttms_file.close();
+ } else { // Could not open TM file. Should never occur.
+ printf("\nERROR: 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;
}
- 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]();
+ // 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
+ );
}
- 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]();
+ // 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];
+ }
+ }
}
- 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++;
+ // 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;
}
- if (isam >= 1) {
- nphs = 1;
- } else {
- if (phsstp == 0.0) nphs = 0;
- else nphs = int((phslst - phs) / phsstp + small);
- nphs++;
+ // 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);
}
- 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]();
- }
- }
+ // 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]);
}
- 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...\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");
+ // label 220
+ if (inpol == 0) {
+ fprintf(output, " LIN\n");
+ } else { // label 222
+ fprintf(output, " CIRC\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];
+ // 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));
+ }
}
- 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];
+ // 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));
}
- 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];
+ 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));
+ }
}
- 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];
+ 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]
+ );
}
- 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());
+ 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(" 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());
}