src/cantilever-calib/test/test_integral.py

   1 import matplotlib
   2 matplotlib.use('Agg') # select backend that doesn't require X Windows
   3 from pylab import *
   4 from numpy import *
   5
   6 N=1000
   7 root_height = 1
   8 i = complex(0,1)
   9
  10 def d(a,b,x) : # know denominator
  11     return (x, ((a**2 - x**2)**2 + (b*x)**2))
  12 def rm(r, height) : # root marker
  13     return ([r, r], float(height) * array([0.0001,1]))
  14 def zr(a,b,p,q) : # p, q \in {-1,1}
  15     return p*i*b/2.0 * (1 + q*sqrt(complex(1,0) - (2*a/b)**2))
  16 def zrc(c,b,p,q) : # p, q \in {-1,1}
  17     return p*i*b/2.0 * (1 + q*sqrt(complex(1,0) - (2*c)**2))
  18 def dr(a,b,x) : # attempt at factoring
  19     c = float(a)/b
  20     return (x, (x-zrc(c,b,+1,+1)) \
  21               *(x-zrc(c,b,-1,+1)) \
  22               *(x-zrc(c,b,+1,-1)) \
  23               *(x-zrc(c,b,-1,-1))   )
  24
  25 def numeric(a,b,x) :
  26     (x, dt) = d(a,b,x)
  27     dx = x[1]-x[0]
  28     return (1/dt).mean()*(x[-1]-x[0])*2 # 2 because we want int from -inf to inf
  29 def theory_residue_thm(a,b) :
  30     c = complex(a,0)/b
  31     zpp = zrc(c,b,+1,+1)
  32     zmp = zrc(c,b,-1,+1)
  33     zpm = zrc(c,b,+1,-1)
  34     zmm = zrc(c,b,-1,-1)
  35     return 2*pi*i* ( 1/((zpp-zpm)*(zpp-zmp)*(zpp-zmm)) + 1/((zpm-zpp)*(zpm-zmp)*(zpm-zmm)) )
  36 def theory_before_plugging_in_roots(a,b) :
  37     c = complex(a,0)/b
  38     zpp = zrc(c,b,+1,+1)
  39     zmp = zrc(c,b,-1,+1)
  40     zpm = zrc(c,b,+1,-1)
  41     zmm = zrc(c,b,-1,-1)
  42     return pi*i/(zpp**2 - zpm**2) * ( 1/zpp - 1/zpm )
  43 def theory_before_mc(a,b) :
  44     c = complex(a,0)/b
  45     zpp = zrc(c,b,+1,+1)
  46     zmp = zrc(c,b,-1,+1)
  47     zpm = zrc(c,b,+1,-1)
  48     zmm = zrc(c,b,-1,-1)
  49     sq = sqrt(1-4*c**2)
  50     return (-4*pi/b**2) \
  51            / ((1+sq)**2 - (1-sq)**2)  *  ((1-sq) - (1+sq)) \
  52                                          /(b/2*(1+sq)*(1-sq))
  53 def theory_final(a,b) :
  54     return pi/(complex(b)* a**2)
  55
  56 theory = theory_final
  57
  58 PLOT_COMPARISONS = False
  59
  60 def test(a) :
  61     x = linspace(0,8*a,N-1)
  62     nA = [] ; tA = []
  63     if PLOT_COMPARISONS == True :
  64         figure()
  65     for b in linspace(a/10.0, 8*a, 500) :
  66         if PLOT_COMPARISONS == True :
  67             (x1, dt) = d(a,b,x)
  68             plot(x1, 1/dt, 'r.')
  69             (x1, dt) = dr(a,b,x)
  70             plot(x1, 1/dt, 'b-')
  71         nA.append(numeric(a,b,x))
  72         tA.append(theory(a,b))
  73     figure()
  74     plot(nA, 'r.', tA, 'b-')
  75
  76 test(10)
  77 test(1)
  78 test(0.1)
  79
  80 show()