Initiate sph work-flow migration to C++

/*! \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;

}

#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."

		);

	}

}