Commit 94f58128 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Introduce FORTRAN/C++ correspondence logs

parent a531bd6b

src/include/sph_subs.h

+133 −0
Original line number Diff line number Diff line
@@ -49,6 +49,139 @@ struct C2 { 
	double vsz[2];

};



/*! \brief Conjugate of a double precision complex number

 *

 * \param z: `std::complex\<double\>` The input complex number.

 * \return result: `std::complex\<double\>` The conjugate of the input number.

 */

std::complex<double> dconjg(std::complex<double> z) {

	double zreal = z.real();

	double zimag = z.imag();

	return std::complex<double>(zreal, -zimag);

}



/*! \brief C++ porting of CG1

 *

 * \param li: `int`

 * \param mu: `int`

 * \param l: `int`

 * \param m: `int`

 * \return result: `double`

 */

double cg1(int lmpml, int mu, int l, int m) {

	double result = 0.0;

	double xd, xn;

	if (lmpml == -1) { // Interpreted as GOTO 30

		xd = 2.0 * l * (2 * l - 1);

		if (mu == -1) {

			xn = 1.0 * (l - 1 - m) * (l - m);

		} else if (mu == 0) {

			xn = 2.0 * (l - m) * (l + m);

		} else if (mu == 1) {

			xn = 1.0 * (l - 1 + m) * (l + m);

		} else {

			throw 111; // Need an exception for unpredicted indices.

		}

		result = sqrt(xn / xd);

	} else if (lmpml == 0) { // Interpreted as GOTO 5

		bool goto10 = (m != 0) || (mu != 0);

		if (!goto10) {

			result = 0.0;

			return result;

		}

		if (mu != 0) {

			xd = 2.0 * l * (l + 1);

			if (mu == -1) {

				xn = 1.0 * (l - m) * (l + m + 1);

				result = -sqrt(xn / xd);

			} else if (mu == 1) { // mu > 0

				xn = 1.0 * (l + m) * (l - m + 1);

				result = sqrt(xn / xd);

			} else {

				throw 111; // Need an exception for unpredicted indices.

			}

		} else { // mu = 0

			xd = 1.0 * l * (l + 1);

			xn = -1.0 * m;

			result = xn / sqrt(xd);

		}

	} else if (lmpml == 1) { // Interpreted as GOTO 60

		xd = 2.0 * (l * 2 + 3) * (l + 1);

		if (mu == -1) {

			xn = 1.0 * (l + 1 + m) * (l + 2 + m);

			result = sqrt(xn / xd);

		} else if (mu == 0) {

			xn = 2.0 * (l + 1 - m) * (l + 1 + m);

			result = -sqrt(xn / xd);

		} else if (mu == 1) {

			xn = 1.0 * (l + 1 - m) * (l + 2 - m);

			result = sqrt(xn / xd);

		} else { // mu was not recognized.

			throw 111; // Need an exception for unpredicted indices.

		}

	} else { // lmpml was not recognized

		throw 111; // Need an exception for unpredicted indices.

	}

	return result;

}



/*! \brief C++ porting of APS

 *

 * \param zpv: `double ****`

 * \param li: `int`

 * \param nsph: `int`

 * \param c1: `C1 *`

 * \param sqk: `double`

 * \param gaps: `double *`

 */

void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {

	std::complex<double> cc0 = std::complex<double>(0.0, 0.0);

	std::complex<double> summ, sume, suem, suee, sum;

	double half_pi = acos(0.0);

	double cofs = half_pi * 2.0 / sqk;

	for (int i40 = 1; i40 <= nsph; i40++) {

		int iogi = c1->iog[i40 - 1];

		if (iogi >= i40) {

			sum = cc0;

			for (int l30 = 1; l30 <= li; l30++) {

				int ltpo = l30 + l30 + 1;

				for (int ilmp = 1; ilmp <= 3; ilmp++) {

					bool goto30 = (l30 == 1 && ilmp == 1) || (l30 == li && ilmp == 3);

					if (!goto30) {

						int lmpml = ilmp - 2;

						int lmp = l30 + lmpml;

						double cofl = sqrt(1.0 * (ltpo * (lmp + lmp + 1)));

						summ = zpv[l30 - 1][ilmp - 1][0][0] /

								(

										dconjg(c1->rmi[l30 - 1][i40 - 1]) *

										c1->rmi[lmp - 1][i40 - 1]

								);

						sume = zpv[l30 - 1][ilmp - 1][0][1] /

								(

										dconjg(c1->rmi[l30 - 1][i40 - 1]) *

										c1->rei[lmp - 1][i40 - 1]

								);

						suem = zpv[l30 - 1][ilmp - 1][1][0] /

								(

										dconjg(c1->rei[l30 - 1][i40 - 1]) *

										c1->rmi[lmp - 1][i40 - 1]

								);

						suee = zpv[l30 - 1][ilmp - 1][1][1] /

								(

										dconjg(c1->rei[l30 - 1][i40 - 1]) *

										c1->rei[lmp - 1][i40 - 1]

								);

						sum += (cg1(lmpml, 0, l30, -1) * (summ - sume - suem + suee) +

								cg1(lmpml, 0, l30, 1) * (summ + sume + suem + suee)) * cofl;

					}

				}

			}

		}

		gaps[i40 - 1] = sum.real() * cofs;

		printf("DEBUG: gaps[%d] = %lE\n", i40, gaps[i40 - 1]);

	}

}



/*! \brief C++ porting of DIEL

 *

 * \param npntmo: `int`

src/sphere/sph.f

+1 −0
Original line number Diff line number Diff line
@@ -645,6 +645,7 @@ CCC 1RMI(LI,NSPH),REI(LI,NSPH),GAPS(NSPH) 
     1CG1(LMPML,0,L,1)*(SUMM+SUME+SUEM+SUEE))*COFL

   30 CONTINUE

      GAPS(I)=DREAL(SUM)*COFS

      PRINT *,"DEBUG: GAPS(",I,") =",GAPS(I)

   40 CONTINUE

      RETURN

      END

src/sphere/sphere.cpp

+4 −3
Original line number Diff line number Diff line
@@ -9,9 +9,8 @@ using namespace std; 


//! \brief C++ implementation of SPH

void sphere() {

	//complex<double> dc0[5];

	//complex<double> dc0m[2][4];

	complex<double> arg, s0, tfsas;

	double gaps[2];

	printf("INFO: making legacy configuration ...\n");

	ScattererConfiguration *conf = ScattererConfiguration::from_dedfb("../../test_data/sphere/DEDFB");

	conf->write_formatted("c_OEDFB");
@@ -204,7 +203,6 @@ void sphere() { 
					int iogi = sconf->iog_vec[i132 - 1];

					if (iogi != i132) {

						for (int l123 = 1; l123 <= gconf->l_max; l123++) {

							// QUESTION: aren't RMI and REI still empty?

							c1->rmi[l123 - 1][i132 - 1] = c1->rmi[l123 - 1][iogi - 1];

							c1->rei[l123 - 1][i132 - 1] = c1->rei[l123 - 1][iogi - 1];

						}
@@ -265,6 +263,9 @@ void sphere() { 
				double cs0 = 0.25 * vk * vk * vk / half_pi;

				sscr0(tfsas, nsph, gconf->l_max, vk, exri, c1);

				printf("DEBUG: TFSAS = (%lE,%lE)\n", tfsas.real(), tfsas.imag());

				double sqk = vk * vk * sconf->exdc;

				aps(zpv, gconf->l_max, nsph, c1, sqk, gaps);

				// Would call RABAS

			} //jxi488

			tppoan.close();

		} else { // In case TPPOAN could not be opened. Should never happen.