Replace C++ complex<double> with C __complex__ double (e062deee) · Commits · Giacomo Mulas / NP_TMcode

src/cluster/cluster.cpp

+85 −87

Original line number	Original line	Diff line number	Diff line
			/* Distributed under the terms of GPLv3 or later. See COPYING for details. */

	/*! \file cluster.cpp		/*! \file cluster.cpp
	*		*
	* \brief Implementation of the calculation for a cluster of spheres.		* \brief Implementation of the calculation for a cluster of spheres.
	*/		*/
	#include <complex>
	#include <cstdio>		#include <cstdio>
	#include <exception>		#include <exception>
	#include <fstream>		#include <fstream>
	#include <string>		#include <string>

			#ifndef INCLUDE_TYPES_H_
			#include "../include/types.h"
			#endif

	#ifndef INCLUDE_ERRORS_H_		#ifndef INCLUDE_ERRORS_H_
	#include "../include/errors.h"		#include "../include/errors.h"
	#endif		#endif
	@@ -47,10 +52,6 @@
	#include "../include/algebraic.h"		#include "../include/algebraic.h"
	#endif		#endif

	//#ifdef LAPACK_ILP64
	//#define USE_LAPACK
	//#endif

	using namespace std;		using namespace std;

	/*! \brief C++ implementation of CLU		/*! \brief C++ implementation of CLU
	@@ -136,17 +137,16 @@ void cluster(string config_file, string data_file, string output_path) {
	C6 *c6 = new C6(c4->lmtpo);		C6 *c6 = new C6(c4->lmtpo);
	FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");		FILE *output = fopen((output_path + "/c_OCLU").c_str(), "w");
	int jer = 0, lcalc = 0;		int jer = 0, lcalc = 0;
	complex<double> arg(0.0, 0.0), ccsam(0.0, 0.0);		dcomplex arg = 0.0 + 0.0 * I;
			dcomplex ccsam = 0.0 + 0.0 * I;
	int configurations = 0;		int configurations = 0;
	for (int ci = 1; ci <= nsph; ci++) {		for (int ci = 1; ci <= nsph; ci++) {
	if (sconf->iog_vec[ci -1] >= ci) configurations++;		if (sconf->iog_vec[ci -1] >= ci) configurations++;
	}		}
	C2 *c2 = new C2(nsph, configurations, npnt, npntts);		C2 *c2 = new C2(nsph, configurations, npnt, npntts);
	complex<double> *am = new complex<double>[mxndm];		dcomplex *am = new dcomplex[mxndm];
	//complex<double> *tam = new complex<double>[mxndm];
	for (int ai = 0; ai < mxndm; ai++) {		for (int ai = 0; ai < mxndm; ai++) {
	am[ai] = new complex<double>[mxndm]();		am[ai] = new dcomplex[mxndm]();
	//tam[ai] = new complex<double>[mxndm]();
	}		}
	const int ndi = c4->nsph * c4->nlim;		const int ndi = c4->nsph * c4->nlim;
	C9 c9 = new C9(ndi, c4->nlem, 2 ndi, 2 * c4->nlem);		C9 c9 = new C9(ndi, c4->nlem, 2 ndi, 2 * c4->nlem);
	@@ -154,35 +154,35 @@ void cluster(string config_file, string data_file, string output_path) {
	double *tqev = new double[3]();		double *tqev = new double[3]();
	double *tqsv = new double[3]();		double *tqsv = new double[3]();
	double tqse, tqss, tqce, tqcs;		double tqse, tqss, tqce, tqcs;
	complex<double> tqspe, tqsps, tqcpe, tqcps;		dcomplex tqspe, tqsps, tqcpe, tqcps;
	tqse = new double*[2];		tqse = new double*[2];
	tqspe = new complex<double>*[2];		tqspe = new dcomplex*[2];
	tqss = new double*[2];		tqss = new double*[2];
	tqsps = new complex<double>*[2];		tqsps = new dcomplex*[2];
	tqce = new double*[2];		tqce = new double*[2];
	tqcpe = new complex<double>*[2];		tqcpe = new dcomplex*[2];
	tqcs = new double*[2];		tqcs = new double*[2];
	tqcps = new complex<double>*[2];		tqcps = new dcomplex*[2];
	for (int ti = 0; ti < 2; ti++) {		for (int ti = 0; ti < 2; ti++) {
	tqse[ti] = new double[nsph]();		tqse[ti] = new double[nsph]();
	tqspe[ti] = new complex<double>[nsph]();		tqspe[ti] = new dcomplex[nsph]();
	tqss[ti] = new double[nsph]();		tqss[ti] = new double[nsph]();
	tqsps[ti] = new complex<double>[nsph]();		tqsps[ti] = new dcomplex[nsph]();
	tqce[ti] = new double[3]();		tqce[ti] = new double[3]();
	tqcpe[ti] = new complex<double>[3]();		tqcpe[ti] = new dcomplex[3]();
	tqcs[ti] = new double[3]();		tqcs[ti] = new double[3]();
	tqcps[ti] = new complex<double>[3]();		tqcps[ti] = new dcomplex[3]();
	}		}
	double *gapv = new double[3]();		double *gapv = new double[3]();
	complex<double> gapp, gappm;		dcomplex gapp, gappm;
	double gap, gapm;		double gap, gapm;
	gapp = new complex<double>*[3];		gapp = new dcomplex*[3];
	gappm = new complex<double>*[3];		gappm = new dcomplex*[3];
	gap = new double*[3];		gap = new double*[3];
	gapm = new double*[3];		gapm = new double*[3];
	for (int gi = 0; gi < 3; gi++) {		for (int gi = 0; gi < 3; gi++) {
	gapp[gi] = new complex<double>[2]();		gapp[gi] = new dcomplex[2]();
	gappm[gi] = new complex<double>[2]();		gappm[gi] = new dcomplex[2]();
	gap[gi] = new double[2]();		gap[gi] = new double[2]();
	gapm[gi] = new double[2]();		gapm[gi] = new double[2]();
	}		}
	@@ -215,7 +215,7 @@ void cluster(string config_file, string data_file, string output_path) {
	double scan, cfmp, sfmp, cfsp, sfsp;		double scan, cfmp, sfmp, cfsp, sfsp;
	// End of global variables for CLU		// End of global variables for CLU
	fprintf(output, " READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM\n");		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",		fprintf(output, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",
	nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam		nsph, c4->li, c4->le, mxndm, inpol, npnt, npntts, iavm, isam
	);		);
	fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");		fprintf(output, " READ(IR,*)RXX(I),RYY(I),RZZ(I)\n");
	@@ -361,7 +361,6 @@ void cluster(string config_file, string data_file, string output_path) {
	cms(am, c1, c1ao, c4, c6);		cms(am, c1, c1ao, c4, c6);
	lapack_int ndit = 2 * nsph * c4->nlim;		lapack_int ndit = 2 * nsph * c4->nlim;
	invert_matrix(am, ndit, jer, mxndm);		invert_matrix(am, ndit, jer, mxndm);
	//zinvert(am, ndit, jer);
	if (jer != 0) break; // jxi488 loop: goes to memory clean		if (jer != 0) break; // jxi488 loop: goes to memory clean
	ztm(am, c1, c1ao, c4, c6, c9);		ztm(am, c1, c1ao, c4, c6, c9);
	if (sconf->idfc >= 0) {		if (sconf->idfc >= 0) {
	@@ -383,9 +382,8 @@ void cluster(string config_file, string data_file, string output_path) {
	// label 160		// label 160
	double cs0 = 0.25 * vk * vk * vk / acos(0.0);		double cs0 = 0.25 * vk * vk * vk / acos(0.0);
	double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;		double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;
	std::complex<double> s0(0.0, 0.0);		dcomplex s0 = 0.0 + 0.0 * I;
	scr0(vk, exri, c1, c1ao, c3, c4);		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;		double sqk = vk * vk * sconf->exdc;
	aps(zpv, c4->li, nsph, c1, sqk, gaps);		aps(zpv, c4->li, nsph, c1, sqk, gaps);
	rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps);		rabas(inpol, c4->li, nsph, c1, tqse, tqspe, tqss, tqsps);
	@@ -401,19 +399,19 @@ void cluster(string config_file, string data_file, string output_path) {
	if (c1->nshl[i] != 1) {		if (c1->nshl[i] != 1) {
	fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i]);		fprintf(output, " SIZE=%15.7lE\n", c2->vsz[i]);
	} else { // label 162		} 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());		fprintf(output, " SIZE=%15.7lE, REFRACTIVE INDEX=%15.7lE%15.7lE\n", c2->vsz[i], real(c2->vkt[i]), imag(c2->vkt[i]));
	}		}
	// label 164		// label 164
	fprintf(output, " ----- SCS ----- ABS ----- EXS ----- ALBEDS --\n");		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, " %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, " ---- 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, " %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());		fprintf(output, " FSAS=%15.7lE%15.7lE\n", real(c1->fsas[i]), imag(c1->fsas[i]));
	csch = 2.0 * vk * sqsfi / c1->gcsv[i];		csch = 2.0 * vk * sqsfi / c1->gcsv[i];
	s0 = c1->fsas[i] * exri;		s0 = c1->fsas[i] * exri;
	qschu = s0.imag() * csch;		qschu = imag(s0) * csch;
	pschu = s0.real() * csch;		pschu = real(s0) * csch;
	s0mag = abs(s0) * cs0;		s0mag = cabs(s0) * cs0;
	fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);		fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
	double rapr = c1->sexs[i] - gaps[i];		double rapr = c1->sexs[i] - gaps[i];
	double cosav = gaps[i] / c1->sscs[i];		double cosav = gaps[i] / c1->sscs[i];
	@@ -422,12 +420,12 @@ void cluster(string config_file, string data_file, string output_path) {
	fprintf(output, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[1][i], tqss[1][i]);		fprintf(output, " IPO= 2, TQEk=%15.7lE, TQSk=%15.7lE\n", tqse[1][i], tqss[1][i]);
	}		}
	} // i170 loop		} // i170 loop
	fprintf(output, " FSAT=%15.7lE%15.7lE\n", c3->tfsas.real(), c3->tfsas.imag());		fprintf(output, " FSAT=%15.7lE%15.7lE\n", real(c3->tfsas), imag(c3->tfsas));
	csch = 2.0 * vk * sqsfi / c3->gcs;		csch = 2.0 * vk * sqsfi / c3->gcs;
	s0 = c3->tfsas * exri;		s0 = c3->tfsas * exri;
	qschu = s0.imag() * csch;		qschu = imag(s0) * csch;
	pschu = s0.real() * csch;		pschu = real(s0) * csch;
	s0mag = abs(s0) * cs0;		s0mag = cabs(s0) * cs0;
	fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);		fprintf(output, " QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschu, pschu, s0mag);
	tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&vk), sizeof(double));
	pcrsm0(vk, exri, inpol, c1, c1ao, c4);		pcrsm0(vk, exri, inpol, c1, c1ao, c4);
	@@ -537,24 +535,24 @@ void cluster(string config_file, string data_file, string output_path) {
	for (int i = 0; i < 2; i++) {		for (int i = 0; i < 2; i++) {
	double value = c1ao->scscm[i];		double value = c1ao->scscm[i];
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->scscpm[i].real();		value = real(c1ao->scscpm[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->scscpm[i].imag();		value = imag(c1ao->scscpm[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->ecscm[i];		value = c1ao->ecscm[i];
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->ecscpm[i].real();		value = real(c1ao->ecscpm[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->ecscpm[i].imag();		value = imag(c1ao->ecscpm[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	}		}
	for (int i = 0; i < 3; i++) {		for (int i = 0; i < 3; i++) {
	for (int j = 0; j < 2; j++) {		for (int j = 0; j < 2; j++) {
	double value = gapm[i][j];		double value = gapm[i][j];
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = gappm[i][j].real();		value = real(gappm[i][j]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = gappm[i][j].imag();		value = imag(gappm[i][j]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	}		}
	}		}
	@@ -575,12 +573,12 @@ void cluster(string config_file, string data_file, string output_path) {
	double absrm = abssm / c3->acs;		double absrm = abssm / c3->acs;
	double acsecs = c3->acs / c3->ecs;		double acsecs = c3->acs / c3->ecs;
	if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;		if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;
	complex<double> s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;		dcomplex s0m = c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;
	double qschum = s0m.imag() * csch;		double qschum = imag(s0m) * csch;
	double pschum = s0m.real() * csch;		double pschum = real(s0m) * csch;
	double s0magm = abs(s0m) * cs0;		double s0magm = cabs(s0m) * cs0;
	double rfinrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].real() / c3->tfsas.real();		double rfinrm = real(c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / real(c3->tfsas);
	double extcrm = c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag() / c3->tfsas.imag();		double extcrm = imag(c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / imag(c3->tfsas);
	if (inpol == 0) {		if (inpol == 0) {
	fprintf(output, " LIN %2d\n", ipol);		fprintf(output, " LIN %2d\n", ipol);
	} else { // label 206		} else { // label 206
	@@ -595,9 +593,9 @@ void cluster(string config_file, string data_file, string output_path) {
	fprintf(output, " %14.7lE%15.7lE%15.7lE\n", scarm, absrm, extrm);		fprintf(output, " %14.7lE%15.7lE%15.7lE\n", scarm, absrm, extrm);
	fprintf(		fprintf(
	output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",		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(),		ilr210, ilr210, real(c1ao->fsacm[ilr210 - 1][ilr210 - 1]),
	c1ao->fsacm[ilr210 - 1][ilr210 - 1].imag(), jlr, ilr210,		imag(c1ao->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,
	c1ao->fsacm[jlr - 1][ilr210 - 1].real(), c1ao->fsacm[jlr - 1][ilr210 - 1].imag()		real(c1ao->fsacm[jlr - 1][ilr210 - 1]), imag(c1ao->fsacm[jlr - 1][ilr210 - 1])
	);		);
	fprintf(		fprintf(
	output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",		output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
	@@ -610,8 +608,8 @@ void cluster(string config_file, string data_file, string output_path) {
	fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);		fprintf(output, " COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
	fprintf(output, " Fk=%15.7lE\n", fz);		fprintf(output, " Fk=%15.7lE\n", fz);
	} // ilr210 loop		} // ilr210 loop
	double rmbrif = (c1ao->fsacm[0][0].real() - c1ao->fsacm[1][1].real()) / c1ao->fsacm[0][0].real();		double rmbrif = (real(c1ao->fsacm[0][0]) - real(c1ao->fsacm[1][1])) / real(c1ao->fsacm[0][0]);
	double rmdchr = (c1ao->fsacm[0][0].imag() - c1ao->fsacm[1][1].imag()) / c1ao->fsacm[0][0].imag();		double rmdchr = (imag(c1ao->fsacm[0][0]) - imag(c1ao->fsacm[1][1])) / imag(c1ao->fsacm[0][0]);
	fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rmbrif);		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);		fprintf(output, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rmdchr);
	}		}
	@@ -641,13 +639,13 @@ void cluster(string config_file, string data_file, string output_path) {
	fprintf(output, " SPHERE %2d\n", i226);		fprintf(output, " SPHERE %2d\n", i226);
	fprintf(		fprintf(
	output, " SAS(1,1)=%15.7lE%15.7lE, SAS(2,1)=%15.7lE%15.7lE\n",		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(),		real(c1->sas[i226 - 1][0][0]), imag(c1->sas[i226 - 1][0][0]),
	c1->sas[i226 - 1][1][0].real(), c1->sas[i226 - 1][1][0].imag()		real(c1->sas[i226 - 1][1][0]), imag(c1->sas[i226 - 1][1][0])
	);		);
	fprintf(		fprintf(
	output, " SAS(1,2)=%15.7lE%15.7lE, SAS(2,2)=%15.7lE%15.7lE\n",		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(),		real(c1->sas[i226 - 1][0][1]), imag(c1->sas[i226 - 1][0][1]),
	c1->sas[i226 - 1][1][1].real(), c1->sas[i226 - 1][1][1].imag()		real(c1->sas[i226 - 1][1][1]), imag(c1->sas[i226 - 1][1][1])
	);		);
	for (int j225 = 0; j225 < 16; j225++) { // QUESTION: check that 16 is a fixed dimension		for (int j225 = 0; j225 < 16; j225++) { // QUESTION: check that 16 is a fixed dimension
	c1ao->vint[j225] = c1ao->vints[i226 - 1][j225];		c1ao->vint[j225] = c1ao->vints[i226 - 1][j225];
	@@ -671,13 +669,13 @@ void cluster(string config_file, string data_file, string output_path) {
	} // i226 loop		} // i226 loop
	fprintf(		fprintf(
	output, " SAT(1,1)=%15.7lE%15.7lE, SAT(2,1)=%15.7lE%15.7lE\n",		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(),		real(c3->tsas[0][0]), imag(c3->tsas[0][0]),
	c3->tsas[1][0].real(), c3->tsas[1][0].imag()		real(c3->tsas[1][0]), imag(c3->tsas[1][0])
	);		);
	fprintf(		fprintf(
	output, " SAT(1,2)=%15.7lE%15.7lE, SAT(2,2)=%15.7lE%15.7lE\n",		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(),		real(c3->tsas[0][1]), imag(c3->tsas[0][1]),
	c3->tsas[1][1].real(), c3->tsas[1][1].imag()		real(c3->tsas[1][1]), imag(c3->tsas[1][1])
	);		);
	fprintf(output, " CLUSTER\n");		fprintf(output, " CLUSTER\n");
	pcros(vk, exri, c1, c1ao, c4);		pcros(vk, exri, c1, c1ao, c4);
	@@ -695,24 +693,24 @@ void cluster(string config_file, string data_file, string output_path) {
	for (int i = 0; i < 2; i++) {		for (int i = 0; i < 2; i++) {
	double value = c1ao->scsc[i];		double value = c1ao->scsc[i];
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->scscp[i].real();		value = real(c1ao->scscp[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->scscp[i].imag();		value = imag(c1ao->scscp[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->ecsc[i];		value = c1ao->ecsc[i];
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->ecscp[i].real();		value = real(c1ao->ecscp[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->ecscp[i].imag();		value = imag(c1ao->ecscp[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	}		}
	for (int i = 0; i < 3; i++) {		for (int i = 0; i < 3; i++) {
	for (int j = 0; j < 2; j++) {		for (int j = 0; j < 2; j++) {
	double value = gap[i][j];		double value = gap[i][j];
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = gapp[i][j].real();		value = real(gapp[i][j]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = gapp[i][j].imag();		value = imag(gapp[i][j]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	}		}
	}		}
	@@ -720,9 +718,9 @@ void cluster(string config_file, string data_file, string output_path) {
	for (int j = 0; j < 3; j++) {		for (int j = 0; j < 3; j++) {
	double value = tqce[i][j];		double value = tqce[i][j];
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = tqcpe[i][j].real();		value = real(tqcpe[i][j]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = tqcpe[i][j].imag();		value = imag(tqcpe[i][j]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	}		}
	}		}
	@@ -730,9 +728,9 @@ void cluster(string config_file, string data_file, string output_path) {
	for (int j = 0; j < 3; j++) {		for (int j = 0; j < 3; j++) {
	double value = tqcs[i][j];		double value = tqcs[i][j];
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = tqcps[i][j].real();		value = real(tqcps[i][j]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = tqcps[i][j].imag();		value = imag(tqcps[i][j]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	}		}
	}		}
	@@ -747,9 +745,9 @@ void cluster(string config_file, string data_file, string output_path) {
	}		}
	// label 254		// label 254
	for (int i = 0; i < 16; i++) {		for (int i = 0; i < 16; i++) {
	double value = c1ao->vint[i].real();		double value = real(c1ao->vint[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->vint[i].imag();		value = imag(c1ao->vint[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	}		}
	for (int i = 0; i < 4; i++) {		for (int i = 0; i < 4; i++) {
	@@ -775,11 +773,11 @@ void cluster(string config_file, string data_file, string output_path) {
	double ratio = c3->acs / c3->ecs;		double ratio = c3->acs / c3->ecs;
	if (ratio < -1.0e-6 \|\| ratio > 1.0e-6) absrat = abssec / c3->acs;		if (ratio < -1.0e-6 \|\| ratio > 1.0e-6) absrat = abssec / c3->acs;
	s0 = c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;		s0 = c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;
	double qschu = s0.imag() * csch;		double qschu = imag(s0) * csch;
	double pschu = s0.real() * csch;		double pschu = real(s0) * csch;
	s0mag = abs(s0) * cs0;		s0mag = cabs(s0) * cs0;
	double refinr = c1ao->fsac[ilr290 - 1][ilr290 - 1].real() / c3->tfsas.real();		double refinr = real(c1ao->fsac[ilr290 - 1][ilr290 - 1]) / real(c3->tfsas);
	double extcor = c1ao->fsac[ilr290 - 1][ilr290 - 1].imag() / c3->tfsas.imag();		double extcor = imag(c1ao->fsac[ilr290 - 1][ilr290 - 1]) / imag(c3->tfsas);
	if (inpol == 0) {		if (inpol == 0) {
	fprintf(output, " LIN %2d\n", ipol);		fprintf(output, " LIN %2d\n", ipol);
	} else { // label 273		} else { // label 273
	@@ -803,13 +801,13 @@ void cluster(string config_file, string data_file, string output_path) {
	);		);
	fprintf(		fprintf(
	output, " FSAC(%1d,%1d)=%15.7lE%15.7lE FSAC(%1d,%1d)=%15.7lE%15.7lE\n",		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(),		ilr290, ilr290, real(c1ao->fsac[ilr290 - 1][ilr290 - 1]), imag(c1ao->fsac[ilr290 - 1][ilr290 - 1]),
	jlr, ilr290, c1ao->fsac[jlr - 1][ilr290 - 1].real(), c1ao->fsac[jlr - 1][ilr290 - 1].imag()		jlr, ilr290, real(c1ao->fsac[jlr - 1][ilr290 - 1]), imag(c1ao->fsac[jlr - 1][ilr290 - 1])
	);		);
	fprintf(		fprintf(
	output, " SAC(%1d,%1d)=%15.7lE%15.7lE SAC(%1d,%1d)=%15.7lE%15.7lE\n",		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(),		ilr290, ilr290, real(c1ao->sac[ilr290 - 1][ilr290 - 1]), imag(c1ao->sac[ilr290 - 1][ilr290 - 1]),
	jlr, ilr290, c1ao->sac[jlr - 1][ilr290 - 1].real(), c1ao->sac[jlr - 1][ilr290 - 1].imag()		jlr, ilr290, real(c1ao->sac[jlr - 1][ilr290 - 1]), imag(c1ao->sac[jlr - 1][ilr290 - 1])
	);		);
	fprintf(		fprintf(
	output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",		output, " RE(FSAC(%1d,%1d))/RE(TFSAS)=%15.7lE, IM(FSAC(%1d,%1d))/IM(TFSAS)=%15.7lE\n",
	@@ -851,8 +849,8 @@ void cluster(string config_file, string data_file, string output_path) {
	);		);
	}		}
	} //ilr290 loop		} //ilr290 loop
	double rbirif = (c1ao->fsac[0][0].real() - c1ao->fsac[1][1].real()) / c1ao->fsac[0][0].real();		double rbirif = (real(c1ao->fsac[0][0]) - real(c1ao->fsac[1][1])) / real(c1ao->fsac[0][0]);
	double rdichr = (c1ao->fsac[0][0].imag() - c1ao->fsac[1][1].imag()) / c1ao->fsac[0][0].imag();		double rdichr = (imag(c1ao->fsac[0][0]) - imag(c1ao->fsac[1][1])) / imag(c1ao->fsac[0][0]);
	fprintf(output, " (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rbirif);		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, " (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rdichr);
	fprintf(output, " MULC\n");		fprintf(output, " MULC\n");
	@@ -874,9 +872,9 @@ void cluster(string config_file, string data_file, string output_path) {
	// Some implicit loops writing to binary.		// Some implicit loops writing to binary.
	for (int i = 0; i < 16; i++) {		for (int i = 0; i < 16; i++) {
	double value;		double value;
	value = c1ao->vintm[i].real();		value = real(c1ao->vintm[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	value = c1ao->vintm[i].imag();		value = imag(c1ao->vintm[i]);
	tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));		tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));
	}		}
	for (int i = 0; i < 4; i++) {		for (int i = 0; i < 4; i++) {

src/cluster/np_cluster.cpp

+6 −1

Original line number	Original line	Diff line number	Diff line
			/* Distributed under the terms of GPLv3 or later. See COPYING for details. */

	/*! \file np_cluster.cpp		/*! \file np_cluster.cpp
	*		*
	* \brief Cluster of spheres scattering problem handler.		* \brief Cluster of spheres scattering problem handler.
	@@ -17,9 +19,12 @@
	*/		*/

	#include <cstdio>		#include <cstdio>
	#include <complex>
	#include <string>		#include <string>

			#ifndef INCLUDE_TYPES_H_
			#include "../include/types.h"
			#endif

	#ifndef INCLUDE_CONFIGURATION_H_		#ifndef INCLUDE_CONFIGURATION_H_
	#include "../include/Configuration.h"		#include "../include/Configuration.h"
	#endif		#endif

src/include/Configuration.h

+5 −3

Original line number	Original line	Diff line number	Diff line
			/* Distributed under the terms of GPLv3 or later. See COPYING for details. */

	/*! \file Configuration.h		/*! \file Configuration.h
	*		*
	* \brief Configuration data structures.		* \brief Configuration data structures.
	@@ -168,7 +170,7 @@ class ScattererConfiguration {
	friend void sphere(std::string, std::string, std::string);		friend void sphere(std::string, std::string, std::string);
	protected:		protected:
	//! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].		//! \brief Matrix of dielectric parameters with size [NON_TRANS_LAYERS x N_SPHERES x N_SCALES].
	std::complex<double> ***dc0_matrix;		dcomplex ***dc0_matrix;
	//! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].		//! \brief Vector of sphere radii expressed in m, with size [N_SPHERES].
	double *radii_of_spheres;		double *radii_of_spheres;
	//! \brief Matrix of fractional transition radii with size [N_SPHERES x LAYERS].		//! \brief Matrix of fractional transition radii with size [N_SPHERES x LAYERS].
	@@ -258,7 +260,7 @@ public:
	* \param dielectric_func_type: `int` Type of dielectric function definition (=0 for constant,		* \param dielectric_func_type: `int` Type of dielectric function definition (=0 for constant,
	* \>0 as function of scale parameter, <0 for functions at XIP value and XI is scale factor		* \>0 as function of scale parameter, <0 for functions at XIP value and XI is scale factor
	* for dimensions).		* for dimensions).
	* \param dc_matrix: Matrix of reference dielectric constants.		* \param dc_matrix: `complex double ***` Matrix of reference dielectric constants.
	* \param has_external: `bool` Flag to set whether to add an external spherical layer.		* \param has_external: `bool` Flag to set whether to add an external spherical layer.
	* \param exdc: `double` External medium dielectric constant.		* \param exdc: `double` External medium dielectric constant.
	* \param wp: `double` wp		* \param wp: `double` wp
	@@ -274,7 +276,7 @@ public:
	int *nshl_vector,		int *nshl_vector,
	double **rcf_vector,		double **rcf_vector,
	int dielectric_func_type,		int dielectric_func_type,
	std::complex<double> ***dc_matrix,		dcomplex ***dc_matrix,
	bool has_external,		bool has_external,
	double exdc,		double exdc,
	double wp,		double wp,

src/include/clu_subs.h

+31 −32

Original line number	Original line	Diff line number	Diff line
			/* Distributed under the terms of GPLv3 or later. See COPYING for details. */

	/*! \file clu_subs.h		/*! \file clu_subs.h
	*		*
	* \brief C++ porting of CLU functions and subroutines.		* \brief C++ porting of CLU functions and subroutines.
	@@ -31,15 +33,15 @@
	*		*
	* \param zpv: `double ****`		* \param zpv: `double ****`
	* \param le: `int`		* \param le: `int`
	* \param am0m: Matrix of complex.		* \param am0m: `complex double **`
	* \param w: Matrix of complex.		* \param w: `complex double **`
	* \param sqk: `double`		* \param sqk: `double`
	* \param gapr: `double **`		* \param gapr: `double **`
	* \param gapp: Matrix of complex.		* \param gapp: `complex double **`
	*/		*/
	void apc(		void apc(
	double **zpv, int le, std::complex<double> am0m, std::complex<double> **w,		double **zpv, int le, dcomplex am0m, dcomplex **w,
	double sqk, double gapr, std::complex<double> gapp		double sqk, double gapr, dcomplex gapp
	);		);

	/*! \brief Compute the asymmetry-corrected scattering cross-section under random average		/*! \brief Compute the asymmetry-corrected scattering cross-section under random average
	@@ -50,32 +52,29 @@ void apc(
	*		*
	* \param zpv: `double ****`		* \param zpv: `double ****`
	* \param le: `int`		* \param le: `int`
	* \param am0m: Matrix of complex.		* \param am0m: `complex double **`
	* \param inpol: `int` Polarization type.		* \param inpol: `int` Polarization type.
	* \param sqk: `double`		* \param sqk: `double`
	* \param gaprm: `double **`		* \param gaprm: `double **`
	* \param gappm: Matrix of complex.		* \param gappm: `complex double **`
	*/		*/
	void apcra(		void apcra(
	double **zpv, const int le, std::complex<double> am0m, int inpol, double sqk,		double **zpv, const int le, dcomplex am0m, int inpol, double sqk,
	double gaprm, std::complex<double> gappm		double gaprm, dcomplex gappm
	);		);

	/*! \brief Complex inner product.		/*! \brief Complex inner product.
	*		*
	* This function performs the complex inner product. It is used by `lucin()`.		* This function performs the complex inner product. It is used by `lucin()`.
	*		*
	* \param z: `complex<double>`		* \param z: `complex double`
	* \param am: Matrix of complex.		* \param am: `complex double **`
	* \param i: `int`		* \param i: `int`
	* \param jf: `int`		* \param jf: `int`
	* \param k: `int`		* \param k: `int`
	* \param nj: `int`		* \param nj: `int`
	*/		*/
	std::complex<double> cdtp(		dcomplex cdtp(dcomplex z, dcomplex **am, int i, int jf, int k, int nj);
	std::complex<double> z, std::complex<double> **am, int i, int jf,
	int k, int nj
	);

	/*! \brief C++ porting of CGEV. QUESTION: description?		/*! \brief C++ porting of CGEV. QUESTION: description?
	*		*
	@@ -92,13 +91,13 @@ double cgev(int ipamo, int mu, int l, int m);
	* This function constructs the multi-centered M-matrix of the cluster, according		* This function constructs the multi-centered M-matrix of the cluster, according
	* to Eq. (5.28) of Borghese, Denti & Saija (2007).		* to Eq. (5.28) of Borghese, Denti & Saija (2007).
	*		*
	* \param am: Matrix of complex.		* \param am: `complex double **`
	* \param c1: `C1 *`		* \param c1: `C1 *`
	* \param c1ao: `C1_AddOns *`		* \param c1ao: `C1_AddOns *`
	* \param c4: `C4 *`		* \param c4: `C4 *`
	* \param c6: `C6 *`		* \param c6: `C6 *`
	*/		*/
	void cms(std::complex<double> *am, C1 c1, C1_AddOns c1ao, C4 c4, C6 *c6);		void cms(dcomplex *am, C1 c1, C1_AddOns c1ao, C4 c4, C6 *c6);

	/*! \brief Compute orientation-averaged scattered field intensity.		/*! \brief Compute orientation-averaged scattered field intensity.
	*		*
	@@ -132,7 +131,7 @@ void crsm1(double vk, double exri, C1 c1, C1_AddOns c1ao, C4 c4, C6 c6);
	* \param c4: `C4 *`		* \param c4: `C4 *`
	* \param c6: `C6 *`		* \param c6: `C6 *`
	*/		*/
	std::complex<double> ghit(		dcomplex ghit(
	int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,		int ihi, int ipamo, int nbl, int l1, int m1, int l2, int m2, C1 *c1,
	C1_AddOns c1ao, C4 c4, C6 *c6		C1_AddOns c1ao, C4 c4, C6 *c6
	);		);
	@@ -152,7 +151,7 @@ std::complex<double> ghit(
	* \param c4: `C4 *`		* \param c4: `C4 *`
	*/		*/
	void hjv(		void hjv(
	double exri, double vk, int &jer, int &lcalc, std::complex<double> &arg,		double exri, double vk, int &jer, int &lcalc, dcomplex &arg,
	C1 c1, C1_AddOns c1ao, C4 *c4		C1 c1, C1_AddOns c1ao, C4 *c4
	);		);

	@@ -161,12 +160,12 @@ void hjv(
	* This function performs the inversion of the multi-centered M-matrix through		* This function performs the inversion of the multi-centered M-matrix through
	* LU decomposition. See Eq. (5.29) in Borghese, Denti & Saija (2007).		* LU decomposition. See Eq. (5.29) in Borghese, Denti & Saija (2007).
	*		*
	* \param am: Matrix of complex.		* \param am: `complex double **`
	* \param nddmst: `const int64_t`		* \param nddmst: `const int64_t`
	* \param n: `int64_t`		* \param n: `int64_t`
	* \param ier: `int &`		* \param ier: `int &`
	*/		*/
	void lucin(std::complex<double> **am, const np_int nddmst, np_int n, int &ier);		void lucin(dcomplex **am, const np_int nddmst, np_int n, int &ier);

	/*! \brief Compute the average extinction cross-section.		/*! \brief Compute the average extinction cross-section.
	*		*
	@@ -180,7 +179,7 @@ void lucin(std::complex<double> **am, const np_int nddmst, np_int n, int &ier);
	* \param cextlr: `double **`		* \param cextlr: `double **`
	* \param cext: `double **`		* \param cext: `double **`
	*/		*/
	void mextc(double vk, double exri, std::complex<double> fsac, double cextlr, double **cext);		void mextc(double vk, double exri, dcomplex fsac, double cextlr, double **cext);

	/*! \brief Compute cross-sections and forward scattering amplitude for the cluster.		/*! \brief Compute cross-sections and forward scattering amplitude for the cluster.
	*		*
	@@ -274,16 +273,16 @@ void r3jmr(int j1, int j2, int j3, int m1, C6 *c6);
	* in Borghese, Denti & Saija (2007).		* in Borghese, Denti & Saija (2007).
	*		*
	* \param le: `int`		* \param le: `int`
	* \param am0m: Matrix of complex.		* \param am0m: `complex double **`
	* \param w: Matrix of complex.		* \param w: `complex double **`
	* \param tqce: `double **`		* \param tqce: `double **`
	* \param tqcpe: Matrix of complex.		* \param tqcpe: `complex double **`
	* \param tqcs: `double **`		* \param tqcs: `double **`
	* \param tqcps: Matrix of complex.		* \param tqcps: `complex double **`
	*/		*/
	void raba(		void raba(
	int le, std::complex<double> am0m, std::complex<double> w, double **tqce,		int le, dcomplex am0m, dcomplex w, double **tqce,
	std::complex<double> tqcpe, double tqcs, std::complex<double> **tqcps		dcomplex tqcpe, double tqcs, dcomplex **tqcps
	);		);

	/*! \brief Compute the radiation force Cartesian components.		/*! \brief Compute the radiation force Cartesian components.
	@@ -391,13 +390,13 @@ void tqr(
	* This function computes the single-centered inverrted M-matrix appearing in Eq. (5.28)		* This function computes the single-centered inverrted M-matrix appearing in Eq. (5.28)
	* of Borghese, Denti & Saija (2007).		* of Borghese, Denti & Saija (2007).
	*		*
	* \param am: Matrix of complex.		* \param am: `complex double **`
	* \param c1: `C1 *` Pointer to a C1 instance.		* \param c1: `C1 *` Pointer to a C1 instance.
	* \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.		* \param c1ao: `C1_AddOns *` Pointer to C1_AddOns instance.
	* \param c4: `C4 *` Pointer to a C4 structure.		* \param c4: `C4 *` Pointer to a C4 structure.
	* \param c6: `C6 *` Pointer to a C6 instance.		* \param c6: `C6 *` Pointer to a C6 instance.
	* \param c9: `C9 *` Pointer to a C9 instance.		* \param c9: `C9 *` Pointer to a C9 instance.
	*/		*/
	void ztm(std::complex<double> *am, C1 c1, C1_AddOns c1ao, C4 c4, C6 c6, C9 c9);		void ztm(dcomplex *am, C1 c1, C1_AddOns c1ao, C4 c4, C6 c6, C9 c9);

	#endif		#endif

src/include/tfrfme.h

+1 −1

Original line number	Original line	Diff line number	Diff line
	@@ -87,7 +87,7 @@ public:
	*		*
	* \return value: `complex double *` Memory address of the WK vector.		* \return value: `complex double *` Memory address of the WK vector.
	*/		*/
	std::complex<double> *get_vector() { return wk; }		dcomplex *get_vector() { return wk; }

	/*! \brief Bring the pointer to the next element at the start of vector.		/*! \brief Bring the pointer to the next element at the start of vector.
	*/		*/