From 50a311fab0b405928ed694f504ccc442f7eba7ba Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Thu, 26 Oct 2023 15:24:27 +0200
Subject: [PATCH] Initiate sph work-flow migration to C++
---
src/nptm.cpp | 23 ++++
src/sphere.cpp | 279 +++++++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 302 insertions(+)
create mode 100644 src/nptm.cpp
create mode 100644 src/sphere.cpp
diff --git a/src/nptm.cpp b/src/nptm.cpp
new file mode 100644
index 00000000..0d74bb33
--- /dev/null
+++ b/src/nptm.cpp
@@ -0,0 +1,23 @@
+/*! \file nptm.cpp
+ */
+
+#include <cstdio>
+#include <string>
+#include "include/Configuration.h"
+
+using namespace std;
+
+extern void sphere();
+
+/*! \brief Main program entry point.
+ *
+ * This is the starting point of the execution flow. Here we may choose
+ * how to configure the code, e.g. by loading a legacy configuration file
+ * or some otherwise formatted configuration data set. The code can be
+ * linked to a luncher script or to a GUI oriented application that performs
+ * the configuration and runs the main program.
+ */
+int main(int argc, char **argv) {
+ sphere();
+ return 0;
+}
diff --git a/src/sphere.cpp b/src/sphere.cpp
new file mode 100644
index 00000000..1220ddbf
--- /dev/null
+++ b/src/sphere.cpp
@@ -0,0 +1,279 @@
+#include <cstdio>
+#include <fstream>
+#include <string>
+#include <complex>
+#include "include/Configuration.h"
+#include "include/sph_subs.h"
+
+using namespace std;
+
+//! \brief C++ implementation of SPH
+void sphere() {
+ //complex<double> dc0[5];
+ //complex<double> dc0m[2][4];
+ complex<double> arg, s0, tfsas;
+ printf("INFO: making legacy configuration ...\n");
+ ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/sphere/DEDFB");
+ conf->write_formatted("c_OEDFB");
+ conf->write_binary("c_TEDF");
+ delete conf;
+ printf("INFO: reading binary configuration ...\n");
+ ScattererConfiguration *sconf = ScattererConfiguration::from_binary("c_TEDF");
+ GeometryConfiguration *gconf = GeometryConfiguration::from_legacy("../../test_data/sphere/DSPH");
+ if (sconf->number_of_spheres == gconf->number_of_spheres) {
+ int nsph = gconf->number_of_spheres;
+ C1 *c1 = new C1;
+ for (int i = 0; i < nsph; i++) {
+ c1->iog[i] = sconf->iog_vec[i];
+ c1->nshl[i] = sconf->nshl_vec[i];
+ c1->ros[i] = sconf->radii_of_spheres[i];
+ }
+ for (int i = 0; i < gconf->l_max; i++) {
+ c1->rmi[i][0] = complex<double>(0.0, 0.0);
+ c1->rmi[i][1] = complex<double>(0.0, 0.0);
+ c1->rei[i][0] = complex<double>(0.0, 0.0);
+ c1->rei[i][1] = complex<double>(0.0, 0.0);
+ }
+ C2 *c2 = new C2;
+ FILE *output = fopen("c_OSPH", "w");
+ fprintf(output, " READ(IR,*)NSPH,LM,INPOL,NPNT,NPNTTS,ISAM\n");
+ fprintf(
+ output,
+ " %5d%5d%5d%5d%5d%5d\n",
+ gconf->number_of_spheres,
+ gconf->l_max,
+ gconf->in_pol,
+ gconf->npnt,
+ gconf->npntts,
+ gconf->meridional_type
+ );
+ fprintf(output, " READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST\n");
+ fprintf(
+ output,
+ " %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",
+ gconf->in_theta_start,
+ gconf->in_theta_step,
+ gconf->in_theta_end,
+ gconf->sc_theta_start,
+ gconf->sc_theta_step,
+ gconf->sc_theta_end
+ );
+ fprintf(output, " READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST\n");
+ fprintf(
+ output,
+ " %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE %9.3lE\n",
+ gconf->in_phi_start,
+ gconf->in_phi_step,
+ gconf->in_phi_end,
+ gconf->sc_phi_start,
+ gconf->sc_phi_step,
+ gconf->sc_phi_end
+ );
+ 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 \n");
+ double sml = 1.0e-3;
+ int nth = 0, nph = 0;
+ if (gconf->in_theta_step != 0.0)
+ nth = int((gconf->in_theta_end - gconf->in_theta_start) / gconf->in_theta_step + sml);
+ nth += 1;
+ if (gconf->in_phi_step != 0.0)
+ nph = int((gconf->in_phi_end - gconf->in_phi_start) / gconf->in_phi_step + sml);
+ nph += 1;
+ int nths = 0, nphs = 0;
+ double thsca;
+ if (gconf->meridional_type > 1) { // ISAM > 1, fixed scattering angle
+ nths = 1;
+ thsca = gconf->sc_theta_start - gconf->in_theta_start;
+ } else { //ISAM <= 1
+ if (gconf->in_theta_step != 0.0)
+ nths = int((gconf->sc_theta_end - gconf->sc_theta_start) / gconf->sc_theta_step + sml);
+ nths += 1;
+ if (gconf->meridional_type == 1) { // ISAM = 1
+ nphs = 1;
+ } else { // ISAM < 1
+ if (gconf->sc_phi_step != 0.0)
+ nphs = int((gconf->sc_phi_end - gconf->sc_phi_start) / gconf->sc_phi_step + sml);
+ nphs += 1;
+ }
+ }
+ int nk = nth * nph;
+ int nks = nths * nphs;
+ int nkks = nk * nks;
+ int th1 = gconf->in_theta_start;
+ int ph1 = gconf->in_phi_start;
+ int ths1 = gconf->sc_theta_start;
+ int phs1 = gconf->sc_phi_start;
+ const double half_pi = acos(0.0);
+ const double pi = 2.0 * half_pi;
+ double gcs = 0.0;
+ for (int i116 = 1; i116 <= nsph; i116++) {
+ int iogi = c1->iog[i116 - 1];
+ if (iogi >= i116) {
+ double gcss = pi * c1->ros[i116 - 1] * c1->ros[i116 - 1];
+ c1->gcsv[i116 - 1] = gcss;
+ int nsh = c1->nshl[i116 - 1];
+ for (int j115 = 1; j115 <= nsh; j115++) {
+ c1->rc[i116 - 1][j115 - 1] = sconf->rcf[i116 - 1][j115 - 1] * c1->ros[i116 - 1];
+ }
+ }
+ gcs += c1->gcsv[iogi - 1];
+ }
+ double ****zpv = new double***[gconf->l_max]; //[gconf->l_max][3][2][2]; // Matrix: dim[LM x 3 x 2 x 2]
+ for (int zi = 0; zi < gconf->l_max; 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(gconf->l_max, zpv);
+ //DEBUG CODE
+ /*
+ for (int zi = 0; zi < gconf->l_max; zi++) {
+ for (int zj = 0; zj < 3; zj++) {
+ for (int zk = 0; zk < 2; zk++) {
+ if (zpv[zi][zj][zk][0] != 0.0)
+ printf(
+ "DEBUG: zpv[%d][%d][%d][0] = %.3lE\n",
+ zi, zj, zk, zpv[zi][zj][zk][0]
+ );
+ if (zpv[zi][zj][zk][1] != 0.0)
+ printf(
+ "DEBUG: zpv[%d][%d][%d][1] = %.3lE\n",
+ zi, zj, zk, zpv[zi][zj][zk][1]
+ );
+ }
+ }
+ }
+ */
+ //END OF DEBUG CODE
+ double exri = sqrt(sconf->exdc);
+ fprintf(output, " REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
+ fstream tppoan;
+ tppoan.open("c_TPPOAN", ios::binary|ios::out);
+ if (tppoan.is_open()) {
+ int imode = 10;
+ tppoan.write(reinterpret_cast<char *>(&imode), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&(gconf->meridional_type)), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&(gconf->in_pol)), sizeof(int));
+ tppoan.write(reinterpret_cast<char *>(&(sconf->number_of_scales)), 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));
+
+ for (int nsi = 0; nsi < nsph; nsi++)
+ tppoan.write(reinterpret_cast<char *>(&(sconf->iog_vec[nsi])), sizeof(int));
+ if (gconf->in_pol == 0) fprintf(output, " LIN\n");
+ else fprintf(output, " CIRC\n");
+ fprintf(output, " \n");
+ double wn = sconf->wp / 3.0e8;
+ double sqsfi = 1.0;
+ double vk, vkarg;
+ 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 <= sconf->number_of_scales; jxi488++) {
+ fprintf(output, "========== JXI =%3d ====================\n", jxi488);
+ double xi = sconf->scale_vec[jxi488 - 1];
+ if (sconf->idfc >= 0) {
+ vk = xi * wn;
+ vkarg = vk;
+ fprintf(output, " VK=%15.7lE, XI=%15.7lE\n", xi, vk);
+ } else { // IDFC < 0
+ vkarg = xi * vk;
+ sqsfi = 1.0 / (xi * xi);
+ fprintf(output, " XI=%15.7lE\n", xi);
+ }
+ tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
+ for (int i132 = 1; i132 <= nsph; i132++) {
+ int iogi = sconf->iog_vec[i132 - 1];
+ if (iogi != i132) {
+ for (int l123 = 1; l123 <= gconf->l_max; l123++) {
+ // QUESTION: aren't RMI and REI still empty?
+ c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];
+ c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];
+ }
+ continue; // i132
+ }
+ 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][0]; // WARNING: IDFC=0 is not tested!
+ } else { // IDFC != 0
+ if (jxi488 == 1) {
+ for (int ic = 0; ic < ici; ic++) {
+ c2->dc0[ic] = sconf->dc0_matrix[ic][i132 - 1][jxi488 - 1];
+ }
+ }
+ }
+ if (nsh % 2 == 0) c2->dc0[ici] = sconf->exdc;
+ int jer = 0;
+ int lcalc = 0;
+ dme(
+ gconf->l_max, i132, gconf->npnt, gconf->npntts, vkarg,
+ sconf->exdc, exri, c1, c2, jer, lcalc, arg
+ );
+ if (jer != 0) {
+ fprintf(output, " STOP IN DME\n");
+ fprintf(output, " AT %1d LCALC=%3d TOO SMALL WITH ARG=%15.7lE+i(%15.7lE)\n", jer, lcalc, arg.real(), arg.imag());
+ tppoan.close();
+ fclose(output);
+ return;
+ }
+ } // i132
+ if (sconf->idfc >= 0 and nsph == 1 and jxi488 == gconf->jwtm) {
+ // This is the condition that writes the transition matrix to output.
+ int is = 1111;
+ fstream ttms;
+ ttms.open("c_TTMS", ios::binary | ios::out);
+ if (ttms.is_open()) {
+ ttms.write(reinterpret_cast<char *>(&is), sizeof(int));
+ ttms.write(reinterpret_cast<char *>(&(gconf->l_max)), sizeof(int));
+ ttms.write(reinterpret_cast<char *>(&vk), sizeof(double));
+ ttms.write(reinterpret_cast<char *>(&exri), sizeof(double));
+ for (int lmi = 0; lmi < gconf->l_max; lmi++) {
+ complex<double> element1 = -1.0 / c1->rmi[0][lmi];
+ complex<double> element2 = -1.0 / c1->rei[0][lmi];
+ ttms.write(reinterpret_cast<char *>(&element1), sizeof(complex<double>));
+ ttms.write(reinterpret_cast<char *>(&element2), sizeof(complex<double>));
+ }
+ ttms.write(reinterpret_cast<char *>(&(sconf->radii_of_spheres[0])), sizeof(double));
+ ttms.close();
+ } else { // Failed to open output file. Should never happen.
+ printf("ERROR: could not open TTMS file.\n");
+ tppoan.close();
+ fclose(output);
+ }
+ }
+ // We are at line 271 of SPH.f
+ double cs0 = 0.25 * vk * vk * vk / half_pi;
+ sscr0(tfsas, nsph, gconf->l_max, vk, exri, c1);
+ printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag());
+ } //jxi488
+ tppoan.close();
+ } else { // In case TPPOAN could not be opened. Should never happen.
+ printf("ERROR: failed to open TPPOAN file.\n");
+ }
+ fclose(output);
+ } else { // NSPH mismatch between geometry and scatterer configurations.
+ throw UnrecognizedConfigurationException(
+ "Inconsistent geometry and scatterer configurations."
+ );
+ }
+}
--
GitLab