src/cantilever-calib/test/test_roots.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=100
   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 rp(a,b,p,q,h) : # plot root
  17     (x,d) = rm(zr(a,b,p,q), h)
  18     if x[0].imag == 0 and x[0] > 0 : plot(x,d,'-x')
  19     print "Root at ", x[0]
  20 def dr(a,b,x) : # attempt at factoring
  21     return (x, (x-zr(a,b,+1,+1)) \
  22               *(x-zr(a,b,-1,+1)) \
  23               *(x-zr(a,b,+1,-1)) \
  24               *(x-zr(a,b,-1,-1))   )
  25 def testroots(a,b,x) :
  26     figure()
  27     hold(True)
  28     (x,dt) = d(a, b, x)
  29     plot(x,dt,'.')
  30     print "For a = %g, b = %g" % (a, b)
  31     rp(a,b,+1,+1, root_height)
  32     rp(a,b,-1,+1, root_height*2)
  33     rp(a,b,+1,-1, root_height*3)
  34     rp(a,b,-1,-1, root_height*4)
  35     title("For a = %g, b = %g" % (a, b))
  36 def testfn(a,b,x) :
  37     figure()
  38     hold(True)
  39     (x,dt) = d(a, b, x)
  40     plot(x,1/dt,'-')
  41     (x,dt) = dr(a, b, x)
  42     plot(x,1/dt,'.')
  43     title("For a = %g, b = %g" % (a, b))
  44
  45 ioff()
  46
  47 plot = semilogy
  48
  49 # critically damped for b^2 = 4a^2
  50 #for b in [0.1, 1, 2, 4, 6] :
  51 #    testroots(a,b,x)
  52 a = 10.0
  53 x = linspace(0,4*a,N-1)
  54 for b in linspace(a/10, 4*a, 5) :
  55     testfn(a,b,x)
  56 a = 1.0
  57 x = linspace(0,4*a,N-1)
  58 for b in linspace(a/10, 4*a, 5) :
  59     testfn(a,b,x)
  60 a = .1
  61 x = linspace(0,4*a,N-1)
  62 for b in linspace(a/10, 4*a, 5) :
  63     testfn(a,b,x)
  64
  65 show()