Solve merge conflicts for OMP implementation (c04023ad) · Commits · Giacomo Mulas / NP_TMcode

src/cluster/cluster.cpp

+1131 −625

File changed.

Preview size limit exceeded, changes collapsed.

src/include/Commons.h

+78 −3

Original line number	Original line	Diff line number	Diff line
	@@ -40,8 +40,6 @@ protected:
	int nsph;		int nsph;
	//! \brief Maximum order of field expansion.		//! \brief Maximum order of field expansion.
	int lm;		int lm;
	//! \brief NLMMT = 2 * LM * (LM + 2)
	int nlmmt;
	//! \brief Contiguous RMI vector.		//! \brief Contiguous RMI vector.
	dcomplex *vec_rmi;		dcomplex *vec_rmi;
	//! \brief Contiguous REI vector.		//! \brief Contiguous REI vector.
	@@ -52,6 +50,10 @@ protected:
	dcomplex *vec_vints;		dcomplex *vec_vints;

	public:		public:
			//! \brief NLMMT = 2 * LM * (LM + 2).
			int nlmmt;
			//! \brief Number of configurations
			int configurations;
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	dcomplex **rmi;		dcomplex **rmi;
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	@@ -104,6 +106,12 @@ public:
	*/		*/
	C1(int ns, int l_max, int nshl, int iog);		C1(int ns, int l_max, int nshl, int iog);

			/*! \brief C1 instance constructor copying all contents from a preexisting template
			*
			* \param rhs: `C1` preexisting template.
			*/
			C1(const C1& rhs);

	//! \brief C1 instance destroyer.		//! \brief C1 instance destroyer.
	~C1();		~C1();
	};		};
	@@ -112,10 +120,14 @@ public:
	*		*
	*/		*/
	class C2 {		class C2 {
			protected:
	//! \brief Number of spheres.		//! \brief Number of spheres.
	int nsph;		int nsph;
	//! \brief Number of required orders.		//! \brief Number of required orders.
	int nhspo;		int nhspo;
			//! \brief QUESTION: what is nl?
			int nl;

	public:		public:
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	dcomplex *ris;		dcomplex *ris;
	@@ -137,6 +149,12 @@ public:
	*/		*/
	C2(int ns, int nl, int npnt, int npntts);		C2(int ns, int nl, int npnt, int npntts);

			/*! \brief C2 instance constructor copying its contents from preexisting instance.
			*
			* \param rhs: `C2` object to copy contents from
			*/
			C2(const C2& rhs);

	//! \brief C2 instance destroyer.		//! \brief C2 instance destroyer.
	~C2();		~C2();
	};		};
	@@ -162,6 +180,10 @@ public:
	*/		*/
	C3();		C3();

			/*! \brief C3 instance constructor copying its contents from a preexisting object.
			*/
			C3(const C3& rhs);

	/*! \brief C3 instance destroyer.		/*! \brief C3 instance destroyer.
	*/		*/
	~C3();		~C3();
	@@ -169,7 +191,8 @@ public:

	/*! \brief Representation of the FORTRAN C4 blocks.		/*! \brief Representation of the FORTRAN C4 blocks.
	*/		*/
	struct C4 {		class C4 {
			public:
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	int litpo;		int litpo;
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	@@ -194,6 +217,17 @@ struct C4 {
	int nsph;		int nsph;
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	int nv3j;		int nv3j;

			/*! \brief C3 instance constructor.
			*/
			C4(int li, int le, int nsph);
			/*! \brief C3 instance constructor copying its contents from a preexisting object.
			*/
			C4(const C4& rhs);

			/*! \brief C3 instance destroyer.
			*/
			~C4();
	};		};

	/*! \brief Vectors and matrices that are specific to cluster C1 blocks.		/*! \brief Vectors and matrices that are specific to cluster C1 blocks.
	@@ -207,6 +241,18 @@ protected:
	int nlemt;		int nlemt;
	//! \brief Maximum expansion order plus one. QUESTION: correct?		//! \brief Maximum expansion order plus one. QUESTION: correct?
	int lmpo;		int lmpo;
			//! \brief QUESTION: definition?
			int litpo;
			//! \brief QUESTION: definition?
			int lmtpo;
			//! \brief QUESTION: definition?
			int litpos;
			//! \brief QUESTION: definition?
			int lmtpos;
			//! \brief QUESTION: definition?
			int nv3j;
			//! \brief QUESTION: definition?
			int lm;

	public:		public:
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	@@ -224,6 +270,10 @@ public:
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	dcomplex *vintm;		dcomplex *vintm;
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
			dcomplex *vint;
			//! \brief QUESTION: definition?
			dcomplex *vintm;
			//! \brief QUESTION: definition?
	dcomplex *vintt;		dcomplex *vintt;
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	dcomplex **fsac;		dcomplex **fsac;
	@@ -260,6 +310,12 @@ public:
	*/		*/
	C1_AddOns(C4 *c4);		C1_AddOns(C4 *c4);

			/*! \brief C1_AddOns instance constructor copying contents from a preexisting object.
			*
			* \param rhs: `C1_AddOns` preexisting object to copy from.
			*/
			C1_AddOns(const C1_AddOns& rhs);

	//! \brief C1_AddOns instance destroyer.		//! \brief C1_AddOns instance destroyer.
	~C1_AddOns();		~C1_AddOns();
	};		};
	@@ -268,6 +324,8 @@ public:
	*/		*/
	class C6 {		class C6 {
	public:		public:
			//! \brief LMTPO = 2 * LM + 1.
			int lmtpo;
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	double *rac3j;		double *rac3j;

	@@ -277,6 +335,12 @@ public:
	*/		*/
	C6(int lmtpo);		C6(int lmtpo);

			/*! \brief C6 instance constructor copying contents from preexisting object.
			*
			* \param lmtpo: `int` QUESTION: definition?
			*/
			C6(const C6& rhs);

	/*! \brief C6 instance destroyer.		/*! \brief C6 instance destroyer.
	*/		*/
	~C6();		~C6();
	@@ -290,6 +354,11 @@ protected:
	int gis_size_0;		int gis_size_0;
	//! \brief Number of rows in the SAM matrix		//! \brief Number of rows in the SAM matrix
	int sam_size_0;		int sam_size_0;
			//! \brief QUESTION: definition?
			int nlem;
			//! \brief QUESTION: definition?
			int nlemt;

	public:		public:
	//! \brief QUESTION: definition?		//! \brief QUESTION: definition?
	dcomplex **gis;		dcomplex **gis;
	@@ -307,6 +376,12 @@ public:
	*/		*/
	C9(int ndi, int nlem, int ndit, int nlemt);		C9(int ndi, int nlem, int ndit, int nlemt);

			/*! \brief C9 instance constructor copying contents from preexisting object.
			*
			* \param rhs: `C9` preexisting object to copy from
			*/
			C9(const C9& rhs);

	/*! \brief C9 instance destroyer.		/*! \brief C9 instance destroyer.
	*/		*/
	~C9();		~C9();

src/include/Configuration.h

+46 −1

Original line number	Original line	Diff line number	Diff line
	@@ -30,6 +30,18 @@
	#ifndef INCLUDE_CONFIGURATION_H_		#ifndef INCLUDE_CONFIGURATION_H_
	#define INCLUDE_CONFIGURATION_H_		#define INCLUDE_CONFIGURATION_H_

			class ScattererConfiguration;
			class GeometryConfiguration;
			class C1;
			class C1_AddOns;
			class C2;
			class C3;
			class C4;
			class C6;
			class C9;

			#include <fstream>

	/**		/**
	* \brief A class to represent the configuration of the scattering geometry.		* \brief A class to represent the configuration of the scattering geometry.
	*		*
	@@ -41,6 +53,7 @@
	class GeometryConfiguration {		class GeometryConfiguration {
	//! Temporary work-around to allow cluster() and sphere() peeking in.		//! Temporary work-around to allow cluster() and sphere() peeking in.
	friend void cluster(std::string, std::string, std::string);		friend void cluster(std::string, std::string, std::string);
			friend int cluster_jxi488_cycle(int, ScattererConfiguration , GeometryConfiguration , C1 , C1_AddOns , C2 , C3 , C4 , C6 , C9 , std::FILE , std::string, double , double , dcomplex , double , dcomplex , double *, double , dcomplex *, int, int, int, int, int, int, int, double , double , double , dcomplex , double , dcomplex , double , dcomplex , double , std::fstream &, double , double , double , double , double , double , double , int, int, np_int, int, int, int, int, int, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double , double , double , double , double , double , double , double , double , double , double &, double &, double &, double &, double &, double, double, int, dcomplex, double, double, np_int, dcomplex *, int, int);
	friend void sphere(std::string, std::string, std::string);		friend void sphere(std::string, std::string, std::string);
	protected:		protected:
	//! \brief Number of spherical components.		//! \brief Number of spherical components.
	@@ -139,6 +152,11 @@ public:
	double sc_ph_start, double sc_ph_step, double sc_ph_end,		double sc_ph_start, double sc_ph_step, double sc_ph_end,
	int jwtm		int jwtm
	);		);
			/*! \brief Build a scattering geometry configuration structure copying it from an existing one.
			*
			* \param rhs: `GeometryConfiguration` preexisting object to copy from.
			*/
			GeometryConfiguration(const GeometryConfiguration& rhs);

	/*! \brief Destroy a GeometryConfiguration instance.		/*! \brief Destroy a GeometryConfiguration instance.
	*/		*/
	@@ -167,6 +185,7 @@ public:
	class ScattererConfiguration {		class ScattererConfiguration {
	//! Temporary work-around to allow cluster() and sphere() peeking in.		//! Temporary work-around to allow cluster() and sphere() peeking in.
	friend void cluster(std::string, std::string, std::string);		friend void cluster(std::string, std::string, std::string);
			friend int cluster_jxi488_cycle(int, ScattererConfiguration , GeometryConfiguration , C1 , C1_AddOns , C2 , C3 , C4 , C6 , C9 , std::FILE , std::string, double , double , dcomplex , double , dcomplex , double *, double , dcomplex *, int, int, int, int, int, int, int, double , double , double , dcomplex , double , dcomplex , double , dcomplex , double , std::fstream &, double , double , double , double , double , double , double , int, int, np_int, int, int, int, int, int, int, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double, double , double , double , double , double , double , double , double , double , double , double &, double &, double &, double &, double &, double, double, int, dcomplex, double, double, np_int, dcomplex *, int, int);
	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].
	@@ -185,9 +204,11 @@ protected:
	std::string reference_variable_name;		std::string reference_variable_name;
	//! \brief Number of spherical components.		//! \brief Number of spherical components.
	int number_of_spheres;		int number_of_spheres;
			//! \brief Number of configurations.
			int configurations;
	//! \brief Number of scales to use in calculation.		//! \brief Number of scales to use in calculation.
	int number_of_scales;		int number_of_scales;
	//! \brief Type of dielectric functions (<0 at XIP, =0 as function of XI, >0 contants).		//! \brief Type of dielectric functions (<0 at XIP, =0 as function of XI, >0 constants).
	int idfc;		int idfc;
	//! \brief External medium dielectric constant. QUESTION: correct?		//! \brief External medium dielectric constant. QUESTION: correct?
	double exdc;		double exdc;
	@@ -268,6 +289,7 @@ public:
	*/		*/
	ScattererConfiguration(		ScattererConfiguration(
	int nsph,		int nsph,
			int configs,
	double *scale_vector,		double *scale_vector,
	int nxi,		int nxi,
	std::string variable_name,		std::string variable_name,
	@@ -282,6 +304,29 @@ public:
	double wp,		double wp,
	double xip		double xip
	);		);
			/*! \brief Build a scatterer configuration structure copying its contents from a preexisting one.
			*
			* Prepare a default configuration structure by allocating the necessary
			* memory structures.
			*
			* \param nsph: `int` The number of spheres in the simulation.
			* \param scale_vector: `double*` The radiation-particle scale vector.
			* \param nxi: `int` The number of radiation-particle scalings.
			* \param variable_name: `string` The name of the radiation-particle scaling type.
			* \param iog_vector: `int*` Array of sphere identification numbers. QUESTION: correct?
			* \param ros_vector: `double*` Sphere radius array.
			* \param nshl_vector: `int*` Array of layer numbers.
			* \param rcf_vector: `double**` Array of fractional break radii. QUESTION: correct?
			* \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
			* for dimensions).
			* \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 exdc: `double` External medium dielectric constant.
			* \param wp: `double` wp
			* \param xip: `double` xip
			*/
			ScattererConfiguration(const ScattererConfiguration& rhs);

	/*! \brief Destroy a scatterer configuration instance.		/*! \brief Destroy a scatterer configuration instance.
	*/		*/

src/libnptm/Commons.cpp

+284 −19

Original line number	Original line	Diff line number	Diff line
	@@ -36,15 +36,15 @@ C1::C1(int ns, int l_max, int _nshl, int _iog) {
	vec_w = new dcomplex[4 * nlmmt]();		vec_w = new dcomplex[4 * nlmmt]();
	w = new dcomplex*[nlmmt];		w = new dcomplex*[nlmmt];
	for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);		for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);
	int configurations = 0;		configurations = 0;
	for (int ci = 1; ci <= nsph; ci++) {		for (int ci = 1; ci <= nsph; ci++) {
	if (_iog[ci - 1] >= ci) configurations++;		if (_iog[ci - 1] >= ci) configurations++;
	}		}
	vint = new dcomplex[16]();		vint = new dcomplex[16]();
	vec_vints = new dcomplex[nsph * 16]();		vec_vints = new dcomplex[nsph * 16]();
	vints = new dcomplex*[nsph];		vints = new dcomplex*[nsph];
	rc = new double*[configurations];		rc = new double*[nsph];
	nshl = new int[configurations]();		nshl = new int[nsph]();
	iog = new int[nsph]();		iog = new int[nsph]();
	int conf_index = 0;		int conf_index = 0;
	for (int vi = 0; vi < nsph; vi++) {		for (int vi = 0; vi < nsph; vi++) {
	@@ -77,6 +77,84 @@ C1::C1(int ns, int l_max, int _nshl, int _iog) {
	}		}
	}		}

			C1::C1(const C1& rhs) {
			nlmmt = rhs.nlmmt;
			nsph = rhs.nsph;
			lm = rhs.lm;

			vec_rmi = new dcomplex[nsph * lm]();
			vec_rei = new dcomplex[nsph * lm]();
			rmi = new dcomplex*[lm];
			rei = new dcomplex*[lm];
			for (int ri = 0; ri < lm; ri++) {
			rmi[ri] = &(vec_rmi[nsph * ri]);
			rei[ri] = &(vec_rei[nsph * ri]);
			/! Copy the contents from the template /
			for (int rj=0; rj<nsph; rj++) {
			rmi[ri][rj] = rhs.rmi[ri][rj];
			rei[ri][rj] = rhs.rei[ri][rj];
			}
			}
			vec_w = new dcomplex[4 * nlmmt]();
			w = new dcomplex*[nlmmt];
			for (int wi = 0; wi < nlmmt; wi++) {
			w[wi] = &(vec_w[4 * wi]);
			for (int wj=0; wj<4; wj++) w[wi][wj] = rhs.w[wi][wj];
			}
			configurations = rhs.configurations;
			vec_vints = new dcomplex[nsph * 16]();
			vints = new dcomplex*[nsph];
			rc = new double*[nsph];
			nshl = new int[nsph]();
			for (int ni = 0; ni < nsph; ni++) nshl[ni] = rhs.nshl[ni];
			iog = new int[nsph]();
			for (int ni = 0; ni < nsph; ni++) iog[ni] = rhs.iog[ni];
			int conf_index = 0;
			for (int vi = 0; vi < nsph; vi++) {
			vints[vi] = &(vec_vints[16 * vi]);
			for (int vj = 0; vj < 16; vj++) vints[vi][vj] = rhs.vints[vi][vj];
			}
			for (int ri=0; ri<configurations; ri++) {
			rc[ri] = new double[nshl[ri]]();
			for (int rj=0; rj<nshl[ri]; rj++) rc[ri][rj] = rhs.rc[ri][rj];
			}
			fsas = new dcomplex[nsph]();
			sscs = new double[nsph]();
			sexs = new double[nsph]();
			sabs = new double[nsph]();
			sqscs = new double[nsph]();
			sqexs = new double[nsph]();
			sqabs = new double[nsph]();
			gcsv = new double[nsph]();
			rxx = new double[nsph]();
			ryy = new double[nsph]();
			rzz = new double[nsph]();
			ros = new double[nsph]();

			sas = new dcomplex**[nsph];
			for (int si = 0; si < nsph; si++) {
			sas[si] = new dcomplex*[2];
			for (int sj=0; sj<2; sj++) {
			sas[si][sj] = new dcomplex[2]();
			for (int sk=0; sk<2; sk++) sas[si][sj][sk] = rhs.sas[si][sj][sk];
			}
			fsas[si] = rhs.fsas[si];
			sscs[si] = rhs.sscs[si];
			sexs[si] = rhs.sexs[si];
			sabs[si] = rhs.sabs[si];
			sqscs[si] = rhs.sqscs[si];
			sqexs[si] = rhs.sqexs[si];
			sqabs[si] = rhs.sqabs[si];
			gcsv[si] = rhs.gcsv[si];
			rxx[si] = rhs.rxx[si];
			ryy[si] = rhs.ryy[si];
			rzz[si] = rhs.rzz[si];
			}
			for (int si = 0; si < configurations; si++) {
			ros[si] = rhs.ros[si];
			}
			}

	C1::~C1() {		C1::~C1() {
	delete[] vec_rmi;		delete[] vec_rmi;
	delete[] vec_rei;		delete[] vec_rei;
	@@ -85,11 +163,8 @@ C1::~C1() {
	delete[] vec_w;		delete[] vec_w;
	delete[] w;		delete[] w;
	int conf_index = 0;		int conf_index = 0;
	for (int ci = 1; ci <= nsph; ci++) {		for (int ci = 0; ci < configurations; ci++) {
	if (iog[ci] >= ci) {		delete[] rc[ci];
	delete[] rc[conf_index];
	conf_index++;
	}
	}		}
	delete[] rc;		delete[] rc;
	delete[] vint;		delete[] vint;
	@@ -121,15 +196,22 @@ C1_AddOns::C1_AddOns(C4 *c4) {
	nsph = c4->nsph;		nsph = c4->nsph;
	lmpo = c4->lmpo;		lmpo = c4->lmpo;
	nlemt = 2 * c4->nlem;		nlemt = 2 * c4->nlem;
	vh = new dcomplex[(nsph * nsph - 1) * c4->litpo]();		litpo = c4->litpo;
	vj0 = new dcomplex[nsph * c4->lmtpo]();		lmtpo = c4->lmtpo;
	vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?		litpos = c4->litpos;
	vyhj = new dcomplex[(nsph * nsph - 1) * c4->litpos]();		lmtpos = c4->lmtpos;
	vyj0 = new dcomplex[nsph * c4->lmtpos]();		nv3j = c4->nv3j;
			lm = c4->lm;
			vh = new dcomplex[(nsph * nsph - 1) * litpo]();
			vj0 = new dcomplex[nsph * lmtpo]();
			vj = new dcomplex[1]();
			vyhj = new dcomplex[(nsph * nsph - 1) * litpos]();
			vyj0 = new dcomplex[nsph * lmtpos]();
	am0m = new dcomplex*[nlemt];		am0m = new dcomplex*[nlemt];
	for (int ai = 0; ai < nlemt; ai++) {		for (int ai = 0; ai < nlemt; ai++) {
	am0m[ai] = new dcomplex[nlemt]();		am0m[ai] = new dcomplex[nlemt]();
	}		}
			vint = new dcomplex[16]();
	vintm = new dcomplex[16]();		vintm = new dcomplex[16]();
	vintt = new dcomplex[16]();		vintt = new dcomplex[16]();
	// This calculates the size of v3j0		// This calculates the size of v3j0
	@@ -149,11 +231,93 @@ C1_AddOns::C1_AddOns(C4 *c4) {
	ecscp = new dcomplex[2]();		ecscp = new dcomplex[2]();
	scscpm = new dcomplex[2]();		scscpm = new dcomplex[2]();
	ecscpm = new dcomplex[2]();		ecscpm = new dcomplex[2]();
			v3j0 = new double[nv3j]();
			ind3j = new int*[lmpo];
			for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[lm]();
			sscs = new double[nsph]();
			ecscm = new double[2]();
			scscm = new double[2]();
			scsc = new double[2]();
			ecsc = new double[2]();
			}

			C1_AddOns::C1_AddOns(const C1_AddOns& rhs) {
			nsph = rhs.nsph;
			lmpo = rhs.lmpo;
			nlemt = rhs.nlemt;
			litpo = rhs.litpo;
			lmtpo = rhs.lmtpo;
			litpos = rhs.litpos;
			lmtpos = rhs.lmtpos;
			nv3j = rhs.nv3j;
			lm = rhs.lm;
			int vhsize=(nsph * nsph - 1) * litpo;
			vh = new dcomplex[vhsize]();
			for (int hi=0; hi<vhsize; hi++) vh[hi] = rhs.vh[hi];
			int vj0size=nsph * lmtpo;
			vj0 = new dcomplex[vj0size]();
			for (int hi=0; hi<vj0size; hi++) vj0[hi] = rhs.vj0[hi];
			vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?
			vj[0] = rhs.vj[0];
			int vyhjsize=(nsph * nsph - 1) * litpos;
			vyhj = new dcomplex[vyhjsize]();
			for (int hi=0; hi<vyhjsize; hi++) vyhj[hi] = rhs.vyhj[hi];
			int vyj0size = nsph * lmtpos;
			vyj0 = new dcomplex[vyj0size]();
			for (int hi=0; hi<vyj0size; hi++) vyj0[hi] = rhs.vyj0[hi];
			am0m = new dcomplex*[nlemt];
			for (int ai = 0; ai < nlemt; ai++) {
			am0m[ai] = new dcomplex[nlemt]();
			for (int aj = 0; aj < nlemt; aj++) am0m[ai][aj] = rhs.am0m[ai][aj];
			}
			vint = new dcomplex[16](); // QUESTION: is dimension 16 generally fixed?
			vintm = new dcomplex[16]();
			vintt = new dcomplex[16]();
			for (int hi = 0; hi < 16; hi++) {
			vint[hi] = rhs.vint[hi];
			vintm[hi] = rhs.vintm[hi];
			vintt[hi] = rhs.vintt[hi];
			}
			fsac = new dcomplex*[2];
			sac = new dcomplex*[2];
			fsacm = new dcomplex*[2];
			for (int fi = 0; fi < 2; fi++) {
			fsac[fi] = new dcomplex[2]();
			sac[fi] = new dcomplex[2]();
			fsacm[fi] = new dcomplex[2]();
			for (int fj = 0; fj < 2; fj++) {
			fsac[fi][fj] = rhs.fsac[fi][fj];
			sac[fi][fj] = rhs.sac[fi][fj];
			fsacm[fi][fj] = rhs.fsacm[fi][fj];
			}
			}
			scscp = new dcomplex[2]();
			ecscp = new dcomplex[2]();
			scscpm = new dcomplex[2]();
			ecscpm = new dcomplex[2]();
			v3j0 = new double[nv3j]();
			for (int hi = 0; hi < nv3j; hi++) v3j0[hi] = rhs.v3j0[hi];
			ind3j = new int*[lmpo];
			for (int ii = 0; ii < lmpo; ii++) {
			ind3j[ii] = new int[lm]();
			for (int ij=0; ij<lm; ij++) ind3j[ii][ij] = rhs.ind3j[ii][ij];
			}
	sscs = new double[nsph]();		sscs = new double[nsph]();
			for (int hi=0; hi<nsph; hi++) sscs[hi] = rhs.sscs[hi];
	ecscm = new double[2]();		ecscm = new double[2]();
	scscm = new double[2]();		scscm = new double[2]();
	scsc = new double[2]();		scsc = new double[2]();
	ecsc = new double[2]();		ecsc = new double[2]();
			for (int hi=0; hi<2; hi++) {
			scscp[hi] = rhs.scscp[hi];
			ecscp[hi] = rhs.ecscp[hi];
			scscpm[hi] = rhs.scscpm[hi];
			ecscpm[hi] = rhs.ecscpm[hi];
			ecscm[hi] = rhs.ecscm[hi];
			scscm[hi] = rhs.scscm[hi];
			scsc[hi] = rhs.scsc[hi];
			ecsc[hi] = rhs.ecsc[hi];
			}
	}		}

	C1_AddOns::~C1_AddOns() {		C1_AddOns::~C1_AddOns() {
	@@ -168,7 +332,8 @@ C1_AddOns::~C1_AddOns() {
	for (int ai = nlemt - 1; ai > -1; ai--) {		for (int ai = nlemt - 1; ai > -1; ai--) {
	delete[] am0m[ai];		delete[] am0m[ai];
	}		}
	delete am0m;		delete[] am0m;
			delete[] vint;
	delete[] vintm;		delete[] vintm;
	delete[] vintt;		delete[] vintt;
	for (int fi = 1; fi > -1; fi--) {		for (int fi = 1; fi > -1; fi--) {
	@@ -189,10 +354,11 @@ C1_AddOns::~C1_AddOns() {
	delete[] ecsc;		delete[] ecsc;
	}		}

	C2::C2(int ns, int nl, int npnt, int npntts) {		C2::C2(int ns, int _nl, int npnt, int npntts) {
	nsph = ns;		nsph = ns;
	int max_n = (npnt > npntts) ? npnt : npntts;		int max_n = (npnt > npntts) ? npnt : npntts;
	nhspo = 2 * max_n - 1;		nhspo = 2 * max_n - 1;
			nl = _nl;
	ris = new dcomplex[nhspo]();		ris = new dcomplex[nhspo]();
	dlri = new dcomplex[nhspo]();		dlri = new dcomplex[nhspo]();
	vkt = new dcomplex[nsph]();		vkt = new dcomplex[nsph]();
	@@ -200,6 +366,26 @@ C2::C2(int ns, int nl, int npnt, int npntts) {
	vsz = new double[nsph]();		vsz = new double[nsph]();
	}		}

			C2::C2(const C2& rhs) {
			nsph = rhs.nsph;
			nhspo = rhs.nhspo;
			nl = rhs.nl;
			ris = new dcomplex[nhspo]();
			dlri = new dcomplex[nhspo]();
			for (int ind=0; ind<nhspo; ind++) {
			ris[ind] = rhs.ris[ind];
			dlri[ind] = rhs.dlri[ind];
			}
			vkt = new dcomplex[nsph]();
			vsz = new double[nsph]();
			for (int ind=0; ind<nsph; ind++) {
			vkt[ind] = rhs.vkt[ind];
			vsz[ind] = rhs.vsz[ind];
			}
			dc0 = new dcomplex[nl]();
			for (int ind=0; ind<nl; ind++) dc0[ind] = rhs.dc0[ind];
			}

	C2::~C2() {		C2::~C2() {
	delete[] ris;		delete[] ris;
	delete[] dlri;		delete[] dlri;
	@@ -219,21 +405,80 @@ C3::C3() {
	acs = 0.0;		acs = 0.0;
	}		}

			C3::C3(const C3& rhs) {
			tsas = new dcomplex*[2];
			for (int ti=0; ti<2; ti++) {
			tsas[ti] = new dcomplex[2];
			for (int tj=0; tj<2; tj++) tsas[ti][tj] = rhs.tsas[ti][tj];
			}
			tfsas = rhs.tfsas;
			gcs = rhs.gcs;
			scs = rhs.scs;
			ecs = rhs.ecs;
			acs = rhs.acs;
			}

	C3::~C3() {		C3::~C3() {
	delete[] tsas[1];		delete[] tsas[1];
	delete[] tsas[0];		delete[] tsas[0];
	delete[] tsas;		delete[] tsas;
	}		}

	C6::C6(int lmtpo) {		C4::C4(int _li, int _le, int _nsph) {
			li = _li;
			le = _le;
			lm = (li > le) ? li : le;
			nv3j = (lm * (lm + 1) * (2 * lm + 7)) / 6;
			nsph = _nsph;
			// The following is needed to initialize C1_AddOns
			litpo = li + li + 1;
			litpos = litpo * litpo;
			lmtpo = li + le + 1;
			lmtpos = lmtpo * lmtpo;
			nlim = li * (li + 2);
			nlem = le * (le + 2);
			lmpo = lm + 1;
			}

			C4::C4(const C4& rhs) {
			li = rhs.li;
			le = rhs.le;
			lm = rhs.lm;
			nv3j = rhs.nv3j;
			nsph = rhs.nsph;
			// The following is needed to initialize C1_AddOns
			litpo = rhs.litpo;
			litpos = rhs.litpos;
			lmtpo = rhs.lmtpo;
			lmtpos = rhs.lmtpos;
			nlim = rhs.nlim;
			nlem = rhs.nlem;
			lmpo = rhs.lmpo;
			}

			C4::~C4() {
			}

			C6::C6(int _lmtpo) {
			lmtpo = _lmtpo;
	rac3j = new double[lmtpo]();		rac3j = new double[lmtpo]();
	}		}

			C6::C6(const C6& rhs) {
			lmtpo = rhs.lmtpo;
			rac3j = new double[lmtpo]();
			for (int ri = 0; ri < lmtpo; ri++) rac3j[ri] = rhs.rac3j[ri];
			}

	C6::~C6() {		C6::~C6() {
	delete[] rac3j;		delete[] rac3j;
	}		}

	C9::C9(int ndi, int nlem, int ndit, int nlemt) {		C9::C9(int ndi, int _nlem, int ndit, int _nlemt) {
			gis_size_0 = ndi;
			sam_size_0 = ndit;
			nlem = _nlem;
			nlemt = _nlemt;
	gis = new dcomplex*[ndi];		gis = new dcomplex*[ndi];
	gls = new dcomplex*[ndi];		gls = new dcomplex*[ndi];
	for (int gi = 0; gi < ndi; gi++) {		for (int gi = 0; gi < ndi; gi++) {
	@@ -242,8 +487,28 @@ C9::C9(int ndi, int nlem, int ndit, int nlemt) {
	}		}
	sam = new dcomplex*[ndit];		sam = new dcomplex*[ndit];
	for (int si = 0; si < ndit; si++) sam[si] = new dcomplex[nlemt]();		for (int si = 0; si < ndit; si++) sam[si] = new dcomplex[nlemt]();
	gis_size_0 = ndi;		}
	sam_size_0 = ndit;
			C9::C9(const C9& rhs) {
			gis_size_0 = rhs.gis_size_0;
			sam_size_0 = rhs.sam_size_0;
			nlem = rhs.nlem;
			nlemt = rhs.nlemt;
			gis = new dcomplex*[gis_size_0];
			gls = new dcomplex*[gis_size_0];
			for (int gi = 0; gi < gis_size_0; gi++) {
			gis[gi] = new dcomplex[nlem]();
			gls[gi] = new dcomplex[nlem]();
			for (int gj=0; gj<nlem; gj++) {
			gis[gi][gj] = rhs.gis[gi][gj];
			gls[gi][gj] = rhs.gls[gi][gj];
			}
			}
			sam = new dcomplex*[sam_size_0];
			for (int si = 0; si < sam_size_0; si++) {
			sam[si] = new dcomplex[nlemt]();
			for (int sj=0; sj<nlemt; sj++) sam[si][sj] = rhs.sam[si][sj];
			}
	}		}

	C9::~C9() {		C9::~C9() {

src/libnptm/Configuration.cpp

+82 −9

Original line number	Original line	Diff line number	Diff line
	@@ -76,6 +76,40 @@ GeometryConfiguration::GeometryConfiguration(
	sph_y = y;		sph_y = y;
	sph_z = z;		sph_z = z;
	}		}
			GeometryConfiguration::GeometryConfiguration(const GeometryConfiguration& rhs)
			{
			number_of_spheres = rhs.number_of_spheres;
			l_max = rhs.l_max;
			in_pol = rhs.in_pol;
			npnt = rhs.npnt;
			npntts = rhs.npntts;
			meridional_type = rhs.meridional_type;
			li = rhs.li;
			le = rhs.le;
			mxndm = rhs.mxndm;
			iavm = rhs.iavm;
			in_theta_start = rhs.in_theta_start;
			in_theta_step = rhs.in_theta_step;
			in_theta_end = rhs.in_theta_end;
			in_phi_start = rhs.in_phi_start;
			in_phi_step = rhs.in_phi_step ;
			in_phi_end = rhs.in_phi_end;
			sc_theta_start = rhs.sc_theta_start;
			sc_theta_step = rhs.sc_theta_step;
			sc_theta_end = rhs.sc_theta_end;
			sc_phi_start = rhs.sc_phi_start;
			sc_phi_step = rhs.sc_phi_step;
			sc_phi_end = rhs.sc_phi_end;
			jwtm = rhs.jwtm;
			sph_x = new double[number_of_spheres]();
			sph_y = new double[number_of_spheres]();
			sph_z = new double[number_of_spheres]();
			for (int ni=0; ni<number_of_spheres; ni++) {
			sph_x[ni] = rhs.sph_x[ni];
			sph_y[ni] = rhs.sph_y[ni];
			sph_z[ni] = rhs.sph_z[ni];
			}
			}

	GeometryConfiguration::~GeometryConfiguration() {		GeometryConfiguration::~GeometryConfiguration() {
	delete[] sph_x;		delete[] sph_x;
	@@ -199,6 +233,7 @@ GeometryConfiguration* GeometryConfiguration::from_legacy(string file_name) {

	ScattererConfiguration::ScattererConfiguration(		ScattererConfiguration::ScattererConfiguration(
	int nsph,		int nsph,
			int configs,
	double *scale_vector,		double *scale_vector,
	int nxi,		int nxi,
	string variable_name,		string variable_name,
	@@ -214,6 +249,7 @@ ScattererConfiguration::ScattererConfiguration(
	double x		double x
	) {		) {
	number_of_spheres = nsph;		number_of_spheres = nsph;
			configurations = configs;
	number_of_scales = nxi;		number_of_scales = nxi;
	reference_variable_name = variable_name;		reference_variable_name = variable_name;
	iog_vec = iog_vector;		iog_vec = iog_vector;
	@@ -241,6 +277,40 @@ ScattererConfiguration::ScattererConfiguration(
	}		}
	}		}
	}		}
			ScattererConfiguration::ScattererConfiguration(const ScattererConfiguration& rhs)
			{
			number_of_spheres = rhs.number_of_spheres;
			configurations = rhs.configurations;
			number_of_scales = rhs.number_of_scales;
			reference_variable_name = rhs.reference_variable_name;
			idfc = rhs.idfc;
			use_external_sphere = rhs.use_external_sphere;
			exdc = rhs.exdc;
			wp = rhs.wp;
			xip = rhs.xip;
			iog_vec = new int[number_of_spheres]();
			radii_of_spheres = new double[number_of_spheres]();
			nshl_vec = new int[number_of_spheres]();
			rcf = new double*[number_of_spheres];
			scale_vec = new double[number_of_scales]();
			dc0_matrix = new dcomplex**[configurations];
			for (int si = 0; si < number_of_scales; si++) scale_vec[si] = rhs.scale_vec[si];
			for (int si=0; si<number_of_spheres; si++) {
			iog_vec[si] = rhs.iog_vec[si];
			}
			int dim3 = (idfc == 0) ? number_of_scales : 1;
			for (int si=0; si<configurations; si++) {
			radii_of_spheres[si] = rhs.radii_of_spheres[si];
			nshl_vec[si] = rhs.nshl_vec[si];
			rcf[si] = new double[nshl_vec[si]]();
			dc0_matrix[si] = new dcomplex*[number_of_spheres];
			for (int sj=0; sj<nshl_vec[si]; sj++) rcf[si][sj] = rhs.rcf[si][sj];
			for (int sj=0; sj<number_of_spheres; sj++) {
			dc0_matrix[si][sj] = new dcomplex[dim3]();
			for (int sk=0; sk<dim3; sk++) dc0_matrix[si][sj][sk] = rhs.dc0_matrix[si][sj][sk];
			}
			}
			}

	ScattererConfiguration::~ScattererConfiguration() {		ScattererConfiguration::~ScattererConfiguration() {
	int configurations = 0;		int configurations = 0;
	@@ -453,9 +523,9 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
	}		}
	}		}
	int index = 0;		int index = 0;
	double *ros_vector = new double[configurations]();		double *ros_vector = new double[nsph]();
	int *nshl_vector = new int[configurations]();		int *nshl_vector = new int[nsph]();
	double *rcf_vector = new double[configurations];		double *rcf_vector = new double[nsph];
	for (int i113 = 1; i113 <= nsph; i113++) {		for (int i113 = 1; i113 <= nsph; i113++) {
	if (iog_vector[i113 - 1] < i113) continue;		if (iog_vector[i113 - 1] < i113) continue;
	str_target = file_lines[++last_read_line];		str_target = file_lines[++last_read_line];
	@@ -519,6 +589,7 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam

	ScattererConfiguration *config = new ScattererConfiguration(		ScattererConfiguration *config = new ScattererConfiguration(
	nsph,		nsph,
			configurations,
	variable_vector,		variable_vector,
	nxi,		nxi,
	variable_name,		variable_name,
	@@ -568,9 +639,9 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
	if (iog[ci] < ci + 1) continue;		if (iog[ci] < ci + 1) continue;
	configuration_count++;		configuration_count++;
	}		}
	nshl_vector = new int[configuration_count]();		nshl_vector = new int[nsph]();
	ros_vector = new double[configuration_count]();		ros_vector = new double[nsph]();
	rcf_vector = new double*[configuration_count];		rcf_vector = new double*[nsph];
	for (int i115 = 1; i115 <= nsph; i115++) {		for (int i115 = 1; i115 <= nsph; i115++) {
	if (iog[i115 - 1] < i115) continue;		if (iog[i115 - 1] < i115) continue;
	str_name = "NSHL_" + to_string(iog[i115 - 1]);		str_name = "NSHL_" + to_string(iog[i115 - 1]);
	@@ -614,6 +685,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
	status = hdf_file->close();		status = hdf_file->close();
	conf = new ScattererConfiguration(		conf = new ScattererConfiguration(
	nsph,		nsph,
			configuration_count,
	xi_vec,		xi_vec,
	nxi,		nxi,
	"XIV",		"XIV",
	@@ -664,9 +736,9 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
	for (int i = 1; i <= _nsph; i++) {		for (int i = 1; i <= _nsph; i++) {
	if (_iog_vec[i - 1] >= i) configurations++;		if (_iog_vec[i - 1] >= i) configurations++;
	}		}
	_nshl_vec = new int[configurations]();		_nshl_vec = new int[_nsph]();
	_ros_vec = new double[configurations]();		_ros_vec = new double[_nsph]();
	_rcf_vec = new double*[configurations];		_rcf_vec = new double*[_nsph];
	for (int i115 = 1; i115 <= _nsph; i115++) {		for (int i115 = 1; i115 <= _nsph; i115++) {
	if (_iog_vec[i115 - 1] < i115) continue;		if (_iog_vec[i115 - 1] < i115) continue;
	input.read(reinterpret_cast<char *>(&(_nshl_vec[i115 - 1])), sizeof(int));		input.read(reinterpret_cast<char *>(&(_nshl_vec[i115 - 1])), sizeof(int));
	@@ -702,6 +774,7 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {

	ScattererConfiguration *conf = new ScattererConfiguration(		ScattererConfiguration *conf = new ScattererConfiguration(
	_nsph,		_nsph,
			configurations,
	_xi_vec,		_xi_vec,
	_nxi,		_nxi,
	"XIV",		"XIV",