diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp index 67dda9a650ac30028f7dda443fd00fb0a1952665..a6d640b8c5d8c16290421e4d68bd46304659f84a 100644 --- a/src/cluster/cluster.cpp +++ b/src/cluster/cluster.cpp @@ -482,7 +482,7 @@ void cluster(const string& config_file, const string& data_file, const string& o #endif #pragma omp barrier - } // close strided loop running on MPI processes + } // close strided loop running on MPI processes, ixi488 loop // delete the shared arrays I used to make available to thread 0 the virtual files of other threads #pragma omp barrier if (myompthread == 0) { @@ -633,769 +633,818 @@ void cluster(const string& config_file, const string& data_file, const string& o fclose(timing_file); delete time_logger; delete logger; +#ifdef MPI_VERSION } +#endif } int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConfiguration *gconf, ScatteringAngles *sa, ClusterIterationData *cid, VirtualAsciiFile *output, const string& output_path, VirtualBinaryFile *vtppoanp) { - int nxi = sconf->number_of_scales; - char virtual_line[256]; - string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n"; - Logger *logger = new Logger(LOG_DEBG); - logger->log(message); - chrono::duration<double> elapsed; - chrono::time_point<chrono::high_resolution_clock> interval_start, interval_end; - int jer = 0; - int lcalc = 0; - int jaw = 1; - int li = gconf->li; - int le = gconf->le; - int lm = 0; - if (le > lm) lm = le; - if (li > lm) lm = li; - int nsph = gconf->number_of_spheres; - np_int mxndm = gconf->mxndm; - int iavm = gconf->iavm; - int inpol = gconf->in_pol; - int npnt = gconf->npnt; - int npntts = gconf->npntts; - int isam = gconf->iavm; - int jwtm = gconf->jwtm; - np_int ndit = 2 * nsph * cid->c4->nlim; - int isq, ibf; + int nxi = sconf->number_of_scales; + char virtual_line[256]; + string message = "INFO: running scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n"; + Logger *logger = new Logger(LOG_DEBG); + logger->log(message); + chrono::duration<double> elapsed; + chrono::time_point<chrono::high_resolution_clock> interval_start, interval_end; + int jer = 0; + int lcalc = 0; + int jaw = 1; + int li = gconf->li; + int le = gconf->le; + int lm = 0; + if (le > lm) lm = le; + if (li > lm) lm = li; + int nsph = gconf->number_of_spheres; + np_int mxndm = gconf->mxndm; + int iavm = gconf->iavm; + int inpol = gconf->in_pol; + int npnt = gconf->npnt; + int npntts = gconf->npntts; + int isam = gconf->iavm; + int jwtm = gconf->jwtm; + np_int ndit = 2 * nsph * cid->c4->nlim; + int isq, ibf; #ifdef USE_NVTX - nvtxRangePush("Prepare matrix calculation"); + nvtxRangePush("Prepare matrix calculation"); #endif - sprintf(virtual_line, "========== JXI =%3d ====================\n\0", jxi488); + sprintf(virtual_line, "========== JXI =%3d ====================\n\0", jxi488); + output->append_line(virtual_line); + double xi = sconf->get_scale(jxi488 - 1); + double exdc = sconf->exdc; + double exri = sqrt(exdc); + int idfc = (int)sconf->idfc; + double vkarg = 0.0; + if (idfc >= 0) { + cid->vk = xi * cid->wn; + vkarg = cid->vk; + sprintf(virtual_line, " VK=%15.7lE, XI=%15.7lE\n\0", cid->vk, xi); output->append_line(virtual_line); - double xi = sconf->get_scale(jxi488 - 1); - double exdc = sconf->exdc; - double exri = sqrt(exdc); - int idfc = (int)sconf->idfc; - double vkarg = 0.0; - if (idfc >= 0) { - cid->vk = xi * cid->wn; - vkarg = cid->vk; - sprintf(virtual_line, " VK=%15.7lE, XI=%15.7lE\n\0", cid->vk, xi); - output->append_line(virtual_line); + } else { + vkarg = xi * cid->vk; + cid->sqsfi = 1.0 / (xi * xi); + sprintf(virtual_line, " XI=%15.7lE\n\0", xi); + output->append_line(virtual_line); + } + hjv(exri, vkarg, jer, lcalc, cid->arg, cid->c1, cid->c1ao, cid->c4); + if (jer != 0) { + sprintf(virtual_line, " STOP IN HJV\n\0"); + output->append_line(virtual_line); + return jer; + // break; // rewrite this to go to the end of the function, to free locally allocated variables and return jer + } + for (int i132 = 1; i132 <= nsph; i132++) { + int iogi = cid->c1->iog[i132 - 1]; + if (iogi != i132) { + for (int l123 = 1; l123 <= li; l123++) { + cid->c1->rmi[l123 - 1][i132 - 1] = cid->c1->rmi[l123 - 1][iogi - 1]; + cid->c1->rei[l123 - 1][i132 - 1] = cid->c1->rei[l123 - 1][iogi - 1]; + } // l123 loop } else { - vkarg = xi * cid->vk; - cid->sqsfi = 1.0 / (xi * xi); - sprintf(virtual_line, " XI=%15.7lE\n\0", xi); - output->append_line(virtual_line); - } - hjv(exri, vkarg, jer, lcalc, cid->arg, cid->c1, cid->c1ao, cid->c4); - if (jer != 0) { - sprintf(virtual_line, " STOP IN HJV\n\0"); - output->append_line(virtual_line); - return jer; - // break; // rewrite this to go to the end of the function, to free locally allocated variables and return jer - } - for (int i132 = 1; i132 <= nsph; i132++) { - int iogi = cid->c1->iog[i132 - 1]; - if (iogi != i132) { - for (int l123 = 1; l123 <= li; l123++) { - cid->c1->rmi[l123 - 1][i132 - 1] = cid->c1->rmi[l123 - 1][iogi - 1]; - cid->c1->rei[l123 - 1][i132 - 1] = cid->c1->rei[l123 - 1][iogi - 1]; - } // l123 loop + int nsh = cid->c1->nshl[i132 - 1]; + int ici = (nsh + 1) / 2; + if (idfc == 0) { + for (int ic = 0; ic < ici; ic++) + cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1); } else { - int nsh = cid->c1->nshl[i132 - 1]; - int ici = (nsh + 1) / 2; - if (idfc == 0) { + if (jxi488 == 1) { for (int ic = 0; ic < ici; ic++) - cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, jxi488 - 1); - } else { - if (jxi488 == 1) { - for (int ic = 0; ic < ici; ic++) - cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0); - } - } - if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc; - dme( - cid->c4->li, i132, npnt, npntts, vkarg, exdc, exri, - cid->c1, cid->c2, jer, lcalc, cid->arg - ); - if (jer != 0) { - sprintf(virtual_line, " STOP IN DME\n\0"); - output->append_line(virtual_line); - return jer; - //break; + cid->c2->dc0[ic] = sconf->get_dielectric_constant(ic, i132 - 1, 0); } } + if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc; + dme( + cid->c4->li, i132, npnt, npntts, vkarg, exdc, exri, + cid->c1, cid->c2, jer, lcalc, cid->arg + ); if (jer != 0) { + sprintf(virtual_line, " STOP IN DME\n\0"); + output->append_line(virtual_line); return jer; //break; } - } // i132 loop + } + if (jer != 0) { + return jer; + //break; + } + } // i132 loop #ifdef USE_NVTX - nvtxRangePop(); + nvtxRangePop(); #endif - interval_start = chrono::high_resolution_clock::now(); + interval_start = chrono::high_resolution_clock::now(); #ifdef USE_NVTX - nvtxRangePush("Calculate inverted matrix"); + nvtxRangePush("Calculate inverted matrix"); #endif - cms(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6); + cms(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6); #ifdef USE_NVTX - nvtxRangePop(); + nvtxRangePop(); #endif - interval_end = chrono::high_resolution_clock::now(); - elapsed = interval_end - interval_start; - message = "INFO: matrix calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n"; - logger->log(message); - interval_start = chrono::high_resolution_clock::now(); + interval_end = chrono::high_resolution_clock::now(); + elapsed = interval_end - interval_start; + message = "INFO: matrix calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n"; + logger->log(message); + interval_start = chrono::high_resolution_clock::now(); #ifdef USE_NVTX - nvtxRangePush("Invert the matrix"); + nvtxRangePush("Invert the matrix"); #endif - invert_matrix(cid->am, ndit, jer, mxndm); + invert_matrix(cid->am, ndit, jer, mxndm); #ifdef USE_NVTX - nvtxRangePop(); + nvtxRangePop(); #endif - interval_end = chrono::high_resolution_clock::now(); - elapsed = interval_end - interval_start; - message = "INFO: matrix inversion for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n"; - logger->log(message); - if (jer != 0) { - message = "ERROR: matrix inversion ended with error code " + to_string(jer) + ".\n"; - logger->err(message); - return jer; - // break; // jxi488 loop: goes to memory clean - } - interval_start = chrono::high_resolution_clock::now(); + interval_end = chrono::high_resolution_clock::now(); + elapsed = interval_end - interval_start; + message = "INFO: matrix inversion for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n"; + logger->log(message); + if (jer != 0) { + message = "ERROR: matrix inversion ended with error code " + to_string(jer) + ".\n"; + logger->err(message); + return jer; + // break; // jxi488 loop: goes to memory clean + } + interval_start = chrono::high_resolution_clock::now(); #ifdef USE_NVTX - nvtxRangePush("Average calculation"); + nvtxRangePush("Average calculation"); #endif - ztm(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6, cid->c9); - if (idfc >= 0) { - if (jxi488 == jwtm) { - int nlemt = 2 * cid->c4->nlem; - string ttms_name = output_path + "/c_TTMS.hd5"; - TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m, "HDF5"); - ttms_name = output_path + "/c_TTMS"; - TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m); - } - } - // label 156: continue from here - if (inpol == 0) { - sprintf(virtual_line, " LIN\n\0"); - output->append_line(virtual_line); - } else { // label 158 - sprintf(virtual_line, " CIRC\n\0"); - output->append_line(virtual_line); + ztm(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6, cid->c9); + if (idfc >= 0) { + if (jxi488 == jwtm) { + int nlemt = 2 * cid->c4->nlem; + string ttms_name = output_path + "/c_TTMS.hd5"; + TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m, "HDF5"); + ttms_name = output_path + "/c_TTMS"; + TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m); } - // label 160 - double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0); - double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0; - dcomplex s0 = 0.0 + 0.0 * I; - scr0(cid->vk, exri, cid->c1, cid->c1ao, cid->c3, cid->c4); - double sqk = cid->vk * cid->vk * exdc; - aps(cid->zpv, cid->c4->li, nsph, cid->c1, sqk, cid->gaps); - rabas(inpol, cid->c4->li, nsph, cid->c1, cid->tqse, cid->tqspe, cid->tqss, cid->tqsps); - if (cid->c4->li != cid->c4->le) { - sprintf(virtual_line, " SPHERES; LMX=LI\n\0"); + } + // label 156: continue from here + if (inpol == 0) { + sprintf(virtual_line, " LIN\n\0"); + output->append_line(virtual_line); + } else { // label 158 + sprintf(virtual_line, " CIRC\n\0"); + output->append_line(virtual_line); + } + // label 160 + double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0); + double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0; + dcomplex s0 = 0.0 + 0.0 * I; + scr0(cid->vk, exri, cid->c1, cid->c1ao, cid->c3, cid->c4); + double sqk = cid->vk * cid->vk * exdc; + aps(cid->zpv, cid->c4->li, nsph, cid->c1, sqk, cid->gaps); + rabas(inpol, cid->c4->li, nsph, cid->c1, cid->tqse, cid->tqspe, cid->tqss, cid->tqsps); + if (cid->c4->li != cid->c4->le) { + sprintf(virtual_line, " SPHERES; LMX=LI\n\0"); + output->append_line(virtual_line); + } + for (int i170 = 1; i170 <= nsph; i170++) { + if (cid->c1->iog[i170 - 1] >= i170) { + int i = i170 - 1; + double albeds = cid->c1->sscs[i] / cid->c1->sexs[i]; + cid->c1->sqscs[i] *= cid->sqsfi; + cid->c1->sqabs[i] *= cid->sqsfi; + cid->c1->sqexs[i] *= cid->sqsfi; + sprintf(virtual_line, " SPHERE %2d\n\0", i170); output->append_line(virtual_line); - } - for (int i170 = 1; i170 <= nsph; i170++) { - if (cid->c1->iog[i170 - 1] >= i170) { - int i = i170 - 1; - double albeds = cid->c1->sscs[i] / cid->c1->sexs[i]; - cid->c1->sqscs[i] *= cid->sqsfi; - cid->c1->sqabs[i] *= cid->sqsfi; - cid->c1->sqexs[i] *= cid->sqsfi; - sprintf(virtual_line, " SPHERE %2d\n\0", i170); - output->append_line(virtual_line); - if (cid->c1->nshl[i] != 1) { - sprintf(virtual_line, " SIZE=%15.7lE\n\0", cid->c2->vsz[i]); - output->append_line(virtual_line); - } else { // label 162 - sprintf(virtual_line, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n\0", cid->c2->vsz[i], real(cid->c2->vkt[i]), imag(cid->c2->vkt[i])); - output->append_line(virtual_line); - } - // label 164 - sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n\0"); - output->append_line(virtual_line); - sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", cid->c1->sscs[i], cid->c1->sabs[i], cid->c1->sexs[i], albeds); - output->append_line(virtual_line); - sprintf(virtual_line, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n\0"); - output->append_line(virtual_line); - sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", cid->c1->sqscs[i], cid->c1->sqabs[i], cid->c1->sqexs[i]); - output->append_line(virtual_line); - sprintf(virtual_line, " FSAS=%15.7lE%15.7lE\n\0", real(cid->c1->fsas[i]), imag(cid->c1->fsas[i])); - output->append_line(virtual_line); - double alamb = 2.0 * 3.141592653589793 / cid->vk; - sprintf(virtual_line, "INSERTION: CS_SPHERE %15.7lE%15.7lE%15.7lE%15.7lE\n\0", alamb, cid->c1->sscs[i], cid->c1->sabs[i], cid->c1->sexs[i]); - output->append_line(virtual_line); - csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcsv[i]; - s0 = cid->c1->fsas[i] * exri; - qschu = imag(s0) * csch; - pschu = real(s0) * csch; - s0mag = cabs(s0) * cs0; - sprintf(virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschu, pschu, s0mag); - output->append_line(virtual_line); - double rapr = cid->c1->sexs[i] - cid->gaps[i]; - double cosav = cid->gaps[i] / cid->c1->sscs[i]; - sprintf(virtual_line, " COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr); + if (cid->c1->nshl[i] != 1) { + sprintf(virtual_line, " SIZE=%15.7lE\n\0", cid->c2->vsz[i]); output->append_line(virtual_line); - sprintf(virtual_line, " IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n\0", cid->tqse[0][i], cid->tqss[0][i]); - output->append_line(virtual_line); - sprintf(virtual_line, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n\0", cid->tqse[1][i], cid->tqss[1][i]); + } else { // label 162 + sprintf(virtual_line, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n\0", cid->c2->vsz[i], real(cid->c2->vkt[i]), imag(cid->c2->vkt[i])); output->append_line(virtual_line); } - } // i170 loop - sprintf(virtual_line, " FSAT=%15.7lE%15.7lE\n\0", real(cid->c3->tfsas), imag(cid->c3->tfsas)); - output->append_line(virtual_line); - csch = 2.0 * cid->vk * cid->sqsfi / cid->c3->gcs; - s0 = cid->c3->tfsas * exri; - qschu = imag(s0) * csch; - pschu = real(s0) * csch; - s0mag = cabs(s0) * cs0; - sprintf(virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschu, pschu, s0mag); - output->append_line(virtual_line); - // tppoan.write(reinterpret_cast<char *>(&(cid->vk)), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(cid->vk)); - pcrsm0(cid->vk, exri, inpol, cid->c1, cid->c1ao, cid->c4); - apcra(cid->zpv, cid->c4->le, cid->c1ao->am0m, inpol, sqk, cid->gapm, cid->gappm); + // label 164 + sprintf(virtual_line, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n\0"); + output->append_line(virtual_line); + sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", cid->c1->sscs[i], cid->c1->sabs[i], cid->c1->sexs[i], albeds); + output->append_line(virtual_line); + sprintf(virtual_line, " ---- SCS/GS -- ABS/GS -- EXS/GS ---\n\0"); + output->append_line(virtual_line); + sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", cid->c1->sqscs[i], cid->c1->sqabs[i], cid->c1->sqexs[i]); + output->append_line(virtual_line); + sprintf(virtual_line, " FSAS=%15.7lE%15.7lE\n\0", real(cid->c1->fsas[i]), imag(cid->c1->fsas[i])); + output->append_line(virtual_line); + double alamb = 2.0 * 3.141592653589793 / cid->vk; + sprintf(virtual_line, "INSERTION: CS_SPHERE %15.7lE%15.7lE%15.7lE%15.7lE\n\0", alamb, cid->c1->sscs[i], cid->c1->sabs[i], cid->c1->sexs[i]); + output->append_line(virtual_line); + csch = 2.0 * cid->vk * cid->sqsfi / cid->c1->gcsv[i]; + s0 = cid->c1->fsas[i] * exri; + qschu = imag(s0) * csch; + pschu = real(s0) * csch; + s0mag = cabs(s0) * cs0; + sprintf(virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschu, pschu, s0mag); + output->append_line(virtual_line); + double rapr = cid->c1->sexs[i] - cid->gaps[i]; + double cosav = cid->gaps[i] / cid->c1->sscs[i]; + sprintf(virtual_line, " COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr); + output->append_line(virtual_line); + sprintf(virtual_line, " IPO= 1, TQEk=%15.7lE, TQSk=%15.7lE\n\0", cid->tqse[0][i], cid->tqss[0][i]); + output->append_line(virtual_line); + sprintf(virtual_line, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n\0", cid->tqse[1][i], cid->tqss[1][i]); + output->append_line(virtual_line); + } + } // i170 loop + sprintf(virtual_line, " FSAT=%15.7lE%15.7lE\n\0", real(cid->c3->tfsas), imag(cid->c3->tfsas)); + output->append_line(virtual_line); + csch = 2.0 * cid->vk * cid->sqsfi / cid->c3->gcs; + s0 = cid->c3->tfsas * exri; + qschu = imag(s0) * csch; + pschu = real(s0) * csch; + s0mag = cabs(s0) * cs0; + sprintf(virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschu, pschu, s0mag); + output->append_line(virtual_line); + // tppoan.write(reinterpret_cast<char *>(&(cid->vk)), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(cid->vk)); + pcrsm0(cid->vk, exri, inpol, cid->c1, cid->c1ao, cid->c4); + apcra(cid->zpv, cid->c4->le, cid->c1ao->am0m, inpol, sqk, cid->gapm, cid->gappm); #ifdef USE_NVTX - nvtxRangePop(); + nvtxRangePop(); #endif - interval_end = chrono::high_resolution_clock::now(); - elapsed = interval_end - interval_start; - message = "INFO: average calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n\0"; - logger->log(message); - interval_start = chrono::high_resolution_clock::now(); + interval_end = chrono::high_resolution_clock::now(); + elapsed = interval_end - interval_start; + message = "INFO: average calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n\0"; + logger->log(message); + interval_start = chrono::high_resolution_clock::now(); #ifdef USE_NVTX - nvtxRangePush("Angle loop"); + nvtxRangePush("Angle loop"); #endif - double th = sa->th; - for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable? - double ph = sa->ph; - double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0; - for (int jph484 = 1; jph484 <= sa->nph; jph484++) { - int jw = 0; - if (sa->nk != 1 || jxi488 <= 1) { - upvmp(th, ph, 0, cost, sint, cosp, sinp, cid->u, cid->upmp, cid->unmp); - if (isam >= 0) { - wmamp( - 0, cost, sint, cosp, sinp, inpol, cid->c4->le, 0, - nsph, cid->argi, cid->u, cid->upmp, cid->unmp, cid->c1 + double th = sa->th; + for (int jth486 = 1; jth486 <= sa->nth; jth486++) { // OpenMP portable? + double ph = sa->ph; + double cost = 0.0, sint = 0.0, cosp = 0.0, sinp = 0.0; + for (int jph484 = 1; jph484 <= sa->nph; jph484++) { + int jw = 0; + if (sa->nk != 1 || jxi488 <= 1) { + upvmp(th, ph, 0, cost, sint, cosp, sinp, cid->u, cid->upmp, cid->unmp); + if (isam >= 0) { + wmamp( + 0, cost, sint, cosp, sinp, inpol, cid->c4->le, 0, + nsph, cid->argi, cid->u, cid->upmp, cid->unmp, cid->c1 + ); + // label 182 + apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp); + raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps); + jw = 1; + } + } else { // label 180, NK == 1 AND JXI488 == 1 + if (isam >= 0) { + // label 182 + apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp); + raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps); + jw = 1; + } + } + // label 184 + double thsl = sa->ths; + double phsph = 0.0; + for (int jths = 1; jths <= sa->nths; jths++) { + double ths = thsl; + int icspnv = 0; + if (isam > 1) ths += sa->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 + double phs = sa->phs; + for (int jphs = 1; jphs <= sa->nphs; jphs++) { + double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0; + if (isam >= 1) { + phs = sa->ph + phsph; + if (phs > 360.0) phs -= 360.0; + } + // label 188 + bool goto190 = (sa->nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1)); + if (!goto190) { + upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, cid->us, cid->upsmp, cid->unsmp); + if (isam >= 0) + wmamp( + 2, costs, sints, cosps, sinps, inpol, cid->c4->le, + 0, nsph, cid->args, cid->us, cid->upsmp, cid->unsmp, cid->c1 + ); + } + // label 190 + if (sa->nkks != 1 || jxi488 <= 1) { + upvsp( + cid->u, cid->upmp, cid->unmp, cid->us, cid->upsmp, cid->unsmp, cid->up, cid->un, cid->ups, cid->uns, + cid->duk, isq, ibf, cid->scan, cid->cfmp, cid->sfmp, cid->cfsp, cid->sfsp ); - // label 182 - apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp); - raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps); - jw = 1; + if (isam < 0) { + wmasp( + cost, sint, cosp, sinp, costs, sints, cosps, sinps, + cid->u, cid->up, cid->un, cid->us, cid->ups, cid->uns, isq, ibf, inpol, cid->c4->le, + 0, nsph, cid->argi, cid->args, cid->c1 + ); + } else { // label 192 + for (int i193 = 0; i193 < 3; i193++) { + cid->up[i193] = cid->upmp[i193]; + cid->un[i193] = cid->unmp[i193]; + cid->ups[i193] = cid->upsmp[i193]; + cid->uns[i193] = cid->unsmp[i193]; + } + } } - } else { // label 180, NK == 1 AND JXI488 == 1 - if (isam >= 0) { - // label 182 + // label 194 + if (iavm == 1) crsm1(cid->vk, exri, cid->c1, cid->c1ao, cid->c4, cid->c6); + if (isam < 0) { apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp); raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps); jw = 1; } - } - // label 184 - double thsl = sa->ths; - double phsph = 0.0; - for (int jths = 1; jths <= sa->nths; jths++) { - double ths = thsl; - int icspnv = 0; - if (isam > 1) ths += sa->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 - double phs = sa->phs; - for (int jphs = 1; jphs <= sa->nphs; jphs++) { - double costs = 0.0, sints = 0.0, cosps = 0.0, sinps = 0.0; - if (isam >= 1) { - phs = sa->ph + phsph; - if (phs > 360.0) phs -= 360.0; - } - // label 188 - bool goto190 = (sa->nks == 1 && (jxi488 > 1 || jth486 > 1 || jph484 > 1)); - if (!goto190) { - upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, cid->us, cid->upsmp, cid->unsmp); - if (isam >= 0) - wmamp( - 2, costs, sints, cosps, sinps, inpol, cid->c4->le, - 0, nsph, cid->args, cid->us, cid->upsmp, cid->unsmp, cid->c1 - ); - } - // label 190 - if (sa->nkks != 1 || jxi488 <= 1) { - upvsp( - cid->u, cid->upmp, cid->unmp, cid->us, cid->upsmp, cid->unsmp, cid->up, cid->un, cid->ups, cid->uns, - cid->duk, isq, ibf, cid->scan, cid->cfmp, cid->sfmp, cid->cfsp, cid->sfsp - ); - if (isam < 0) { - wmasp( - cost, sint, cosp, sinp, costs, sints, cosps, sinps, - cid->u, cid->up, cid->un, cid->us, cid->ups, cid->uns, isq, ibf, inpol, cid->c4->le, - 0, nsph, cid->argi, cid->args, cid->c1 - ); - } else { // label 192 - for (int i193 = 0; i193 < 3; i193++) { - cid->up[i193] = cid->upmp[i193]; - cid->un[i193] = cid->unmp[i193]; - cid->ups[i193] = cid->upsmp[i193]; - cid->uns[i193] = cid->unsmp[i193]; - } + // label 196 + // tppoan.write(reinterpret_cast<char *>(&th), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(th)); + // tppoan.write(reinterpret_cast<char *>(&ph), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(ph)); + // tppoan.write(reinterpret_cast<char *>(&ths), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(ths)); + // tppoan.write(reinterpret_cast<char *>(&phs), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(phs)); + // tppoan.write(reinterpret_cast<char *>(&(cid->scan)), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(cid->scan)); + if (jaw != 0) { + jaw = 0; + mextc(cid->vk, exri, cid->c1ao->fsacm, cid->cextlr, cid->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 = cid->cext[i][j]; + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); } } - // label 194 - if (iavm == 1) crsm1(cid->vk, exri, cid->c1, cid->c1ao, cid->c4, cid->c6); - if (isam < 0) { - apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp); - raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps); - jw = 1; + for (int i = 0; i < 2; i++) { + double value = cid->c1ao->scscm[i]; + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = real(cid->c1ao->scscpm[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = imag(cid->c1ao->scscpm[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = cid->c1ao->ecscm[i]; + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = real(cid->c1ao->ecscpm[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = imag(cid->c1ao->ecscpm[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); } - // label 196 - // tppoan.write(reinterpret_cast<char *>(&th), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(th)); - // tppoan.write(reinterpret_cast<char *>(&ph), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(ph)); - // tppoan.write(reinterpret_cast<char *>(&ths), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(ths)); - // tppoan.write(reinterpret_cast<char *>(&phs), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(phs)); - // tppoan.write(reinterpret_cast<char *>(&(cid->scan)), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(cid->scan)); - if (jaw != 0) { - jaw = 0; - mextc(cid->vk, exri, cid->c1ao->fsacm, cid->cextlr, cid->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 = cid->cext[i][j]; - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - } - } - for (int i = 0; i < 2; i++) { - double value = cid->c1ao->scscm[i]; - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = real(cid->c1ao->scscpm[i]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = imag(cid->c1ao->scscpm[i]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = cid->c1ao->ecscm[i]; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 2; j++) { + double value = cid->gapm[i][j]; // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); - value = real(cid->c1ao->ecscpm[i]); + value = real(cid->gappm[i][j]); // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); - value = imag(cid->c1ao->ecscpm[i]); + value = imag(cid->gappm[i][j]); // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); } - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 2; j++) { - double value = cid->gapm[i][j]; - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = real(cid->gappm[i][j]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = imag(cid->gappm[i][j]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - } - } - sprintf(virtual_line, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n\0", iavm); - output->append_line(virtual_line); - int jlr = 2; - for (int ilr210 = 1; ilr210 <= 2; ilr210++) { - int ipol = (ilr210 % 2 == 0) ? 1 : -1; - if (ilr210 == 2) jlr = 1; - double extsm = cid->c1ao->ecscm[ilr210 - 1]; - double qextm = extsm * cid->sqsfi / cid->c3->gcs; - double extrm = extsm / cid->c3->ecs; - double scasm = cid->c1ao->scscm[ilr210 - 1]; - double albdm = scasm / extsm; - double qscam = scasm * cid->sqsfi / cid->c3->gcs; - double scarm = scasm / cid->c3->scs; - double abssm = extsm - scasm; - double qabsm = abssm * cid->sqsfi / cid->c3->gcs; - double absrm = abssm / cid->c3->acs; - double acsecs = cid->c3->acs / cid->c3->ecs; - if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0; - dcomplex s0m = cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri; - double qschum = imag(s0m) * csch; - double pschum = real(s0m) * csch; - double s0magm = cabs(s0m) * cs0; - double rfinrm = real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / real(cid->c3->tfsas); - double extcrm = imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / imag(cid->c3->tfsas); - if (inpol == 0) { - sprintf(virtual_line, " LIN %2d\n\0", ipol); - output->append_line(virtual_line); - } else { // label 206 - sprintf(virtual_line, " CIRC %2d\n\0", ipol); - output->append_line(virtual_line); - } - // label 208 - sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n\0"); - output->append_line(virtual_line); - sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", scasm, abssm, extsm, albdm); - output->append_line(virtual_line); - sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n\0"); - output->append_line(virtual_line); - sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", qscam, qabsm, qextm); - output->append_line(virtual_line); - sprintf(virtual_line, " ---- SCCRT --- ABCRT --- EXCRT ----\n\0"); - output->append_line(virtual_line); - sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", scarm, absrm, extrm); - output->append_line(virtual_line); - sprintf( - virtual_line, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n\0", - ilr210, ilr210, real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]), - imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210, - real(cid->c1ao->fsacm[jlr - 1][ilr210 - 1]), imag(cid->c1ao->fsacm[jlr - 1][ilr210 - 1]) - ); - output->append_line(virtual_line); - sprintf( - virtual_line, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n\0", - ilr210, ilr210, rfinrm, ilr210, ilr210, extcrm - ); - output->append_line(virtual_line); - sprintf(virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschum, pschum, s0magm); - output->append_line(virtual_line); - double rapr = cid->c1ao->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1]; - double cosav = cid->gapm[2][ilr210 - 1] / cid->c1ao->scscm[ilr210 - 1]; - double fz = rapr; - sprintf(virtual_line, " COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr); - output->append_line(virtual_line); - sprintf(virtual_line, " Fk=%15.7lE\n\0", fz); + } + sprintf(virtual_line, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n\0", iavm); + output->append_line(virtual_line); + int jlr = 2; + for (int ilr210 = 1; ilr210 <= 2; ilr210++) { + int ipol = (ilr210 % 2 == 0) ? 1 : -1; + if (ilr210 == 2) jlr = 1; + double extsm = cid->c1ao->ecscm[ilr210 - 1]; + double qextm = extsm * cid->sqsfi / cid->c3->gcs; + double extrm = extsm / cid->c3->ecs; + double scasm = cid->c1ao->scscm[ilr210 - 1]; + double albdm = scasm / extsm; + double qscam = scasm * cid->sqsfi / cid->c3->gcs; + double scarm = scasm / cid->c3->scs; + double abssm = extsm - scasm; + double qabsm = abssm * cid->sqsfi / cid->c3->gcs; + double absrm = abssm / cid->c3->acs; + double acsecs = cid->c3->acs / cid->c3->ecs; + if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0; + dcomplex s0m = cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri; + double qschum = imag(s0m) * csch; + double pschum = real(s0m) * csch; + double s0magm = cabs(s0m) * cs0; + double rfinrm = real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / real(cid->c3->tfsas); + double extcrm = imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / imag(cid->c3->tfsas); + if (inpol == 0) { + sprintf(virtual_line, " LIN %2d\n\0", ipol); output->append_line(virtual_line); - double alamb = 2.0 * 3.141592653589793 / cid->vk; - sprintf(virtual_line, "INSERTION: CSM_CLUSTER %15.7lE%15.7lE%15.7lE%15.7lE\n\0", alamb, scasm, abssm, extsm); + } else { // label 206 + sprintf(virtual_line, " CIRC %2d\n\0", ipol); output->append_line(virtual_line); - } // ilr210 loop - double rmbrif = (real(cid->c1ao->fsacm[0][0]) - real(cid->c1ao->fsacm[1][1])) / real(cid->c1ao->fsacm[0][0]); - double rmdchr = (imag(cid->c1ao->fsacm[0][0]) - imag(cid->c1ao->fsacm[1][1])) / imag(cid->c1ao->fsacm[0][0]); - sprintf(virtual_line, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n\0", rmbrif); + } + // label 208 + sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n\0"); output->append_line(virtual_line); - sprintf(virtual_line, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n\0", rmdchr); + sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", scasm, abssm, extsm, albdm); output->append_line(virtual_line); - } - // label 212 - sprintf(virtual_line, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n\0", jth486, jph484, jths, jphs); + sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n\0"); + output->append_line(virtual_line); + sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", qscam, qabsm, qextm); + output->append_line(virtual_line); + sprintf(virtual_line, " ---- SCCRT --- ABCRT --- EXCRT ----\n\0"); + output->append_line(virtual_line); + sprintf(virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", scarm, absrm, extrm); + output->append_line(virtual_line); + sprintf( + virtual_line, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n\0", + ilr210, ilr210, real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]), + imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210, + real(cid->c1ao->fsacm[jlr - 1][ilr210 - 1]), imag(cid->c1ao->fsacm[jlr - 1][ilr210 - 1]) + ); + output->append_line(virtual_line); + sprintf( + virtual_line, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n\0", + ilr210, ilr210, rfinrm, ilr210, ilr210, extcrm + ); + output->append_line(virtual_line); + sprintf(virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", qschum, pschum, s0magm); + output->append_line(virtual_line); + double rapr = cid->c1ao->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1]; + double cosav = cid->gapm[2][ilr210 - 1] / cid->c1ao->scscm[ilr210 - 1]; + double fz = rapr; + sprintf(virtual_line, " COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr); + output->append_line(virtual_line); + sprintf(virtual_line, " Fk=%15.7lE\n\0", fz); + output->append_line(virtual_line); + double alamb = 2.0 * 3.141592653589793 / cid->vk; + sprintf(virtual_line, "INSERTION: CSM_CLUSTER %15.7lE%15.7lE%15.7lE%15.7lE\n\0", alamb, scasm, abssm, extsm); + output->append_line(virtual_line); + } // ilr210 loop + double rmbrif = (real(cid->c1ao->fsacm[0][0]) - real(cid->c1ao->fsacm[1][1])) / real(cid->c1ao->fsacm[0][0]); + double rmdchr = (imag(cid->c1ao->fsacm[0][0]) - imag(cid->c1ao->fsacm[1][1])) / imag(cid->c1ao->fsacm[0][0]); + sprintf(virtual_line, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n\0", rmbrif); output->append_line(virtual_line); - sprintf(virtual_line, " TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n\0", th, ph, ths, phs); + sprintf(virtual_line, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n\0", rmdchr); output->append_line(virtual_line); - sprintf(virtual_line, " SCAND=%10.3lE\n\0", cid->scan); + } + // label 212 + sprintf(virtual_line, "********** JTH =%3d, JPH =%3d, JTHS =%3d, JPHS =%3d ********************\n\0", jth486, jph484, jths, jphs); + output->append_line(virtual_line); + sprintf(virtual_line, " TIDG=%10.3lE, PIDG=%10.3lE, TSDG=%10.3lE, PSDG=%10.3lE\n\0", th, ph, ths, phs); + output->append_line(virtual_line); + sprintf(virtual_line, " SCAND=%10.3lE\n\0", cid->scan); + output->append_line(virtual_line); + sprintf(virtual_line, " CFMP=%15.7lE, SFMP=%15.7lE\n\0", cid->cfmp, cid->sfmp); + output->append_line(virtual_line); + sprintf(virtual_line, " CFSP=%15.7lE, SFSP=%15.7lE\n\0", cid->cfsp, cid->sfsp); + output->append_line(virtual_line); + if (isam >= 0) { + sprintf(virtual_line, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n\0", cid->un[0], cid->un[1], cid->un[2]); output->append_line(virtual_line); - sprintf(virtual_line, " CFMP=%15.7lE, SFMP=%15.7lE\n\0", cid->cfmp, cid->sfmp); + sprintf(virtual_line, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n\0", cid->uns[0], cid->uns[1], cid->uns[2]); output->append_line(virtual_line); - sprintf(virtual_line, " CFSP=%15.7lE, SFSP=%15.7lE\n\0", cid->cfsp, cid->sfsp); + } else { // label 214 + sprintf(virtual_line, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n\0", cid->un[0], cid->un[1], cid->un[2]); output->append_line(virtual_line); - if (isam >= 0) { - sprintf(virtual_line, " UNI=(%12.5lE,%12.5lE,%12.5lE)\n\0", cid->un[0], cid->un[1], cid->un[2]); - output->append_line(virtual_line); - sprintf(virtual_line, " UNS=(%12.5lE,%12.5lE,%12.5lE)\n\0", cid->uns[0], cid->uns[1], cid->uns[2]); - output->append_line(virtual_line); - } else { // label 214 - sprintf(virtual_line, " UN=(%12.5lE,%12.5lE,%12.5lE)\n\n\0", cid->un[0], cid->un[1], cid->un[2]); + } + // label 220 + if (inpol == 0) { + sprintf(virtual_line, " LIN\n\0"); + output->append_line(virtual_line); + } else { // label 222 + sprintf(virtual_line, " CIRC\n\0"); + output->append_line(virtual_line); + } + // label 224 + scr2(cid->vk, vkarg, exri, cid->duk, cid->c1, cid->c1ao, cid->c3, cid->c4); + if (cid->c4->li != cid->c4->le) { + sprintf(virtual_line, " SPHERES; LMX=MIN0(LI,LE)\n\0"); + output->append_line(virtual_line); + } + for (int i226 = 1; i226 <= nsph; i226++) { + if (cid->c1->iog[i226 - 1] >= i226) { + sprintf(virtual_line, " SPHERE %2d\n\0", i226); output->append_line(virtual_line); - } - // label 220 - if (inpol == 0) { - sprintf(virtual_line, " LIN\n\0"); + sprintf( + virtual_line, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n\0", + real(cid->c1->sas[i226 - 1][0][0]), imag(cid->c1->sas[i226 - 1][0][0]), + real(cid->c1->sas[i226 - 1][1][0]), imag(cid->c1->sas[i226 - 1][1][0]) + ); output->append_line(virtual_line); - } else { // label 222 - sprintf(virtual_line, " CIRC\n\0"); + sprintf( + virtual_line, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n\0", + real(cid->c1->sas[i226 - 1][0][1]), imag(cid->c1->sas[i226 - 1][0][1]), + real(cid->c1->sas[i226 - 1][1][1]), imag(cid->c1->sas[i226 - 1][1][1]) + ); output->append_line(virtual_line); - } - // label 224 - scr2(cid->vk, vkarg, exri, cid->duk, cid->c1, cid->c1ao, cid->c3, cid->c4); - if (cid->c4->li != cid->c4->le) { - sprintf(virtual_line, " SPHERES; LMX=MIN0(LI,LE)\n\0"); + for (int j225 = 0; j225 < 16; j225++) { + cid->c1->vint[j225] = cid->c1->vints[i226 - 1][j225]; + } // j225 loop + mmulc(cid->c1->vint, cid->cmullr, cid->cmul); + sprintf(virtual_line, " MULS\n\0"); output->append_line(virtual_line); - } - for (int i226 = 1; i226 <= nsph; i226++) { - if (cid->c1->iog[i226 - 1] >= i226) { - sprintf(virtual_line, " SPHERE %2d\n\0", i226); - output->append_line(virtual_line); + for (int i1 = 0; i1 < 4; i1++) { sprintf( - virtual_line, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n\0", - real(cid->c1->sas[i226 - 1][0][0]), imag(cid->c1->sas[i226 - 1][0][0]), - real(cid->c1->sas[i226 - 1][1][0]), imag(cid->c1->sas[i226 - 1][1][0]) + virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n\0", + cid->cmul[i1][0], cid->cmul[i1][1], cid->cmul[i1][2], cid->cmul[i1][3] ); output->append_line(virtual_line); + } // i1 loop + sprintf(virtual_line, " MULSLR\n\0"); + output->append_line(virtual_line); + for (int i1 = 0; i1 < 4; i1++) { sprintf( - virtual_line, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n\0", - real(cid->c1->sas[i226 - 1][0][1]), imag(cid->c1->sas[i226 - 1][0][1]), - real(cid->c1->sas[i226 - 1][1][1]), imag(cid->c1->sas[i226 - 1][1][1]) + virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n\0", + cid->cmullr[i1][0], cid->cmullr[i1][1], cid->cmullr[i1][2], cid->cmullr[i1][3] ); output->append_line(virtual_line); - for (int j225 = 0; j225 < 16; j225++) { - cid->c1->vint[j225] = cid->c1->vints[i226 - 1][j225]; - } // j225 loop - mmulc(cid->c1->vint, cid->cmullr, cid->cmul); - sprintf(virtual_line, " MULS\n\0"); - output->append_line(virtual_line); - for (int i1 = 0; i1 < 4; i1++) { - sprintf( - virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n\0", - cid->cmul[i1][0], cid->cmul[i1][1], cid->cmul[i1][2], cid->cmul[i1][3] - ); - output->append_line(virtual_line); - } // i1 loop - sprintf(virtual_line, " MULSLR\n\0"); - output->append_line(virtual_line); - for (int i1 = 0; i1 < 4; i1++) { - sprintf( - virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n\0", - cid->cmullr[i1][0], cid->cmullr[i1][1], cid->cmullr[i1][2], cid->cmullr[i1][3] - ); - output->append_line(virtual_line); - } // i1 loop - } - } // i226 loop - sprintf( - virtual_line, " SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n\0", - real(cid->c3->tsas[0][0]), imag(cid->c3->tsas[0][0]), - real(cid->c3->tsas[1][0]), imag(cid->c3->tsas[1][0]) - ); - output->append_line(virtual_line); - sprintf( - virtual_line, " SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n\0", - real(cid->c3->tsas[0][1]), imag(cid->c3->tsas[0][1]), - real(cid->c3->tsas[1][1]), imag(cid->c3->tsas[1][1]) - ); - output->append_line(virtual_line); - sprintf(virtual_line, " CLUSTER\n\0"); - output->append_line(virtual_line); - pcros(cid->vk, exri, cid->c1, cid->c1ao, cid->c4); - mextc(cid->vk, exri, cid->c1ao->fsac, cid->cextlr, cid->cext); - mmulc(cid->c1->vint, cid->cmullr, cid->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 = cid->cext[i][j]; - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - } + } // i1 loop + } + } // i226 loop + sprintf( + virtual_line, " SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n\0", + real(cid->c3->tsas[0][0]), imag(cid->c3->tsas[0][0]), + real(cid->c3->tsas[1][0]), imag(cid->c3->tsas[1][0]) + ); + output->append_line(virtual_line); + sprintf( + virtual_line, " SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n\0", + real(cid->c3->tsas[0][1]), imag(cid->c3->tsas[0][1]), + real(cid->c3->tsas[1][1]), imag(cid->c3->tsas[1][1]) + ); + output->append_line(virtual_line); + sprintf(virtual_line, " CLUSTER\n\0"); + output->append_line(virtual_line); + pcros(cid->vk, exri, cid->c1, cid->c1ao, cid->c4); + mextc(cid->vk, exri, cid->c1ao->fsac, cid->cextlr, cid->cext); + mmulc(cid->c1->vint, cid->cmullr, cid->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 = cid->cext[i][j]; + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); } - for (int i = 0; i < 2; i++) { - double value = cid->c1ao->scsc[i]; + } + for (int i = 0; i < 2; i++) { + double value = cid->c1ao->scsc[i]; + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = real(cid->c1ao->scscp[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = imag(cid->c1ao->scscp[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = cid->c1ao->ecsc[i]; + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = real(cid->c1ao->ecscp[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = imag(cid->c1ao->ecscp[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + } + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 2; j++) { + double value = cid->gap[i][j]; // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); - value = real(cid->c1ao->scscp[i]); + value = real(cid->gapp[i][j]); // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); - value = imag(cid->c1ao->scscp[i]); + value = imag(cid->gapp[i][j]); // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); - value = cid->c1ao->ecsc[i]; + } + } + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 3; j++) { + double value = cid->tqce[i][j]; // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); - value = real(cid->c1ao->ecscp[i]); + value = real(cid->tqcpe[i][j]); // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); - value = imag(cid->c1ao->ecscp[i]); + value = imag(cid->tqcpe[i][j]); // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); } - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 2; j++) { - double value = cid->gap[i][j]; - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = real(cid->gapp[i][j]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = imag(cid->gapp[i][j]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - } - } - for (int i = 0; i < 2; i++) { - for (int j = 0; j < 3; j++) { - double value = cid->tqce[i][j]; - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = real(cid->tqcpe[i][j]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = imag(cid->tqcpe[i][j]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - } - } - for (int i = 0; i < 2; i++) { - for (int j = 0; j < 3; j++) { - double value = cid->tqcs[i][j]; - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = real(cid->tqcps[i][j]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = imag(cid->tqcps[i][j]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - } - } - for (int i = 0; i < 3; i++) { - double value = cid->u[i]; + } + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 3; j++) { + double value = cid->tqcs[i][j]; // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); - value = cid->up[i]; + value = real(cid->tqcps[i][j]); // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); - value = cid->un[i]; + value = imag(cid->tqcps[i][j]); // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); } } - // label 254 - for (int i = 0; i < 16; i++) { - double value = real(cid->c1->vint[i]); + for (int i = 0; i < 3; i++) { + double value = cid->u[i]; + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = cid->up[i]; // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); - value = imag(cid->c1->vint[i]); + value = cid->un[i]; // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); vtppoanp->append_line(VirtualBinaryLine(value)); } - for (int i = 0; i < 4; i++) { - for (int j = 0; j < 4; j++) { - double value = cid->cmul[i][j]; - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - } + } + // label 254 + for (int i = 0; i < 16; i++) { + double value = real(cid->c1->vint[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = imag(cid->c1->vint[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + } + for (int i = 0; i < 4; i++) { + for (int j = 0; j < 4; j++) { + double value = cid->cmul[i][j]; + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); } - int jlr = 2; - for (int ilr290 = 1; ilr290 <= 2; ilr290++) { - int ipol = (ilr290 % 2 == 0) ? 1 : -1; - if (ilr290 == 2) jlr = 1; - double extsec = cid->c1ao->ecsc[ilr290 - 1]; - double qext = extsec * cid->sqsfi / cid->c3->gcs; - double extrat = extsec / cid->c3->ecs; - double scasec = cid->c1ao->scsc[ilr290 - 1]; - double albedc = scasec / extsec; - double qsca = scasec * cid->sqsfi / cid->c3->gcs; - double scarat = scasec / cid->c3->scs; - double abssec = extsec - scasec; - double qabs = abssec * cid->sqsfi / cid->c3->gcs; - double absrat = 1.0; - double ratio = cid->c3->acs / cid->c3->ecs; - if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / cid->c3->acs; - s0 = cid->c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri; - double qschu = imag(s0) * csch; - double pschu = real(s0) * csch; - s0mag = cabs(s0) * cs0; - double refinr = real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / real(cid->c3->tfsas); - double extcor = imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / imag(cid->c3->tfsas); - if (inpol == 0) { - sprintf(virtual_line, " LIN %2d\n\0", ipol); - output->append_line(virtual_line); - } else { // label 273 - sprintf(virtual_line, " CIRC %2d\n\0", ipol); - output->append_line(virtual_line); - } - // label 275 - sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n\0"); - output->append_line(virtual_line); - sprintf( - virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", - scasec, abssec, extsec, albedc - ); + } + int jlr = 2; + for (int ilr290 = 1; ilr290 <= 2; ilr290++) { + int ipol = (ilr290 % 2 == 0) ? 1 : -1; + if (ilr290 == 2) jlr = 1; + double extsec = cid->c1ao->ecsc[ilr290 - 1]; + double qext = extsec * cid->sqsfi / cid->c3->gcs; + double extrat = extsec / cid->c3->ecs; + double scasec = cid->c1ao->scsc[ilr290 - 1]; + double albedc = scasec / extsec; + double qsca = scasec * cid->sqsfi / cid->c3->gcs; + double scarat = scasec / cid->c3->scs; + double abssec = extsec - scasec; + double qabs = abssec * cid->sqsfi / cid->c3->gcs; + double absrat = 1.0; + double ratio = cid->c3->acs / cid->c3->ecs; + if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / cid->c3->acs; + s0 = cid->c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri; + double qschu = imag(s0) * csch; + double pschu = real(s0) * csch; + s0mag = cabs(s0) * cs0; + double refinr = real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / real(cid->c3->tfsas); + double extcor = imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / imag(cid->c3->tfsas); + if (inpol == 0) { + sprintf(virtual_line, " LIN %2d\n\0", ipol); output->append_line(virtual_line); - sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n\0"); + } else { // label 273 + sprintf(virtual_line, " CIRC %2d\n\0", ipol); output->append_line(virtual_line); - sprintf( - virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", - qsca, qabs, qext - ); + } + // label 275 + sprintf(virtual_line, " ----- SCC ----- ABC ----- EXC ----- ALBEDC --\n\0"); + output->append_line(virtual_line); + sprintf( + virtual_line, " %14.7lE%15.7lE%15.7lE%15.7lE\n\0", + scasec, abssec, extsec, albedc + ); + output->append_line(virtual_line); + sprintf(virtual_line, " --- SCC/TGS - ABC/TGS - EXC/TGS ---\n\0"); + output->append_line(virtual_line); + sprintf( + virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", + qsca, qabs, qext + ); + output->append_line(virtual_line); + sprintf(virtual_line, " ---- SCCRT --- ABCRT --- EXCRT ----\n\0"); + output->append_line(virtual_line); + sprintf( + virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", + scarat, absrat, extrat + ); + output->append_line(virtual_line); + sprintf( + virtual_line, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n\0", + ilr290, ilr290, real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]), + jlr, ilr290, real(cid->c1ao->fsac[jlr - 1][ilr290 - 1]), imag(cid->c1ao->fsac[jlr - 1][ilr290 - 1]) + ); + output->append_line(virtual_line); + sprintf( + virtual_line, " SAC(%1d,%1d)=%15.7lE%15.7lE SAC(%1d,%1d)=%15.7lE%15.7lE\n\0", + ilr290, ilr290, real(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]), + jlr, ilr290, real(cid->c1ao->sac[jlr - 1][ilr290 - 1]), imag(cid->c1ao->sac[jlr - 1][ilr290 - 1]) + ); + output->append_line(virtual_line); + sprintf( + virtual_line, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n\0", + ilr290, ilr290, refinr, ilr290, ilr290, extcor + ); + output->append_line(virtual_line); + sprintf( + virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", + qschu, pschu, s0mag + ); + output->append_line(virtual_line); + double alamb = 2.0 * 3.141592653589793 / cid->vk; + if (ilr290 == 1) { + sprintf(virtual_line, "INSERTION: CS1_CLUSTER %9.2lf%9.2lf%9.2lf%9.2lf%9.2lf%13.5lf%13.5lf%13.5lf\n\0", alamb, th, ph, ths, phs, scasec, abssec, extsec); + } else if (ilr290 == 2) { + sprintf(virtual_line, "INSERTION: CS2_CLUSTER %9.2lf%9.2lf%9.2lf%9.2lf%9.2lf%13.5lf%13.5lf%13.5lf\n\0", alamb, th, ph, ths, phs, scasec, abssec, extsec); + } + output->append_line(virtual_line); + bool goto190 = isam >= 0 && (jths > 1 || jphs > 1); + if (!goto190) { + cid->gapv[0] = cid->gap[0][ilr290 - 1]; + cid->gapv[1] = cid->gap[1][ilr290 - 1]; + cid->gapv[2] = cid->gap[2][ilr290 - 1]; + double extins = cid->c1ao->ecsc[ilr290 - 1]; + double scatts = cid->c1ao->scsc[ilr290 - 1]; + double rapr, cosav, fp, fn, fk, fx, fy, fz; + rftr(cid->u, cid->up, cid->un, cid->gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz); + sprintf(virtual_line, " COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr); output->append_line(virtual_line); - sprintf(virtual_line, " ---- SCCRT --- ABCRT --- EXCRT ----\n\0"); + sprintf(virtual_line, " Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n\0", fp, fn, fk); output->append_line(virtual_line); - sprintf( - virtual_line, " %14.7lE%15.7lE%15.7lE\n\0", - scarat, absrat, extrat - ); + sprintf(virtual_line, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n\0", fx, fy, fz); output->append_line(virtual_line); - sprintf( - virtual_line, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n\0", - ilr290, ilr290, real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]), - jlr, ilr290, real(cid->c1ao->fsac[jlr - 1][ilr290 - 1]), imag(cid->c1ao->fsac[jlr - 1][ilr290 - 1]) - ); + cid->tqev[0] = cid->tqce[ilr290 - 1][0]; + cid->tqev[1] = cid->tqce[ilr290 - 1][1]; + cid->tqev[2] = cid->tqce[ilr290 - 1][2]; + cid->tqsv[0] = cid->tqcs[ilr290 - 1][0]; + cid->tqsv[1] = cid->tqcs[ilr290 - 1][1]; + cid->tqsv[2] = cid->tqcs[ilr290 - 1][2]; + double tep, ten, tek, tsp, tsn, tsk; + tqr(cid->u, cid->up, cid->un, cid->tqev, cid->tqsv, tep, ten, tek, tsp, tsn, tsk); + sprintf(virtual_line, " TQEl=%15.7lE, TQEr=%15.7lE, TQEk=%15.7lE\n\0", tep, ten, tek); output->append_line(virtual_line); - sprintf( - virtual_line, " SAC(%1d,%1d)=%15.7lE%15.7lE SAC(%1d,%1d)=%15.7lE%15.7lE\n\0", - ilr290, ilr290, real(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]), - jlr, ilr290, real(cid->c1ao->sac[jlr - 1][ilr290 - 1]), imag(cid->c1ao->sac[jlr - 1][ilr290 - 1]) - ); + sprintf(virtual_line, " TQSl=%15.7lE, TQSr=%15.7lE, TQSk=%15.7lE\n\0", tsp, tsn, tsk); output->append_line(virtual_line); sprintf( - virtual_line, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n\0", - ilr290, ilr290, refinr, ilr290, ilr290, extcor + virtual_line, " TQEx=%15.7lE, TQEy=%15.7lE, TQEz=%15.7lE\n\0", + cid->tqce[ilr290 - 1][0], cid->tqce[ilr290 - 1][1], cid->tqce[ilr290 - 1][2] ); output->append_line(virtual_line); sprintf( - virtual_line, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n\0", - qschu, pschu, s0mag + virtual_line, " TQSx=%15.7lE, TQSy=%15.7lE, TQSz=%15.7lE\n\0", + cid->tqcs[ilr290 - 1][0], cid->tqcs[ilr290 - 1][1], cid->tqcs[ilr290 - 1][2] ); output->append_line(virtual_line); - double alamb = 2.0 * 3.141592653589793 / cid->vk; - if (ilr290 == 1) { - sprintf(virtual_line, "INSERTION: CS1_CLUSTER %9.2lf%9.2lf%9.2lf%9.2lf%9.2lf%13.5lf%13.5lf%13.5lf\n\0", alamb, th, ph, ths, phs, scasec, abssec, extsec); - } else if (ilr290 == 2) { - sprintf(virtual_line, "INSERTION: CS2_CLUSTER %9.2lf%9.2lf%9.2lf%9.2lf%9.2lf%13.5lf%13.5lf%13.5lf\n\0", alamb, th, ph, ths, phs, scasec, abssec, extsec); - } - output->append_line(virtual_line); - bool goto190 = isam >= 0 && (jths > 1 || jphs > 1); - if (!goto190) { - cid->gapv[0] = cid->gap[0][ilr290 - 1]; - cid->gapv[1] = cid->gap[1][ilr290 - 1]; - cid->gapv[2] = cid->gap[2][ilr290 - 1]; - double extins = cid->c1ao->ecsc[ilr290 - 1]; - double scatts = cid->c1ao->scsc[ilr290 - 1]; - double rapr, cosav, fp, fn, fk, fx, fy, fz; - rftr(cid->u, cid->up, cid->un, cid->gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz); - sprintf(virtual_line, " COSAV=%15.7lE, RAPRS=%15.7lE\n\0", cosav, rapr); - output->append_line(virtual_line); - sprintf(virtual_line, " Fl=%15.7lE, Fr=%15.7lE, Fk=%15.7lE\n\0", fp, fn, fk); - output->append_line(virtual_line); - sprintf(virtual_line, " Fx=%15.7lE, Fy=%15.7lE, Fz=%15.7lE\n\0", fx, fy, fz); - output->append_line(virtual_line); - cid->tqev[0] = cid->tqce[ilr290 - 1][0]; - cid->tqev[1] = cid->tqce[ilr290 - 1][1]; - cid->tqev[2] = cid->tqce[ilr290 - 1][2]; - cid->tqsv[0] = cid->tqcs[ilr290 - 1][0]; - cid->tqsv[1] = cid->tqcs[ilr290 - 1][1]; - cid->tqsv[2] = cid->tqcs[ilr290 - 1][2]; - double tep, ten, tek, tsp, tsn, tsk; - tqr(cid->u, cid->up, cid->un, cid->tqev, cid->tqsv, tep, ten, tek, tsp, tsn, tsk); - sprintf(virtual_line, " TQEl=%15.7lE, TQEr=%15.7lE, TQEk=%15.7lE\n\0", tep, ten, tek); - output->append_line(virtual_line); - sprintf(virtual_line, " TQSl=%15.7lE, TQSr=%15.7lE, TQSk=%15.7lE\n\0", tsp, tsn, tsk); - output->append_line(virtual_line); - sprintf( - virtual_line, " TQEx=%15.7lE, TQEy=%15.7lE, TQEz=%15.7lE\n\0", - cid->tqce[ilr290 - 1][0], cid->tqce[ilr290 - 1][1], cid->tqce[ilr290 - 1][2] - ); - output->append_line(virtual_line); - sprintf( - virtual_line, " TQSx=%15.7lE, TQSy=%15.7lE, TQSz=%15.7lE\n\0", - cid->tqcs[ilr290 - 1][0], cid->tqcs[ilr290 - 1][1], cid->tqcs[ilr290 - 1][2] - ); - output->append_line(virtual_line); - } - } //ilr290 loop - double rbirif = (real(cid->c1ao->fsac[0][0]) - real(cid->c1ao->fsac[1][1])) / real(cid->c1ao->fsac[0][0]); - double rdichr = (imag(cid->c1ao->fsac[0][0]) - imag(cid->c1ao->fsac[1][1])) / imag(cid->c1ao->fsac[0][0]); - sprintf(virtual_line, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n\0", rbirif); + } + } //ilr290 loop + double rbirif = (real(cid->c1ao->fsac[0][0]) - real(cid->c1ao->fsac[1][1])) / real(cid->c1ao->fsac[0][0]); + double rdichr = (imag(cid->c1ao->fsac[0][0]) - imag(cid->c1ao->fsac[1][1])) / imag(cid->c1ao->fsac[0][0]); + sprintf(virtual_line, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n\0", rbirif); + output->append_line(virtual_line); + sprintf(virtual_line, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n\0", rdichr); + output->append_line(virtual_line); + sprintf(virtual_line, " MULC\n\0"); + output->append_line(virtual_line); + for (int i = 0; i < 4; i++) { + sprintf( + virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n\0", + cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->cmul[i][3] + ); output->append_line(virtual_line); - sprintf(virtual_line, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n\0", rdichr); + } + sprintf(virtual_line, " MULCLR\n\0"); + output->append_line(virtual_line); + for (int i = 0; i < 4; i++) { + sprintf( + virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n\0", + cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3] + ); output->append_line(virtual_line); + } + if (iavm != 0) { + mmulc(cid->c1ao->vintm, cid->cmullr, cid->cmul); + // Some implicit loops writing to binary. + for (int i = 0; i < 16; i++) { + double value; + value = real(cid->c1ao->vintm[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + value = imag(cid->c1ao->vintm[i]); + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + } + for (int i = 0; i < 4; i++) { + for (int j = 0; j < 4; j++) { + double value = cid->cmul[i][j]; + // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); + vtppoanp->append_line(VirtualBinaryLine(value)); + } + } + sprintf(virtual_line, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n\0", iavm); + output->append_line(virtual_line); + if (inpol == 0) { + sprintf(virtual_line, " LIN\n\0"); + output->append_line(virtual_line); + } else { // label 316 + sprintf(virtual_line, " CIRC\n\0"); + output->append_line(virtual_line); + } + // label 318 sprintf(virtual_line, " MULC\n\0"); output->append_line(virtual_line); for (int i = 0; i < 4; i++) { @@ -1414,74 +1463,27 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf ); output->append_line(virtual_line); } - if (iavm != 0) { - mmulc(cid->c1ao->vintm, cid->cmullr, cid->cmul); - // Some implicit loops writing to binary. - for (int i = 0; i < 16; i++) { - double value; - value = real(cid->c1ao->vintm[i]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - value = imag(cid->c1ao->vintm[i]); - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - } - for (int i = 0; i < 4; i++) { - for (int j = 0; j < 4; j++) { - double value = cid->cmul[i][j]; - // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double)); - vtppoanp->append_line(VirtualBinaryLine(value)); - } - } - sprintf(virtual_line, " CLUSTER (ENSEMBLE AVERAGE, MODE%2d)\n\0", iavm); - output->append_line(virtual_line); - if (inpol == 0) { - sprintf(virtual_line, " LIN\n\0"); - output->append_line(virtual_line); - } else { // label 316 - sprintf(virtual_line, " CIRC\n\0"); - output->append_line(virtual_line); - } - // label 318 - sprintf(virtual_line, " MULC\n\0"); - output->append_line(virtual_line); - for (int i = 0; i < 4; i++) { - sprintf( - virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n\0", - cid->cmul[i][0], cid->cmul[i][1], cid->cmul[i][2], cid->cmul[i][3] - ); - output->append_line(virtual_line); - } - sprintf(virtual_line, " MULCLR\n\0"); - output->append_line(virtual_line); - for (int i = 0; i < 4; i++) { - sprintf( - virtual_line, " %15.7lE%15.7lE%15.7lE%15.7lE\n\0", - cid->cmullr[i][0], cid->cmullr[i][1], cid->cmullr[i][2], cid->cmullr[i][3] - ); - output->append_line(virtual_line); - } - } - // label 420, continues jphs loop - if (isam < 1) phs += sa->phsstp; - } // jphs loop, labeled 480 - if (isam <= 1) thsl += sa->thsstp; - } // jths loop, labeled 482 - ph += sa->phstp; - } // jph484 loop - th += sa->thstp; - } // jth486 loop + } + // label 420, continues jphs loop + if (isam < 1) phs += sa->phsstp; + } // jphs loop, labeled 480 + if (isam <= 1) thsl += sa->thsstp; + } // jths loop, labeled 482 + ph += sa->phstp; + } // jph484 loop + th += sa->thstp; + } // jth486 loop #ifdef USE_NVTX - nvtxRangePop(); + nvtxRangePop(); #endif - interval_end = chrono::high_resolution_clock::now(); - elapsed = interval_end - interval_start; - message = "INFO: angle loop for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n"; - logger->log(message); + interval_end = chrono::high_resolution_clock::now(); + elapsed = interval_end - interval_start; + message = "INFO: angle loop for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n"; + logger->log(message); - logger->log("INFO: finished scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n"); + logger->log("INFO: finished scale iteration " + to_string(jxi488) + " of " + to_string(nxi) + ".\n"); - delete logger; + delete logger; - return jer; - } + return jer; +}