Fix swapping progression on Swap1 (cfd0fa52) · Commits · Giacomo Mulas / NP_TMcode

src/include/tfrfme.h

+36 −84

Original line number	Diff line number	Diff line
		@@ -11,6 +11,8 @@
		class Swap1 {
		protected:
		//! NLMMT = 2 * LM * (LM + 2)
		int last_index;
		int nkv;
		int nlmmt;

		//! QUESTION: definition?
		@@ -46,13 +48,18 @@ public:
		/*! \brief Swap1 instance constructor.
		*
		* \param lm: `int` Maximum field expansion order.
		* \param _nkv: `int` Number of vector coordinates. QUESTION: correct?
		*/
		Swap1(int lm);
		Swap1(int lm, int _nkv);

		/*! \brief Swap1 instance destroyer.
		*/
		~Swap1() { delete[] wk; }

		/*! \brief Append an element at the end of the vector.
		*/
		void append(std::complex<double> value) { wk[last_index++] = value; }

		/*! \brief Load a Swap1 instance from binary file.
		*
		* \param file_name: `string` Name of the file.
		@@ -62,26 +69,23 @@ public:
		*/
		static Swap1* from_binary(std::string file_name, std::string mode="LEGACY");

		/*! \brief Get an element from the WK vector.
		*
		* \param index: `int` Index of the desired element.
		* \return value: `complex<double>` The value of the requested element.
		*/
		std::complex<double> get_element(int index) { return wk[index]; }

		/*! \brief Calculate the necessary amount of memory to create a new instance.
		*
		* \param lm: `int` Maximum field expansion order.
		* \param _nkv: `int` Number of vector coordinates. QUESTION: correct?
		* \return size: `long` The necessary memory size in bytes.
		*/
		static long get_memory_requirement(int lm);
		static long get_memory_requirement(int lm, int _nkv);

		/*! \brief Set an element in the WK vector.
		/*! \brief Get the pointer to the WK vector.
		*
		* \param index: `int` Index of the desired element.
		* \param value: `complex<double>` The value to be stored in the vector.
		* \return value: `complex<double> *` Memory address of the WK vector.
		*/
		std::complex<double> *get_vector() { return wk; }

		/*! \brief Bring the pointer to the next element at the start of vector.
		*/
		void set_element(int index, std::complex<double> value) { wk[index] = value; }
		void reset() { last_index = 0; }

		/*! \brief Write a Swap1 instance to binary file.
		*
		@@ -184,13 +188,11 @@ public:
		*/
		static Swap2* from_binary(std::string file_name, std::string mode="LEGACY");

		/*! \brief Get an element from the VKZM matrix.
		/*! \brief Get the pointer to the VKZM matrix.
		*
		* \param row: `int` Index of the desired element row.
		* \param col: `int` Index of the desired element column.
		* \return value: `double` The value of the requested element.
		* \return value: `double **` Pointer to the VKZM matrix.
		*/
		double get_matrix_element(int row, int col) { return vkzm[row][col]; }
		double **get_matrix() { return vkzm; }

		/*! \brief Calculate the necessary amount of memory to create a new instance.
		*
		@@ -206,20 +208,11 @@ public:
		*/
		double get_param(std::string param_name);

		/*! \brief Get an element from the VKV vector.
		/*! \brief Get the pointer to the VKV vector.
		*
		* \param index: `int` Index of the desired element.
		* \return value: `double` The value of the requested element.
		* \return value: `double *` Pointer to the VKV vector.
		*/
		double get_vector_element(int index) { return vkv[index]; }

		/*! \brief Set an element in the VKZM matrix.
		*
		* \param row: `int` Index of the desired element row.
		* \param col: `int` Index of the desired element column.
		* \param value: `double` The value of the element to be stored.
		*/
		void set_matrix_element(int row, int col, double value) { vkzm[row][col] = value; }
		double *get_vector() { return vkv; }

		/*! \brief Set a parameter by its name and value.
		*
		@@ -228,13 +221,6 @@ public:
		*/
		void set_param(std::string param_name, double value);

		/*! \brief Get an element from the VKV vector.
		*
		* \param index: `int` Index of the desired element.
		* \param value: `double` The value of element to be stored.
		*/
		void set_vector_element(int index, double value) { vkv[index] = value; }

		/*! \brief Write a Swap2 instance to binary file.
		*
		* \param file_name: `string` Name of the file.
		@@ -355,13 +341,11 @@ public:
		*/
		static TFRFME* from_binary(std::string file_name, std::string mode="LEGACY");

		/*! \brief Get an element from the WSUM matrix.
		/*! \brief Get the pointer to the WSUM matrix.
		*
		* \param row: `int` Row index of the element to be read.
		* \param col: `int` Column index of the element to be read.
		* \return value: `complex<double>` The value of the requested element.
		* \return value: `complex<double> **` Pointer to the WSUM matrix.
		*/
		std::complex<double> get_matrix_element(int row, int col);
		std::complex<double> **get_matrix() { return wsum; }

		/*! \brief Calculate the necessary amount of memory to create a new instance.
		*
		@@ -385,34 +369,23 @@ public:
		*/
		double get_param(std::string param_name);

		/*! \brief Get the X coordinate of the computed point at given index.
		/*! \brief Get the pointer to the X coordinates vector.
		*
		* \param index: `int` Index of the requested point.
		* \return x: `double` X coordinate of the requested point.
		* \return x: `double *` Pointer to X coordinates vector.
		*/
		double get_x(int index) { return xv[index]; }
		double *get_x() { return xv; }

		/*! \brief Get the Y coordinate of the computed point at given index.
		/*! \brief Get the pointer to the Y coordinates vector.
		*
		* \param index: `int` Index of the requested point.
		* \return y: `double` Y coordinate of the requested point.
		* \return y: `double *` Pointer to Y coordinates vector.
		*/
		double get_y(int index) { return yv[index]; }
		double *get_y() { return yv; }

		/*! \brief Get the Z coordinate of the computed point at given index.
		/*! \brief Get the pointer to the Z coordinates vector.
		*
		* \param index: `int` Index of the requested point.
		* \return z: `double` Z coordinate of the requested point.
		* \return z: `double *` Pointer to Z coordinates vector.
		*/
		double get_z(int index) { return zv[index]; }

		/*! \brief Set an element in the WSUM matrix.
		*
		* \param row: `int` Row index of the element to be read.
		* \param col: `int` Column index of the element to be read.
		* \param value: `complex<double>` The value to be placed in the matrix.
		*/
		void set_matrix_element(int row, int col, std::complex<double> value);
		double *get_z() { return zv; }

		/*! \brief Set a configuration parameter.
		*
		@@ -421,27 +394,6 @@ public:
		*/
		void set_param(std::string param_name, double value);

		/*! \brief Set the X coordinate of the computed point at given index.
		*
		* \param index: `int` Index of the requested point.
		* \param value: `double` X coordinate of the requested point.
		*/
		void set_x(int index, double value) { xv[index] = value; }

		/*! \brief Set the Y coordinate of the computed point at given index.
		*
		* \param index: `int` Index of the requested point.
		* \param value: `double` Y coordinate of the requested point.
		*/
		void set_y(int index, double value) { yv[index] = value; }

		/*! \brief Set the Z coordinate of the computed point at given index.
		*
		* \param index: `int` Index of the requested point.
		* \param value: `double` Z coordinate of the requested point.
		*/
		void set_z(int index, double value) { zv[index] = value; }

		/*! \brief Write a trapping configuration instance to binary file.
		*
		* \param file_name: `string` Name of the file.

src/include/tra_subs.h

+16 −16

Original line number	Diff line number	Diff line
		@@ -90,19 +90,6 @@ void ffrf(
		double fffs, CIL cil, CCR *ccr
		);

		/*! C++ porting of FFRT
		*
		* \param ac: Vector of complex. QUESTION: definition?
		* \param ws: Vector of complex. QUESTION: definition?
		* \param ffte: `double *`. QUESTION: definition?
		* \param ffts: `double *`. QUESTION: definition?
		* \param cil: `CIL *` Pointer to a CIL structure.
		*/
		void ffrt(
		std::complex<double> ac, std::complex<double> ws, double ffte, double ffts,
		CIL *cil
		);

		/*! C++ porting of FRFMER
		*
		* \param nkv: `int` QUESTION: definition?
		@@ -119,12 +106,25 @@ void ffrt(
		* \param tt1: `Swap1 *` Pointer to first swap object.
		* \param tt2: `Swap2 *` Pointer to second swap object.
		*/
		void frfmer(
		std::complex<double> *frfmer(
		int nkv, double vkm, double vknmx, double apfafa, double tra,
		double spd, double rir, double ftcn, int le, int lmode, double pmf,
		Swap1 tt1, Swap2 tt2
		);

		/*! C++ porting of FFRT
		*
		* \param ac: Vector of complex. QUESTION: definition?
		* \param ws: Vector of complex. QUESTION: definition?
		* \param ffte: `double *`. QUESTION: definition?
		* \param ffts: `double *`. QUESTION: definition?
		* \param cil: `CIL *` Pointer to a CIL structure.
		*/
		void ffrt(
		std::complex<double> ac, std::complex<double> ws, double ffte, double ffts,
		CIL *cil
		);

		/*! C++ porting of PWMALP
		*
		* \param w: Matrix of complex. QUESTION: definition?
		@@ -166,7 +166,7 @@ void sampoa(

		/*! C++ porting of WAMFF
		*
		* \param tt1: `Swap1 *` Pointer to the swap object containing the WK vector.
		* \param wk: `complex<double> *` QUESTION: definition?
		* \param x: `double`
		* \param y: `double`
		* \param z: `double`
		@@ -180,6 +180,6 @@ void sampoa(
		* \param pmf: `double` QUESTION: definition?
		*/
		void wamff(
		Swap1 *tt1, double x, double y, double z, int lm, double apfafa,
		std::complex<double> *wk, double x, double y, double z, int lm, double apfafa,
		double tra, double spd, double rir, double ftcn, int lmode, double pmf
		);

src/libnptm/tfrfme.cpp

+61 −50

Original line number	Diff line number	Diff line
		@@ -28,9 +28,12 @@
		using namespace std;

		// >>> START OF Swap1 CLASS IMPLEMENTATION <<<
		Swap1::Swap1(int lm) {
		Swap1::Swap1(int lm, int _nkv) {
		nkv = _nkv;
		nlmmt = 2 * lm * (lm + 2);
		wk = new complex<double>[nlmmt]();
		const int size = nkv * nkv * nlmmt;
		wk = new complex<double>[size]();
		last_index = 0;
		}

		Swap1* Swap1::from_binary(string file_name, string mode) {
		@@ -53,20 +56,23 @@ Swap1* Swap1::from_hdf5(string file_name) {
		herr_t status = hdf_file->get_status();
		double *elements;
		string str_type;
		int _nlmmt, lm, num_elements, index;
		int _nlmmt, _nkv, lm, num_elements, index;
		complex<double> value;
		complex<double> *_wk = NULL;
		if (status == 0) {
		status = hdf_file->read("NLMMT", "INT32", &_nlmmt);
		status = hdf_file->read("NKV", "INT32", &_nkv);
		lm = (int)((-2.0 + sqrt(4.0 + 2.0 * _nlmmt)) / 2.0);
		num_elements = 2 * _nlmmt;
		instance = new Swap1(lm);
		num_elements = 2 * _nlmmt * _nkv * _nkv;
		instance = new Swap1(lm, _nkv);
		_wk = instance->get_vector();
		elements = new double[num_elements]();
		str_type = "FLOAT64_(" + to_string(num_elements) + ")";
		status = hdf_file->read("WK", str_type, elements);
		for (int wi = 0; wi < _nlmmt; wi++) {
		for (int wi = 0; wi < num_elements / 2; wi++) {
		index = 2 * wi;
		value = complex<double>(elements[index], elements[index + 1]);
		instance->set_element(wi, value);
		_wk[wi] = value;
		} // wi loop
		delete[] elements;
		status = hdf_file->close();
		@@ -78,17 +84,21 @@ Swap1* Swap1::from_hdf5(string file_name) {
		Swap1* Swap1::from_legacy(string file_name) {
		fstream input;
		Swap1 *instance = NULL;
		int _nlmmt, lm;
		complex<double> *_wk = NULL;
		int _nlmmt, _nkv, lm;
		double rval, ival;
		input.open(file_name.c_str(), ios::in \| ios::binary);
		if (input.is_open()) {
		input.read(reinterpret_cast<char *>(&_nlmmt), sizeof(int));
		lm = (int)((-2.0 + sqrt(4.0 + 2.0 * _nlmmt)) / 2.0);
		instance = new Swap1(lm);
		for (int j = 0; j < _nlmmt; j++) {
		input.read(reinterpret_cast<char *>(&_nkv), sizeof(int));
		instance = new Swap1(lm, _nkv);
		_wk = instance->get_vector();
		int num_elements = _nlmmt * _nkv * _nkv;
		for (int j = 0; j < num_elements; j++) {
		input.read(reinterpret_cast<char *>(&rval), sizeof(double));
		input.read(reinterpret_cast<char *>(&ival), sizeof(double));
		instance->set_element(j, complex<double>(rval, ival));
		_wk[j] = complex<double>(rval, ival);
		}
		input.close();
		} else {
		@@ -97,9 +107,9 @@ Swap1* Swap1::from_legacy(string file_name) {
		return instance;
		}

		long Swap1::get_memory_requirement(int lm) {
		long size = (long)sizeof(int);
		size += (long)(sizeof(complex<double>) * 2 * lm * (lm + 2));
		long Swap1::get_memory_requirement(int lm, int _nkv) {
		long size = (long)(3 * sizeof(int));
		size += (long)(sizeof(complex<double>) * 2 * lm * (lm + 2) * _nkv * _nkv);
		return size;
		}

		@@ -120,15 +130,18 @@ void Swap1::write_hdf5(string file_name) {
		List<void *> rec_ptr_list(1);
		herr_t status;
		string str_type;
		int num_elements = 2 * nlmmt;
		int num_elements = 2 * nlmmt * nkv * nkv;
		rec_name_list.set(0, "NLMMT");
		rec_type_list.set(0, "INT32_(1)");
		rec_ptr_list.set(0, &nlmmt);
		rec_name_list.append("NKV");
		rec_type_list.append("INT32_(1)");
		rec_ptr_list.append(&nkv);
		rec_name_list.append("WK");
		str_type = "FLOAT64_(" + to_string(num_elements) + ")";
		rec_type_list.append(str_type);
		double *ptr_elements = new double[num_elements]();
		for (int wi = 0; wi < nlmmt; wi++) {
		for (int wi = 0; wi < num_elements / 2; wi++) {
		ptr_elements[2 * wi] = wk[wi].real();
		ptr_elements[2 * wi + 1] = wk[wi].imag();
		}
		@@ -156,8 +169,10 @@ void Swap1::write_legacy(string file_name) {
		double rval, ival;
		output.open(file_name.c_str(), ios::out \| ios::binary);
		if (output.is_open()) {
		int num_elements = nlmmt * nkv * nkv;
		output.write(reinterpret_cast<char *>(&nlmmt), sizeof(int));
		for (int j = 0; j < nlmmt; j++) {
		output.write(reinterpret_cast<char *>(&nkv), sizeof(int));
		for (int j = 0; j < num_elements; j++) {
		rval = wk[j].real();
		ival = wk[j].imag();
		output.write(reinterpret_cast<char *>(&rval), sizeof(double));
		@@ -173,7 +188,11 @@ bool Swap1::operator ==(Swap1 &other) {
		if (nlmmt != other.nlmmt) {
		return false;
		}
		for (int i = 0; i < nlmmt; i++) {
		if (nkv != other.nkv) {
		return false;
		}
		const int num_elements = nlmmt * nkv * nkv;
		for (int i = 0; i < num_elements; i++) {
		if (wk[i] != other.wk[i]) {
		return false;
		}
		@@ -260,19 +279,23 @@ Swap2* Swap2::from_legacy(string file_name) {
		fstream input;
		Swap2 *instance = NULL;
		int _nkv, _nlmmt, _nrvc;
		double **_vkzm = NULL;
		double *_vkv = NULL;
		double value;
		input.open(file_name.c_str(), ios::in \| ios::binary);
		if (input.is_open()) {
		input.read(reinterpret_cast<char *>(&_nkv), sizeof(int));
		instance = new Swap2(_nkv);
		_vkzm = instance->get_matrix();
		_vkv = instance->get_vector();
		for (int vj = 0; vj < _nkv; vj++) {
		input.read(reinterpret_cast<char *>(&value), sizeof(double));
		instance->set_vector_element(vj, value);
		_vkv[vj] = value;
		}
		for (int mi = 0; mi < _nkv; mi++) {
		for (int mj = 0; mj < _nkv; mj++) {
		input.read(reinterpret_cast<char *>(&value), sizeof(double));
		instance->set_matrix_element(mi, mj, value);
		_vkzm[mi][mj] = value;
		}
		}
		input.read(reinterpret_cast<char *>(&value), sizeof(double));
		@@ -589,10 +612,11 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
		unsigned int flags = H5F_ACC_RDONLY;
		HDFFile *hdf_file = new HDFFile(file_name, flags);
		herr_t status = hdf_file->get_status();
		double coord_vec, elements;
		double *elements;
		string str_type;
		int _nlmmt, _nrvc, num_elements;
		complex<double> value;
		complex<double> **_wsum = NULL;
		if (status == 0) {
		int lmode, lm, nkv, nxv, nyv, nzv;
		double vk, exri, an, ff, tra, spd, frsh, exril;
		@@ -611,6 +635,7 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
		status = hdf_file->read("FRSH", "FLOAT64", &frsh);
		status = hdf_file->read("EXRIL", "FLOAT64", &exril);
		instance = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
		_wsum = instance->get_matrix();
		instance->set_param("vk", vk);
		instance->set_param("exri", exri);
		instance->set_param("an", an);
		@@ -621,24 +646,12 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
		instance->set_param("exril", exril);
		// Vectors and matrix need to be copied element-wise and
		// subsequently deleted.
		coord_vec = new double[nxv]();
		str_type = "FLOAT64_(" + to_string(nxv) + ")";
		status = hdf_file->read("XVEC", str_type, coord_vec);
		for (int xi = 0; xi < nxv; xi++)
		instance->set_x(xi, coord_vec[xi]);
		delete[] coord_vec;
		coord_vec = new double[nyv]();
		status = hdf_file->read("XVEC", str_type, instance->xv);
		str_type = "FLOAT64_(" + to_string(nyv) + ")";
		status = hdf_file->read("YVEC", str_type, coord_vec);
		for (int yi = 0; yi < nyv; yi++)
		instance->set_y(yi, coord_vec[yi]);
		delete[] coord_vec;
		coord_vec = new double[nzv]();
		status = hdf_file->read("YVEC", str_type, instance->yv);
		str_type = "FLOAT64_(" + to_string(nzv) + ")";
		status = hdf_file->read("ZVEC", str_type, coord_vec);
		for (int zi = 0; zi < nzv; zi++)
		instance->set_z(zi, coord_vec[zi]);
		delete[] coord_vec;
		status = hdf_file->read("ZVEC", str_type, instance->zv);
		_nlmmt = 2 * lm * (lm + 2);
		_nrvc = nxv * nyv * nzv;
		num_elements = 2 * _nlmmt * _nrvc;
		@@ -649,7 +662,7 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
		for (int wj = 0; wj < _nrvc; wj++) {
		int index = (2 * _nrvc) * wi + 2 * wj;
		value = complex<double>(elements[index], elements[index + 1]);
		instance->set_matrix_element(wi, wj, value);
		_wsum[wi][wj] = value;
		} // wj loop
		} // wi loop
		delete[] elements;
		@@ -662,6 +675,8 @@ TFRFME* TFRFME::from_hdf5(string file_name) {
		TFRFME* TFRFME::from_legacy(string file_name) {
		fstream input;
		TFRFME *instance = NULL;
		complex<double> **_wsum = NULL;
		double *coord_vec = NULL;
		input.open(file_name.c_str(), ios::in \| ios::binary);
		if (input.is_open()) {
		int lmode, lm, nkv, nxv, nyv, nzv;
		@@ -682,6 +697,7 @@ TFRFME* TFRFME::from_legacy(string file_name) {
		input.read(reinterpret_cast<char *>(&frsh), sizeof(double));
		input.read(reinterpret_cast<char *>(&exril), sizeof(double));
		instance = new TFRFME(lmode, lm, nkv, nxv, nyv, nzv);
		_wsum = instance->get_matrix();
		instance->set_param("vk", vk);
		instance->set_param("exri", exri);
		instance->set_param("an", an);
		@@ -690,17 +706,20 @@ TFRFME* TFRFME::from_legacy(string file_name) {
		instance->set_param("spd", spd);
		instance->set_param("frsh", frsh);
		instance->set_param("exril", exril);
		coord_vec = instance->get_x();
		for (int xi = 0; xi < nxv; xi++) {
		input.read(reinterpret_cast<char *>(&dval), sizeof(double));
		instance->set_x(xi, dval);
		coord_vec[xi] = dval;
		}
		coord_vec = instance->get_y();
		for (int yi = 0; yi < nyv; yi++) {
		input.read(reinterpret_cast<char *>(&dval), sizeof(double));
		instance->set_y(yi, dval);
		coord_vec[yi] = dval;
		}
		coord_vec = instance->get_z();
		for (int zi = 0; zi < nzv; zi++) {
		input.read(reinterpret_cast<char *>(&dval), sizeof(double));
		instance->set_z(zi, dval);
		coord_vec[zi] = dval;
		}
		int _nlmmt = 2 * lm * (lm + 2);
		int _nrvc = nxv * nyv * nzv;
		@@ -709,7 +728,7 @@ TFRFME* TFRFME::from_legacy(string file_name) {
		input.read(reinterpret_cast<char *>(&rval), sizeof(double));
		input.read(reinterpret_cast<char *>(&ival), sizeof(double));
		complex<double> value(rval, ival);
		instance->set_matrix_element(wi, wj, value);
		_wsum[wi][wj] = value;
		} // wj loop
		} // wi loop
		input.close();
		@@ -719,10 +738,6 @@ TFRFME* TFRFME::from_legacy(string file_name) {
		return instance;
		}

		std::complex<double> TFRFME::get_matrix_element(int row, int col) {
		return wsum[row][col];
		}

		long TFRFME::get_memory_requirement(
		int _lmode, int _lm, int _nkv, int _nxv,
		int _nyv, int _nzv
		@@ -758,10 +773,6 @@ double TFRFME::get_param(string param_name) {
		return value;
		}

		void TFRFME::set_matrix_element(int row, int col, complex<double> value) {
		wsum[row][col] = value;
		}

		void TFRFME::set_param(string param_name, double value) {
		if (param_name.compare("vk") == 0) vk = value;
		else if (param_name.compare("exri") == 0) exri = value;
		@@ -814,7 +825,7 @@ void TFRFME::write_hdf5(string file_name) {
		rec_ptr_list.append(&nzv);
		rec_name_list.append("VK");
		rec_type_list.append("FLOAT64_(1)");
		rec_ptr_list.append(&lm);
		rec_ptr_list.append(&vk);
		rec_name_list.append("EXRI");
		rec_type_list.append("FLOAT64_(1)");
		rec_ptr_list.append(&exri);

src/libnptm/tra_subs.cpp

+40 −45

Original line number	Diff line number	Diff line
		@@ -244,7 +244,7 @@ void ffrt(
		delete[] ctqcs;
		}

		void frfmer(
		complex<double> *frfmer(
		int nkv, double vkm, double vknmx, double apfafa, double tra,
		double spd, double rir, double ftcn, int le, int lmode, double pmf,
		Swap1 tt1, Swap2 tt2
		@@ -252,39 +252,32 @@ void frfmer(
		const int nlemt = le * (le + 2) * 2;
		const complex<double> cc0(0.0, 0.0);
		double sq = vkm * vkm;
		double *_vkv = tt2->get_vector();
		double **_vkzm = tt2->get_matrix();
		complex<double> *wk = new complex<double>[nlemt]();
		for (int jy90 = 0; jy90 < nkv; jy90++) {
		double vky = tt2->get_vector_element(jy90);
		double vky = _vkv[jy90];
		double sqy = vky * vky;
		for (int jx80 = 0; jx80 < nkv; jx80++) {
		double vkx = tt2->get_vector_element(jx80);
		double vkx = _vkv[jx80];
		double sqx = vkx * vkx;
		double sqn = sqx + sqy;
		double vkn = sqrt(sqn);
		if (vkn <= vknmx) {
		double vkz = sqrt(sq - sqn);
		wamff(tt1, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf);
		/*for (int j = 0; j < nlemt; j++) {
		double vreal = wk[j].real();
		double vimag = wk[j].imag();
		tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
		tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
		}*/
		wamff(wk, vkx, vky, vkz, le, apfafa, tra, spd, rir, ftcn, lmode, pmf);
		for (int j = 0; j < nlemt; j++) tt1->append(wk[j]);
		//tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
		tt2->set_matrix_element(jx80, jy90, vkz);
		_vkzm[jx80][jy90] = vkz;
		} else { // label 50
		for (int j = 0; j < nlemt; j++) {
		//double vreal = 0.0;
		//double vimag = 0.0;
		//tt1.write(reinterpret_cast<char *>(&vreal), sizeof(double));
		//tt1.write(reinterpret_cast<char *>(&vimag), sizeof(double));
		tt1->set_element(j, cc0);
		}
		for (int j = 0; j < nlemt; j++) tt1->append(cc0);
		//double vkz = 0.0;
		//tt2.write(reinterpret_cast<char *>(&vkz), sizeof(double));
		tt2->set_matrix_element(jx80, jy90, 0.0);
		_vkzm[jx80][jy90] = 0.0;
		}
		} // jx80 loop
		} // jy90 loop
		return wk;
		}

		void pwmalp(std::complex<double> *w, double up, double un, std::complex<double> ylm, int lw) {
		@@ -369,7 +362,7 @@ void sampoa(
		}

		void wamff(
		Swap1 *tt1, double x, double y, double z, int lm, double apfafa,
		complex<double> *wk, double x, double y, double z, int lm, double apfafa,
		double tra, double spd, double rir, double ftcn, int lmode, double pmf
		) {
		const int nlmm = lm * (lm + 2);
		@@ -389,27 +382,29 @@ void wamff(
		bool onx = (y == 0.0);
		bool ony = (x == 0.0);
		bool onz = (onx && ony);
		if (!(onz && onx && ony)) {
		if (!onz) {
		if (!onx) {
		if (!ony) {
		rho = sqrt(x * x + y * y);
		cph = x / rho;
		sph = y / rho;
		} else {
		if (onz) { // label 10
		cph = 1.0;
		sph = 0.0;
		} else {
		if (onx) { // label 12
		rho = sqrt(x * x);
		cph = (x < 0.0)? -1.0 : 1.0;
		sph = 0.0;
		} else {
		if (ony) { // label 13
		rho = sqrt(y * y);
		// goes to 15
		} else { // label 13
		rho = (y < 0.0) ? -y : y;
		cph = 0.0;
		sph = (y < 0.0) ? -1.0 : 1.0;
		// goes to 15
		}
		} else { // label 12
		rho = (x < 0.0) ? -x : x;
		cph = (x < 0.0) ? -1.0 : 1.0;
		sph = 0.0;
		// goes to 15
		}
		}
		} else { // label 10
		cph = 1.0;
		sph = 0.0;
		// goes to 15
		}
		// label 15
		if (z == 0.0) {
		@@ -464,13 +459,13 @@ void wamff(
		} else if (lmode == 2) { // label 50
		ftc1 = 2.0 * cthl / (cthl * rir + cth);
		cfam = (sth pmf * ftc1);
		for (int i52 = 0; i52 < nlmmt; i52++) tt1->set_element(i52, w[i52][0] * cfam);
		for (int i52 = 0; i52 < nlmmt; i52++) wk[i52] = w[i52][0] * cfam;
		// returns
		skip62 = true;
		} else if (lmode == 3) { // label 53
		ftc2 = 2.0 * cthl / (cthl + cth * rir);
		cfam = (sth pmf * ftc2);
		for (int i55 = 0; i55 < nlmmt; i55++) tt1->set_element(i55, w[i55][1] * cfam);
		for (int i55 = 0; i55 < nlmmt; i55++) wk[i55] = w[i55][1] * cfam;
		// returns
		skip62 = true;
		}
		@@ -485,7 +480,7 @@ void wamff(
		if (lmode == 0) {
		double f1 = cph * fam;
		double f2 = -sph * fam;
		for (int i58 = 0; i58 < nlmmt; i58++) tt1->set_element(i58, w[i58][0] * f1 + w[i58][1] * f2);
		for (int i58 = 0; i58 < nlmmt; i58++) wk[i58] = w[i58][0] * f1 + w[i58][1] * f2;
		// returns
		skip62 = true;
		} else if (lmode == 1) { // label 60
		@@ -496,19 +491,19 @@ void wamff(
		skip62 = false;
		} else if (lmode == 2) { // label 65
		fam = (pmf sth);
		for (int i67 = 0; i67 < nlmmt; i67++) tt1->set_element(i67, w[i67][0] * fam);
		for (int i67 = 0; i67 < nlmmt; i67++) wk[i67] = w[i67][0] * fam;
		// returns
		skip62 = true;
		} else if (lmode == 3) { // label 68
		fam = (pmf sth);
		for (int i70 = 0; i70 < nlmmt; i70++) tt1->set_element(i70, w[i70][1] * fam);
		for (int i70 = 0; i70 < nlmmt; i70++) wk[i70] = w[i70][1] * fam;
		// returns
		skip62 = true;
		}
		}
		if (!skip62) {
		if (lmode == 0 \|\| lmode == 1) { // label 62
		for (int i63 = 0; i63 < nlmmt; i63++) tt1->set_element(i63, w[i63][0] * cf1 + w[i63][1] * cf2);
		for (int i63 = 0; i63 < nlmmt; i63++) wk[i63] = w[i63][0] * cf1 + w[i63][1] * cf2;
		}
		}
		}

src/trapping/cfrfme.cpp

+36 −25

File changed.

Preview size limit exceeded, changes collapsed.