Commit 73f64cc5 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Add cluster specific funtions GHIT, CGEV and R3JRR

parent dea86173

src/include/clu_subs.h

+121 −5
Original line number Diff line number Diff line
@@ -61,7 +61,7 @@ double cgev(int ipamo, int mu, int l, int m) { 
		}

	} else { // label 30

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

		if (mu < 0) { // label 35. CHECK: is clu.f code line 2466 a switch?

		if (mu < 0) { // label 35. CHECK: is clu.f code line 2458 a switch?

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

		} else if (mu == 0) { // label 40

			xn = 2.0 * (l - m) * (l + m);
@@ -82,7 +82,123 @@ double cgev(int ipamo, int mu, int l, int m) { 
 * \param c6: `C6 *`

 */

void r3jjr(int j2, int j3, int m2, int m3, C6 *c6) {



	int jmx = j3 + j2;

	int jdf = j3 - j2;

	int m1 = m2 - m3;

	int abs_jdf = (jdf >= 0) ? jdf : -jdf;

	int abs_m1 = (m1 >= 0) ? m1 : -m1;

	int jmn = (abs_jdf > abs_m1) ? abs_jdf : abs_m1;

	int njmo = jmx - jmn;

	int jf = jmx + jmx + 1;

	int isn = 1;

	if ((jdf + m1) % 2 != 0) isn = -1;

	if (njmo <= 0) {

		double sj = 1.0 * jf;

		double cnr = (1.0 / sqrt(sj)) * isn;

		c6->rac3j[0] = cnr;

		return;

	} else { // label 15

		double sjt = 1.0;

		double sjr = 1.0 * jf;

		int jsmpos = (jmx + 1) * (jmx + 1);

		int jdfs = jdf * jdf;

		int m1s = m1 * m1;

		int mdf = m3 - m2;

		int idjc = m1 * (j3 * (j3 + 1) - j2 * (j2 +1));

		int j1 = jmx;

		int j1s = j1 * j1;

		int j1po = j1 + 1;

		double ccj = 1.0 * (j1s - jdfs) * (j1s - m1s);

		double cj = sqrt(ccj * (jsmpos - j1s));

		double dj = 1.0 * jf * (j1 * j1po * mdf + idjc);

		double cjp = 0.0; // WARNING: not sure it should be really defined here.

		if (njmo <= 1) {

			c6->rac3j[0] = -dj /(cj * j1po);

			double sj = sjr + (c6->rac3j[0] * c6->rac3j[0]) * (jf - 2);

			double cnr = (1.0 / sqrt(sj)) * isn;

			c6->rac3j[1] = cnr;

			c6->rac3j[0] *= cnr;

			return;

		} else { // label 20

			int nj = njmo + 1;

			int nmat = (nj + 1) / 2;

			c6->rac3j[nj - 1] = 1.0;

			c6->rac3j[njmo - 1] = -dj / (cj * j1po);

			if (nmat != njmo) {

				int nbr = njmo - nmat;

				for (int ibr45 = 1; ibr45 <= nbr; ibr45++) {

					int irr = nj - ibr45;

					jf -= 2;

					j1--;

					j1s = j1 * j1;

					j1po = j1 + 1;

					cjp = cj;

					ccj = 1.0 * (j1s - jdfs) * (j1s - m1s);

					cj = sqrt(ccj * (jsmpos - j1s));

					sjt = c6->rac3j[irr - 1] * c6->rac3j[irr - 1];

					dj = 1.0 * jf * (j1 * j1po * mdf + idjc);

					c6->rac3j[irr - 2] = -c6->rac3j[irr - 1] * dj

							+ c6->rac3j[irr] * cjp * j1 / (cj * j1po);

					sjr += (sjt * jf);

				} // ibr45 loop

			}

			// label 50

			double osjt = sjt;

			sjt = c6->rac3j[nmat - 1] * c6->rac3j[nmat - 1];

			if (sjt >= osjt) {

				sjr += (sjt * (jf - 2));

			} else { // label 55

				nmat++;

			}

			// label 60

			double racmat = c6->rac3j[nmat - 1];

			c6->rac3j[0] = 1.0;

			jf = jmn + jmn + 1;

			double sjl = 1.0 * jf;

			j1 = jmn;

			if (j1 != 0) {

				j1po = j1 + 1;

				int j1pos = j1po * j1po;

				double ccjp = 1.0 * (j1pos - jdfs) * (j1pos - m1s);

				cjp = sqrt(ccjp * (jsmpos - j1pos));

				dj = 1.0 * jf * (j1 * j1po * mdf + idjc);

				c6->rac3j[1] = - dj / (cjp * j1);

			} else { // label 62

				cjp = sqrt(1.0 * (jsmpos - 1));

				dj = 1.0 * mdf;

				c6->rac3j[1] = -dj / cjp;

			}

			// label 63

			int nmatmo = nmat - 1;

			if (nmatmo >= 2) {

				for (int irl70 = 2; irl70 <= nmatmo; irl70++) {

					jf += 2;

					j1++;

					j1po = j1 + 1;

					int j1pos = j1po * j1po;

					cj = cjp;

					double ccjp = 1.0 * (j1pos - jdfs) * (j1pos - m1s);

					cjp = sqrt(ccjp * (jsmpos - j1pos));

					sjt = c6->rac3j[irl70 - 1] * c6->rac3j[irl70 - 1];

					dj = 1.0 * jf * (j1 * j1po * mdf + idjc);

					c6->rac3j[irl70] = -(

							c6->rac3j[irl70 - 1] * dj

							+ c6->rac3j[irl70 - 2] * cj * j1po

					) / (cjp * j1);

					sjl += (sjt * jf);

				}

			}

			// label 75

			double ratrac = racmat / c6->rac3j[nmat - 1];

			double rats = ratrac * ratrac;

			double sj = sjr + sjl * rats;

			c6->rac3j[nmat - 1] = racmat;

			double cnr = (1.0 / sqrt(sj)) * isn;

			for (int irr80 = nmat; irr80 <= nj; irr80++) c6->rac3j[irr80 - 1] *= cnr;

			double cnl = cnr * ratrac;

			for (int irl85 = 1; irl85 <= nmatmo; irl85++) c6->rac3j[irl85 - 1] *= cnl;

		}

	}

}
@@ -167,7 +283,7 @@ std::complex<double> ghit( 
						lt14 += 2;

					}

				} else { // label 16

					// Would call R3JJR(L1MP, L2, -MUPM1, MUPM2)

					r3jjr(l1mp, l2, -mupm1, mupm2, c6);

					int il = ilin;

					int lt20 = lminpo;

					while (lt20 <= lmaxpo) {
@@ -221,7 +337,7 @@ std::complex<double> ghit( 
							lt34 += 2;

						}

					} else { // label 36

						// Would call R3JJR

						r3jjr(l1mp, l2, -mupm1, mupm2, c6);

						int il = ilin;

						int lt40 = lminpo;

						while (lt40 <= lmaxpo) {
@@ -280,7 +396,7 @@ std::complex<double> ghit( 
						lt54 += 2;

					}

				} else { // label 56

					// Would call R3JJR(L1MP, L2, -MUPM1, MUPM2)

					r3jjr(l1mp, l2, -mupm1, mupm2, c6);

					int il = ilin;

					int lt60 = lminpo;

					while (lt60 <= lmaxpo) {