From e98a4656ca9b958c38ceaa54399c2c14e2d0796d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 17 Mar 2017 19:58:50 +0200 Subject: [PATCH] =?UTF-8?q?gaunt=20pokra=C4=8Dov=C3=A1n=C3=AD=20p=C5=99epi?= =?UTF-8?q?su?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Former-commit-id: aed3bfd9a4652e79e1619c22b873c82ea3ed17ea --- qpms/gaunt.c | 192 +++++++++++++++++++++++++++++++++++++++++++++++- qpms/qpms_c.pyx | 6 ++ 2 files changed, 196 insertions(+), 2 deletions(-) diff --git a/qpms/gaunt.c b/qpms/gaunt.c index 449c690..551827a 100644 --- a/qpms/gaunt.c +++ b/qpms/gaunt.c @@ -2,6 +2,10 @@ #include #include +#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); +} -/* -double gaunt_gevero_direct_convert(int m, int n, int mu, int nu, int qmax, double *v_aq, int *error) { +// 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) { + 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 */ + + diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 1a1771b..e9c5e66 100644 --- a/qpms/qpms_c.pyx +++ b/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(