diff --git a/src/Makefile.am b/src/Makefile.am index 8657170..f91b0a0 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -34,7 +34,6 @@ libconviqt_la_SOURCES = \ conviqt_math.cpp \ conviqt_io.cpp \ conviqt_util.cpp \ - wignercalc.cpp \ cconviqt.cpp libconviqt_la_LIBADD = \ diff --git a/src/conviqt.hpp b/src/conviqt.hpp index c2faca1..b240847 100644 --- a/src/conviqt.hpp +++ b/src/conviqt.hpp @@ -268,15 +268,6 @@ private : void timing_line(std::string label, double timer, long counter); }; -// wignercalc.cpp - -double xpow(int expo, double val); -void wignerCalc(tsize n, tsize mmax, double theta, levels::arr2 &d); -void wignerCalcGeneral(tsize n, tsize mmax, double theta, levels::arr2 &d); -void wignerCalcHalfpi(tsize n, tsize mmax, levels::arr2 &d); -double wignerCalc00_halfpi(tsize n); -double wignerCalc00(double theta, tsize n); - // conviqt_util.cpp void sift_down_DL(levels::arr &ra, levels::arr &brr, const int l, const int r); diff --git a/src/libconviqt_io.cpp b/src/libconviqt_io.cpp deleted file mode 100644 index e69de29..0000000 diff --git a/src/wignercalc.cpp b/src/wignercalc.cpp deleted file mode 100644 index 57378df..0000000 --- a/src/wignercalc.cpp +++ /dev/null @@ -1,275 +0,0 @@ -#include "conviqt.hpp" - -namespace conviqt { - -double wigner_xpow (int expo, double val) - { return (expo&1) ? -val : val; } - -double wignerCalc00_halfpi( tsize n ) -{ - if (n&1) return 0.; - double rec=1.; - for (tsize k=1; k &d ) -{ - - try { - d.fast_alloc(2*n+1,2*mmax+1); - } catch ( std::bad_alloc & e ) { - std::cerr << " Out of memory allocating " << (2*n+1)*(2*mmax+1)*8./1024/1024 << "MB for wignerCalcGeneral" << std::endl; - throw; - } - - double cth=cos(theta), sth=sin(theta); - double xsth=1./sth; - double fact1 = std::sqrt(2.*n)*cos(theta/2.)*sin(theta/2.); - levels::arr prod1(n+1), prod2(n+1); - - levels::arr sqtab(n); - for (tsize i=0; i0) - { - d00 = d[n][mmax+1]; - d00_2 = d[n][mmax-1]; - } - for (tsize i=1; i<=mmax; ++i) - { - //I STILL NEED TO WORRY ABOUT THE PROD1 ARRAY HAVING VANISHING DENOMINATORS - bool flag1=false, flag2=false; - if (abs(d00) < 1e-14) - { - flag1 = true; - d[n-i][mmax-i] = -sqtab[n-(i-1)]/sqtab[n-i]*d[n-i+2][mmax-i]; - } - if (abs(d00_2) < 1e-14) - { - flag2 = true; - d[n-i][mmax+i] = -sqtab[n-(i-1)]/sqtab[n-i]*d[n-i+2][mmax+i]; - } - rec1 = prod1[n-i] = fact1/(n*cth-i); - double rec2 = prod2[n-i] = -fact1/(n*cth+i); - for (tsize k=1; k<=n-i; ++k) - { - if ((flag1) && (k==n-i)) continue; - if (abs(prod1[n-i-k+1]) >= 1e14) - prod1[n-i-k+1] = rec1 = 1e16; - prod1[n-i-k] = rec1 = -sqtab[k]/( (i-(n-k)*cth)*xsth + sqtab[k-1]*rec1 ); - } - for (tsize k=1; k<=n-i; ++k) - { - if ((flag2) && (k==n-i)) continue; - if (abs(prod2[n-i-k+1]) >= 1e14) - prod2[n-i-k+1] = rec2 = 1e16; - prod2[n-i-k] = rec2 = -sqtab[k]/( (i+(n-k)*cth)*xsth + sqtab[k-1]*rec2 ); - } - ratio1 = d00; - double ratio2 = d00_2; - if (flag1) ratio1 = d[n-i][mmax-i]; - if (flag2) ratio2 = d[n-i][mmax+i]; - for (tsize l=i; l<=n; ++l) - { - if ((flag1) && (l==i)) continue; - ratio1 *= prod1[l-i]; - double p_li = wigner_xpow(l+i,1.); - d[n-l][mmax-i] = p_li*ratio1; - if (l<=mmax) - d[n-i][mmax-l] = ratio1; - } - for (tsize l=i; l<=n; ++l) - { - if ((flag2) && (l==i)) continue; - ratio2 *= -prod2[l-i]; - double p_li = wigner_xpow(l+i,1.); - d[n-l][mmax+i] = p_li*ratio2; - if (l<=mmax) - d[n-i][mmax+l] = p_li*ratio2; - } - if (i &d ) -{ - - try { - d.fast_alloc(2*n+1,2*mmax+1); - } catch ( std::bad_alloc & e ) { - std::cerr << " Out of memory allocating " << (2*n+1)*(2*mmax+1)*8./1024/1024 << "MB for wignerCalcHalfpi" << std::endl; - throw; - } - - if (n==0) - { - d[0][0]=1; - return; - } - - double d00; - if ((n&1)==0) - { - double rec = d00 = d[n][mmax] = wignerCalc00_halfpi(n); - if (1<=mmax) d[n-1][mmax-1] = rec; - for (tsize k=2; k<=n; k+=2) - { - rec *= -std::sqrt(((n+k-1)*(n-k+2))/double((n-k+1)*(n+k))); - d[n-k][mmax] = rec; - if (k<=mmax) d[n][mmax-k] = rec; - } - for (tsize k=1; k<=n; k+=2) - { - d[n-k][mmax] = 0; - if (k<=mmax) d[n][mmax-k] = 0; - } - } - else - { - double rec = wignerCalc00_halfpi(n-1)*std::sqrt(n/(n+1.)); - d[n-1][mmax] = -rec; - if (1<=mmax) d[n][mmax-1] = rec; - for (tsize k=3; k<=n; k+=2) - { - rec *= -std::sqrt((n+k-1)*(n-k+2)/double((n-k+1)*(n+k))); - d[n-k][mmax] = -rec; - if (k<=mmax) d[n][mmax-k] = rec; - } - for (tsize k=0; k<=n; k+=2) - { - d[n-k][mmax] = 0; - if (k<=mmax) d[n][mmax-k] = 0; - } - d00 = d[n-1][mmax]; - } - - levels:: arr prod(n+1), sqtab(n); - for (tsize i=0; i &d ) -{ - try { - d.fast_alloc (2*n+1,2*mmax+1); - } catch ( std::bad_alloc & e ) { - std::cerr << " Out of memory allocating " << (2*n+1)*(2*mmax+1)*8./1024/1024 << "MB for wignerCalc" << std::endl; - throw; - } - - tdiff mm = mmax; - if (n == 0) - { - d[0][0]=1.; - return; - } - if (abs(theta) >= twopi) - { - double thetasign = 1.; - if (theta < 0) thetasign = -1.; - long anglefac(0); - double anglediff = 2.*twopi; - while( abs(anglediff) > twopi ) - { - anglefac++; - anglediff = theta - thetasign*anglefac*twopi; - } - theta = theta - thetasign*anglefac*twopi; - } - if ( (levels::abs_approx(theta, 0., 1e-14)) || (levels::abs_approx(fabs(theta), twopi, conv_acc) )) - { - d.fill(0.); - for (tdiff k=-mm; k<=mm; ++k) - d[n+k][mmax+k] = 1.; - return; - } - if (levels::abs_approx(theta, pi, 1e-14)) - { - d.fill(0.); - for (tdiff k=-mm; k<=mm; ++k) - d[n+k][mmax+k] = wigner_xpow(n+k,1.); - return; - } - if ( levels::abs_approx(theta, halfpi, 1e-14) ) - { - wignerCalcHalfpi( n, mmax, d ); - return; - } - wignerCalcGeneral( n, mmax, theta, d ); -} - -} // namespace conviqt