gaunt pokračování přepisu

Former-commit-id: aed3bfd9a4652e79e1619c22b873c82ea3ed17ea

qpms/gaunt.c

@@ -2,6 +2,10 @@
 #include <math.h>
 #include <stdio.h>
 
+#define ZERO_THRESHOLD 1.e-8
+
+//  "Besides, the determined Real Programmer can write FORTRAN programs in any language."
+//  -- Ed Post, Real Programmers Don't Use Pascal, 1982.
 
 // logarithm of factorial (from basicsubs.f90)
 double lnf (double z) {
@@ -37,9 +41,93 @@ double f_a0 (int m, int n, int mu, int nu) {
 	return exp(logw+logp);
 }
 
+// coefficient Ap(m,n,mu,nu,p) (from vec_trans.f90)
+int Ap(int m, int n, int mu, int nu, int p) {
+	return p*(p-1)*(m-mu)-(m+mu)*(n-nu)*(n+nu+1);
+}
+
+// coefficient a(m,n,mu,nu,1) normalized by the backward recursion (from vec_trans.f90)
+
+double f_a1norm(int m, int n, int mu, int nu) {
+	int n4 = n + nu - m - mu;
+	return ((2.*n + 2.*nu - 3.) / 2.) * (1. - ((2.*n + 2.*nu - 1.) / (n4 * (n4-1.))) 
+			* ((m - n) * (m - n + 1.) / (2.*n - 1.) + (mu-nu) * (mu-nu+1.)/(2.nu-1.)));
+}
+
+// coefficient a(m,n,mu,nu,2) normalized by the backward recursion (From vec_trans.f90)
+
+double f_a2norm(int m, int n, int mu, int nu) {
+	double n4 = n + nu - m - mu;
+	double n2nu = 2*n + 2*nu;
+	double mn = m - n;
+	double munu = mu - nu;
+
+	return ((n2nu-1.)*(n2nu-7.)/4.) * \
+                 ( ((n2nu-3.)/(n4*(n4-1.))) * \
+                 ( ((n2nu-5.)/(2.*(n4-2.)*(n4-3.))) * \
+                 ( mn*(mn+1.)*(mn+2.)*(mn+3.)/((2.*n-1.)*(2.*n-3.)) + \
+                   2.*mn*(mn+1.)*munu*(munu+1.)/((2.*n-1.)*(2.*nu-1.)) + \
+                   munu*(munu+1.)*(munu+2.)*(munu+3.)/((2.*nu-1.)*(2.*nu-3.)) \
+                   ) - mn*(mn+1.)/(2.*n-1.) - munu*(munu+1.)/(2.*nu-1.) ) +0.5);
+}
+
+// just for convenience – square of ints
+int isq(int x) {return x * x;}
+
+double f_alpha(int n, int nu, int p) {
+	return (isq(p) - isq(n-nu))*(isq(p)-isq(n+nu+1)) / (double)(4*isq(p)-1);
+}
+
+// coeff a(m,n,mu,nu,2) normalizzato per la backward recursion calcolato questa volta per ricorsione
+// (from vec_trans.f90)
+double f_a2normr(int m, int n, int mu, int nu, double a1norm) {
+	int p = n + nu - 4;
+	int p1 = p - m - mu;
+	int p2 = p + m + mu;
+	int Ap4 = f_Ap(m,n,mu,nu,p+4);
+	int Ap6 = f_Ap(m,n,mu,nu,p+6);
+	
+	double alphap1=f_alpha(n,nu,p+1);
+	double alphap2=f_alpha(n,nu,p+2);
+	
+	if (!Ap4) {
+		if(!Ap6) {
+	
+	                double c0=(p+2)*(p1+1)*alphap1;
+	                double c1=(p+1)*(p2+2)*alphap2;
+	
+	                return (c1/c0)*a1norm;
+		} else /* Ap6 != 0 */ {
+	                int Ap2=f_Ap(m,n,mu,nu,p+2);
+	                int Ap3=f_Ap(m,n,mu,nu,p+3);
+	                int Ap5=f_Ap(m,n,mu,nu,p+5);
+	                double alphap5=f_alpha(n,nu,p+5);
+	
+	                double c0=(p+2)*(p+3)*(p+5)*(p1+1)*(p1+2)*(p1+4)*Ap6*alphap1;
+	                double c1=(p+5)*(p1+4)*Ap6*(Ap2*Ap3+(p+1)*(p+3)*(p1+2)*(p2+2)*alphap2);
+	                double c2=(p+2)*(p2+3)*Ap2*(Ap5*Ap6+(p+4)*(p+6)*(p1+5)*(p2+5)*alphap5);
+	
+	                return (c1/c0)*a1norm+(c2/c0);
+		}
+	} else /* Ap4 != 0 */ {
+	                int Ap2=f_Ap(m,n,mu,nu,p+2);
+	                int Ap3=f_Ap(m,n,mu,nu,p+3);
+	                double alphap3=f_alpha(n,nu,p+3);
+	                double alphap4=f_alpha(n,nu,p+4);
+	
+	                double c0=(p+2)*(p+3)*(p1+1)*(p1+2)*Ap4*alphap1;
+	                double c1=Ap2*Ap3*Ap4+(p+1)*(p+3)*(p1+2)*(p2+2)*Ap4*alphap2 \
+	                 + (p+2)*(p+4)*(p1+3)*(p2+3)*Ap2*alphap3;
+	                double c2=-(p+2)*(p+3)*(p2+3)*(p2+4)*Ap2*alphap4;
+	                
+			return (c1/c0)*a1norm+(c2/c0);
+	}
+}
+
+// gaunt_xu from vec_trans.f90
+// btw, what is the difference from gaunt_xu2?
+double gaunt_xu2(int m, int n, int mu, int nu, int qmax, double *v_aq, int *error) {
 	
-/*
-double gaunt_gevero_direct_convert(int m, int n, int mu, int nu, int qmax, double *v_aq, int *error) {
 	int v_zero[qmax] = {0};
 	*error = 0;
 
@@ -55,8 +143,108 @@ double gaunt_gevero_direct_convert(int m, int n, int mu, int nu, int qmax, doubl
 			break;
 		case 1:
 			v_aq[0] = f_a0(m,n,mu,nu);
+			v_aq[1] = v_aq[0] + f_a1norm(m,n,mu,nu);
+			break;
+		case 2:
+			v_aq[0] = f_a0(m,n,mu,nu);
+			v_aq[1] = v_aq[0] + f_a1norm(m,n,mu,nu);
+			v_aq[2] = v_aq[0] * f_a2normr(m,n,mu,nu,v_aq[1]/v_aq[0]);
+			break;
+		default:
+			if (m == 0 && mu == 0) {
+				v_aq[0] = f_a0(m,n,mu,nu);
+
+				// backward recursion
+				for (int q = 1; q <= qmax; ++q) {
+					int p = n + nu - 2*q;
+					double c0 = f_alpha(n,nu,p+1);
+					double c1 = f_alpha(n,nu,p+2);
+					v_aq[q] = (c1/c0) * v_aq[q-1];
+				}
+			else if (mu == m && nu == n) {
+				v_aq[0] = f_a0(m,n,mu,nu);
+
+				// backward recursion
+				for (int q = 1; q <= qmax; ++q) {
+					// calculate pre-coefficients
+					int p = n + nu - 2*q;
+					int p1 = p - m - mu;
+					int p2 = p + m + mu;
+
+					// calculate recursion coefficients
+					double c0 = (p+2) * (p1+1) * f_alpha(n,nu,p+1);
+					double c1 = (p+1) * (p2+2) * f_alpha(n,nu,p+2);
+
+					//recursion
+					v_aq[q] = (c1/c0) * v_aq[q-1];
+				}
+			} else if (mu == -m) {
+				v_aq[0] = f_a0(m,n,mu,nu);
+				v_aq[1] = f_a1norm(m,n,mu,nu)*v_aq[0];
+
+				// backward recursion
+				for (int q = 2; q <= qmaq; ++q) {
+					// calculate pre-coefficient
+					int p = n + nu - 2*q;
+
+					// calculate recursion coefficients
+					double c0 = f_alpha(n, nu, p+1);
+					double c1 = 4*isq(m) + f_alpha(n,nu,p+2) + f_alpha(n,nu,p+3);
+					double c2 = - f_alpha(n, nu, p+4);
+
+					// recursion
+					v_aq[q] = (c1/c0)*v_aq[q-1] + (c2/c0)*v_aq[q-2];
+				}
+			}  else if  // TODO vec_trans.f90:1233
+
+		
+
 			!!!!!!!!!TODO CONTINUE HERE
 
 
+/*
+// gaunt_xu2 from vec_trans.f90
+// btw, what is the difference from gaunt_xu?
+double gaunt_xu2(int m, int n, int mu, int nu, int qmax, double *v_aq, int *error) {
 	
+	int v_zero[qmax] = {0};
+	*error = 0;
+
+	if(abs(m)>n || abs(mu)=nu) {
+		*error = 1;
+		fprintf(stderr, "invalid values for m, n, mu or nu\n")
+		return NAN;
+	}
+
+	switch(qmax) {
+		case 0:
+			v_aq[0]  = f_a0(m,n,mu,nu);
+			break;
+		case 1:
+			v_aq[0] = f_a0(m,n,mu,nu);
+			v_aq[1] = v_aq[0] + f_a1norm(m,n,mu,nu);
+
+			// Check for zeros
+			if (fabs(v_aq[1]/v_aq[0]) < ZERO_THRESHOLD) {
+				v_aq[1] = 0.;
+				v_zero[1] = 0.;
+			}
+			break;
+		case 2:
+			v_aq[0] = f_a0(m,n,mu,nu);
+			v_aq[1] = v_aq[0] + f_a1norm(m,n,mu,nu);
+	
+			// Check for zeros
+			if (fabs(v_aq[1]/v_aq[0]) < ZERO_THRESHOLD) {
+				v_aq[1] = 0.;
+				v_zero[1] = 0.;
+			}
+
+			v_aq[2] = v_aq[0] * f_a2norm_
+			break;
+		
+
+			!!!!!!!!!TODO CONTINUE HERE
 */
+
+

qpms/qpms_c.pyx

@@ -66,3 +66,9 @@ def get_y_mn_unsigned(int nmax):
             ymn_plus[m,n] = i
             i = i + 1
     return(ymn_plus, ymn_minus)
+
+cdef int q_max(int m, int n, int mu, int nu):
+    return min(n,nu,(n+nu-abs(m+mu)//2)
+
+
+#cdef translation_coefficient_A(