oldtests/bessel_precision/sagebesselgen.py

   1 from __future__ import print_function
   2
   3 def printhankelrow(lMax, x, file=sys.stdout):
   4     print(lMax, N(x), end=' ', file=file);
   5     for l in range(lMax+1):
   6         print(N(spherical_hankel1(l,x)), end = ' ', file=file)
   7     print('', file=file)
   8
   9
  10 def printbesselJrow(lMax, x, file=sys.stdout):
  11     print(lMax, N(x), end=' ', file=file);
  12     for l in range(lMax+1):
  13         print(N(spherical_bessel_J(l,x)), end = ' ', file=file)
  14     print('', file=file)
  15
  16
  17 def printbesselYrow(lMax, x, file=sys.stdout):
  18     print(lMax, N(x), end=' ', file=file);
  19     for l in range(lMax+1):
  20         print(N(spherical_bessel_Y(l,x)), end = ' ', file=file)
  21     print('', file=file)
  22
  23
  24 #cf DLMF 10.51.2
  25
  26 def printbesselDJrow(lMax, x, file=sys.stdout):
  27     print(lMax, N(x), end=' ', file=file);
  28     for l in range(lMax+1):
  29         print(N(-spherical_bessel_J(l+1,x)
  30             + l/x*spherical_bessel_J(l,x)), end = ' ', file=file)
  31     print('', file=file)
  32
  33
  34 def printbesselDYrow(lMax, x, file=sys.stdout):
  35     print(lMax, N(x), end=' ', file=file);
  36     for l in range(lMax+1):
  37         print(N(-spherical_bessel_Y(l+1,x)
  38             + l/x*spherical_bessel_Y(l,x)), end = ' ', file=file)
  39     print('', file=file)
  40
  41
  42 def ank(n, k):
  43     return factorial(n+k)/2**k/factorial(k)/factorial(n-k)
  44
  45 def genall(lMax):
  46     f = open('besselDJcases', 'w')
  47     for o in IntegerRange(1,100):
  48         printbesselDJrow(lMax, o, file=f)
  49         printbesselDJrow(lMax, 1/o, file=f)
  50         printbesselDJrow(lMax, o/sqrt(3), file=f)
  51     f = open('besselDYcases', 'w')
  52     for o in IntegerRange(1,100):
  53         printbesselDYrow(lMax, o, file=f)
  54         printbesselDYrow(lMax, 1/o, file=f)
  55         printbesselDYrow(lMax, o/sqrt(3), file=f)
  56     f = open('besselJcases', 'w')
  57     for o in IntegerRange(1,100):
  58         printbesselJrow(lMax, o, file=f)
  59         printbesselJrow(lMax, 1/o, file=f)
  60         printbesselJrow(lMax, o/sqrt(3), file=f)
  61     f = open('besselYcases', 'w')
  62     for o in IntegerRange(1,100):
  63         printbesselYrow(lMax, o, file=f)
  64         printbesselYrow(lMax, 1/o, file=f)
  65         printbesselYrow(lMax, o/sqrt(3), file=f)
  66
  67
  68 import math
  69 M_LN2 = math.log(2)
  70
  71 def ankf(n,k):
  72     n = float(n)
  73     k = float(k)
  74     return math.exp(math.lgamma(n+k+1) - k * M_LN2 - math.lgamma(k+1) - math.lgamma(n-k+1))
  75
  76 def ankrelerr(n,k):
  77     a = ank(n,k)
  78     b = ankf(n,k)
  79     return 2 * abs(a - b)/(abs(a)+abs(b))