C implementation of the B translation coefficient.

Former-commit-id: f7c8930de88c1d62ac99e29e9adda4ca78aa216e

qpms/BUGS.rst

@@ -0,0 +1,6 @@
+gaunt.c
+=======
+abort při určitých vstupech, např. -6 10 8 10
+(ačkoliv fortran originál paří bez problémů)
+
+

qpms/test_translations.c

@@ -16,9 +16,19 @@ testcase_single_trans_t testcases_Taylor[] = {
 
 int main() {
 	for(testcase_single_trans_t *tc = testcases_Taylor; tc->J != QPMS_BESSEL_UNDEF; tc++) {
-	complex double A = qpms_trans_single_A_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J);
+		if (tc->n > 12 || tc->nu > 12 || !tc->n || !tc->nu ) continue;
-	printf("m=%d, n=%d, mu=%d, nu=%d, relerr=%.16f\n", tc->m,tc->n,tc->mu,tc->nu,
+
-			cabs(tc->result_A - A)/((cabs(A) < cabs(tc->result_A)) ? cabs(A) : cabs(tc->result_A)));
+		printf("m=%d, n=%d, mu=%d, nu=%d,\n", tc->m,tc->n,tc->mu,tc->nu);
+		complex double A = qpms_trans_single_A_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J);
+		complex double B = qpms_trans_single_B_Taylor(tc->m, tc->n, tc->mu, tc->nu, tc->kdlj, true, tc->J);
+		printf("A = %.16f+%.16fj, relerr=%.16f, J=%d\n",
+			creal(A), cimag(A),
+			cabs(tc->result_A - A)/((cabs(A) < cabs(tc->result_A)) ? cabs(A) : cabs(tc->result_A)),
+			tc->J);
+		printf("B = %.16f+%.16fj, relerr=%.16f, J=%d\n",
+			creal(B), cimag(B),
+			cabs(tc->result_B - B)/((cabs(B) < cabs(tc->result_B)) ? cabs(B) : cabs(tc->result_B)),
+			tc->J);
 	}
 }
 

qpms/translations.c

@@ -65,10 +65,9 @@ int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double
 complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kdlj,
 		bool r_ge_d, qpms_bessel_t J) { // TODO make J enum
     if(r_ge_d) J = QPMS_BESSEL_REGULAR;
-    double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
+
-                +lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1)
-                +lgamma(n+nu+1) - lgamma(2*(n+nu)+1));
     double costheta = cos(kdlj.theta);
+
     int qmax = gaunt_q_max(-m,n,mu,nu); // nemá tu být +m?
     // N.B. -m !!!!!!
     double a1q[qmax+1];
@@ -76,8 +75,7 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd
     gaunt_xu(-m,n,mu,nu,qmax,a1q,&err);
     double a1q0 = a1q[0];
     if (err) abort();
-    //double *leg = malloc(sizeof(double)*gsl_sf_legendre_array_n(n+nu));
+
-    //if (!leg) abort();
     double leg[gsl_sf_legendre_array_n(n+nu)];
     if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort();
     complex double bes[n+nu+1];
@@ -94,8 +92,10 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd
 	    complex double summandq = (n*(n+1) + nu*(nu+1) - p*(p+1)) * min1pow(q) * a1q_n * zp * Pp;
 	    sum += summandq;
     }
-    //free(leg);
+    
-    //free(bes);
+    double exponent=(lgamma(2*n+1)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
+                +lgamma(n+nu+m-mu+1)-lgamma(n-m+1)-lgamma(nu+mu+1)
+                +lgamma(n+nu+1) - lgamma(2*(n+nu)+1));
     complex double presum = exp(exponent);
     presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n) / (4*n);
  
@@ -104,3 +104,61 @@ complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kd
     return (presum / prenormratio) * sum;
 }
 
+complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kdlj,
+		bool r_ge_d, qpms_bessel_t J) { // TODO make J enum
+    if(r_ge_d) J = QPMS_BESSEL_REGULAR;
+    double costheta = cos(kdlj.theta);
+    
+    int q2max = gaunt_q_max(-m-1,n+1,mu+1,nu);
+    int Qmax = gaunt_q_max(-m,n+1,mu,nu);
+    double a2q[q2max+1], a3q[Qmax+1], a2q0, a3q0;
+    int err;
+    if (mu == nu) {
+	    for (int q = 0; q <= q2max; ++q) 
+	    	a2q[q] = 0;
+	    a2q0 = 1;
+    }
+    else { 
+	    gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort();
+	    a2q0 = a2q[0];
+    }
+    gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err); if (err) abort();
+    a3q0 = a3q[0];
+
+    double leg[gsl_sf_legendre_array_n(n+nu+1)];
+    if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort();
+    complex double bes[n+nu+2];
+    if (qpms_sph_bessel_array(J, n+nu+1, kdlj.r, bes)) abort();
+ 
+    complex double sum = 0;
+    for (int q = 0; q <= Qmax; ++q) {
+	    int p = n+nu-2*q;
+	    double a2q_n = a2q[q]/a2q0;
+	    double a3q_n = a3q[q]/a3q0;
+	    complex double zp_ = bes[p+1];
+	    int Pp_order_ = mu-m;
+	    if(p+1 < abs(Pp_order_)) continue; // FIXME raději nastav lépe meze
+	    double Pp_ = leg[gsl_sf_legendre_array_index(p+1, abs(Pp_order_))];
+	    if (Pp_order_ < 0) Pp_ *= min1pow(mu-m) * exp(lgamma(1+1+p+Pp_order_)-lgamma(1+1+p-Pp_order_));
+	    complex double summandq = ((2*(n+1)*(nu-mu)*a2q_n
+                 -(-nu*(nu+1) - n*(n+3) - 2*mu*(n+1)+p*(p+3))* a3q_n)
+                *min1pow(q) * zp_ * Pp_);
+	    sum += summandq;
+    }
+
+    double exponent=(lgamma(2*n+3)-lgamma(n+2)+lgamma(2*nu+3)-lgamma(nu+2)
+                +lgamma(n+nu+m-mu+2)-lgamma(n-m+1)-lgamma(nu+mu+1)
+                +lgamma(n+nu+2) - lgamma(2*(n+nu)+3));
+    complex double presum = exp(exponent);
+    presum *= cexp(I*(mu-m)*kdlj.phi) * min1pow(m) * ipow(nu+n+1) / (
+        (4*n)*(n+1)*(n+m+1));
+
+    // Taylor normalisation v2, proven to be equivalent
+    complex double prenormratio = ipow(nu-n) * sqrt(((2.*nu+1)/(2.*n+1))* exp(
+        lgamma(n+m+1)-lgamma(n-m+1)+lgamma(nu-mu+1)-lgamma(nu+mu+1)));
+    
+    return (presum / prenormratio) * sum;
+}
+
+
+

qpms/translations.h

@@ -20,5 +20,7 @@ int qpms_sph_bessel_array(qpms_bessel_t typ, int lmax, double x, complex double
 complex double qpms_trans_single_A_Taylor(int m, int n, int mu, int nu, sph_t kdlj,
                 bool r_ge_d, qpms_bessel_t J);
 
+complex double qpms_trans_single_B_Taylor(int m, int n, int mu, int nu, sph_t kdlj,
+                bool r_ge_d, qpms_bessel_t J);
 
 #endif // QPMS_TRANSLATIONS_H