oldtests/transcoeff_cruzan.py

   1 # [Xu] = Journal of computational physics 139, 137–165
   2
   3 from __future__ import print_function
   4
   5 def p_q(q, n, nu):
   6     return n + nu - 2*q
   7
   8 def qmax(M, n, mu, nu):
   9     return floor(min(n, nu, (n+nu-abs(M+mu))/2))
  10
  11 def Qmax(M, n, mu, nu): # [Xu](60)
  12     return floor(min(n, nu, (n+nu+1-abs(M+mu))/2))
  13
  14 def gaunta_p(M, n, mu, nu, p): # [Xu](5)
  15     #print (M,n,mu,nu,p, file=sys.stderr)
  16     return (-1)**(M+mu) * (2*p +1) * sqrt(
  17         factorial(n+M) * factorial(nu+mu) * factorial(p-M-mu)
  18       / factorial(n-M) / factorial(nu-mu) / factorial(p+M+mu)) * (
  19         wigner_3j(n, nu, p, 0, 0, 0) * wigner_3j(n, nu, p, M, mu, -M-mu))
  20
  21 def bCXcoeff(M, n, mu, nu, p): # [Xu](61)
  22     #print(M,n,mu,nu,p,file=sys.stderr)
  23     return (-1)**(M+mu) * (2*p + 3) * sqrt(
  24         factorial(n+M) * factorial(nu+mu) * factorial(p+1-M-mu)
  25       / factorial(n-M) / factorial(nu-mu) / factorial(p+1+M+mu)) * (
  26         wigner_3j(n, nu, p, 0, 0, 0) * wigner_3j(n, nu, p+1, M, mu, -M-mu))
  27
  28 def ACXcoeff(m, n, mu, nu, q): # [Xu](58)
  29     p = p_q(q,n,nu)
  30     return ((-1)**m * (2*nu + 1) * factorial(n+m) * factorial(nu-mu) / (
  31         2 * n * (n+1) * factorial(n-m) * factorial(nu+mu)) * I**p *
  32         (n*(n+1) + nu*(nu+1) - p*(p+1)) * gaunta_p(-m,n,mu,nu,p))
  33
  34 def BCXcoeff(m, n, mu, nu, q): # [Xu](59)
  35     p = p_q(q,n,nu)
  36     return ((-1)**(m+1) * (2*nu + 1) * factorial(n+m) * factorial(nu-mu) / (
  37         2 * n * (n+1) * factorial(n-m) * factorial(nu+mu)) * I**(p+1) *
  38         sqrt(((p+1)**2-(n-nu)**2) * ((n+nu+1)**2-(p+1)**2))
  39         * bCXcoeff(-m,n,mu,nu,p))
  40
  41 def printACXcoeffs(lMax, file=sys.stdout):
  42     for n in IntegerRange(lMax+1):
  43         for nu in IntegerRange(lMax+1):
  44             for m in IntegerRange(-n, n+1):
  45                 for mu in IntegerRange(-nu, nu+1):
  46                     for q in IntegerRange(qmax(-m,n,mu,nu)):
  47                         #print(m, n, mu, nu, q, p_q(q,n,nu), file=sys.stderr)
  48                         coeff= ACXcoeff(m, n, mu, nu, q);
  49                         print(N(coeff, prec=53),
  50                             ", // %d, %d, %d, %d, %d," % (m,n,mu,nu,q),
  51                             coeff,
  52                             file=file)
  53     return
  54
  55 def printBCXcoeffs(lMax, file=sys.stdout):
  56     for n in IntegerRange(lMax+1):
  57         for nu in IntegerRange(lMax+1):
  58             for m in IntegerRange(-n, n+1):
  59                 for mu in IntegerRange(-nu, nu+1):
  60                     for q in IntegerRange(1, Qmax(-m,n,mu,nu) +1 ):
  61                         #print(m, n, mu, nu, q, p_q(q,n,nu), file=sys.stderr)
  62                         coeff= BCXcoeff(m, n, mu, nu, q);
  63                         print(N(coeff, prec=53),
  64                             ", // %d, %d, %d, %d, %d," % (m,n,mu,nu,q),
  65                             coeff,
  66                             file=file)
  67     return
  68
  69 sphericalBessels = (None,
  70         spherical_bessel_J,
  71         spherical_bessel_Y,
  72         spherical_hankel1,
  73         spherical_hankel2
  74         )
  75
  76 # N.B. sage's gen_legendre_P _does_ include (-1)**m Condon-Shortley phase
  77 # whereas formulae in [Xu] do not.
  78 def trcoeff_ACX(m, n, mu, nu, besseltype, kd, th, fi, csphase=1): # [Xu](58)
  79     res = 0
  80     for q in range(qmax(-m,n,mu,nu)+1):
  81         p = p_q(q,n,nu)
  82         res += ACXcoeff(m,n,mu,nu,q) * sphericalBessels[besseltype](p,kd) * gen_legendre_P(p, mu-m, cos(th)) * (-csphase)**(mu-m) # compensate for csphase
  83     res *= exp(I*(mu-m)*fi)
  84     return res
  85
  86 def trcoeff_BCX(m, n, mu, nu, besseltype, kd, th, fi, csphase=1): # [Xu](59)
  87     res = 0
  88     for q in IntegerRange(1,Qmax(-m,n,mu,nu)+1):
  89         p = p_q(q,n,nu)
  90         res += BCXcoeff(m,n,mu,nu,q) * sphericalBessels[besseltype](p+1,kd) * gen_legendre_P(p+1, mu-m, cos(th)) * (-csphase)**(mu-m)
  91     res *= exp(I*(mu-m)*fi)
  92     return res
  93
  94
  95 def legpi_xu(n, m, fi): # momentálně neošetřeny okraje (cos(fi) == +- 1)
  96     return m/sin(fi) * gen_legendre_P(n, m, cos(fi))
  97
  98 def legtau_xu(n, m, fi):
  99     locx = var('locx')
 100     return -sin(fi)*derivative(gen_legendre_P(n,m,locx), locx).substitute(locx = cos(fi))
 101
 102 def vswf_M_xu(besseltype, n, m, kr, th, fi):
 103     postpart = sphericalBessels[besseltype](n, kr) * exp(I * m * fi)
 104     tc = I*legpi_xu(n,m,th) * postpart
 105     fc = -legtau_xu(n,m,th) * postpart
 106     return (0, tc, fc)
 107
 108 def vswf_N_xu(besseltype, n, m, kr, th, fi):
 109     eimf = exp(I * m * fi)
 110     rc = n * (n+1) * gen_legendre_P(n, m, cos(th)) * sphericalBessels[besseltype](n, kr)/kr * eimf
 111     krv = var('krv')
 112     radpart = derivative(krv * sphericalBessels[besseltype](n, krv), krv).substitute(krv=kr)/kr
 113     tc = legtau_xu(n,m,th) * radpart * eimf
 114     fc = I*legpi_xu(n,m,th)  * radpart * eimf
 115     return (rc,tc,fc)
 116
 117 def cart2sph(v):
 118     (x, y, z) = v
 119     r = sqrt(x**2 + y**2 + z**2)
 120     th = arccos(z/r) if r else 0
 121     fi = arctan2(y,x)
 122     return (r, th, fi)
 123
 124 def sph2cart(s):
 125     (r, th, fi) = s
 126     sinth = sin(th)
 127     x = r * sinth * cos(fi)
 128     y = r * sinth * sin(fi)
 129     z = r * cos(th)
 130     return (x,y,z)
 131
 132 def sphvec2cart(loccart, sph):
 133     r, th, fi = sph
 134     sinth = sin(th)
 135     costh = cos(th)
 136     sinfi = sin(fi)
 137     cosfi = cos(fi)
 138
 139     rx = sinth * cosfi
 140     ry = sinth * sinfi
 141     rz = costh
 142     tx = costh * cosfi
 143     ty = costh * sinfi
 144     tz = -sinth
 145     fx = -sinfi
 146     fy = cosfi
 147     fz = 0
 148
 149     rc, tc, fc = loccart
 150     x = rx * rc + tx * tc + fx * fc
 151     y = ry * rc + ty * tc + fy * fc
 152     z = rz * rc + tz * tc + fz * fc
 153     return (x, y, z)
 154
 155 def cart2sphvec(cart, sph):
 156     _, th, fi = sph
 157     x, y, z = cart
 158
 159     rx = sinth * cosfi
 160     ry = sinth * sinfi
 161     rz = costh
 162     tx = costh * cosfi
 163     ty = costh * sinfi
 164     tz = -sinth
 165     fx = -sinfi
 166     fy = cosfi
 167     fz = 0
 168
 169     rc = rx * x + ry * y + rz * z
 170     tc = tx * x + ty * y + tz * z
 171     fc = fx * x + fy * y + fz * z
 172     return (rc, tc, fc)
 173
 174 def test_M_translation_xu(lMax, origl, origm, origcartat, cartshift):
 175     ox, oy, oz = origcartat
 176     sx, sy, sz = cartshift
 177     newcartat = (ox - sx, oy - sy, oz - sz)
 178     w1s = cart2sph(origcartat)
 179     w2s = cart2sph(newcartat)
 180     pass # TODO
 181
 182
 183