From 73f64cc5db9234239787d301358d895eb4a31ecc Mon Sep 17 00:00:00 2001
From: Giovanni La Mura <giovanni.lamura@inaf.it>
Date: Mon, 13 Nov 2023 14:12:44 +0100
Subject: [PATCH] Add cluster specific funtions GHIT, CGEV and R3JRR
---
src/include/clu_subs.h | 126 +++++++++++++++++++++++++++++++++++++++--
1 file changed, 121 insertions(+), 5 deletions(-)
diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 52568582..a0bfc929 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -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) {
--
GitLab