eic-opticks/ana/piecewise2.py

0001 #!/usr/bin/env python
0002 """
0003
0004 Hmm maybe its impossible to get anywhere with integral without first fixing the BetaInverse
0005
0006
0007 In [22]: pw.subs(b, 1.55)
0008 Out[22]: Piecewise((Max(0, 0.542862145685886*e - 3.9544951109741), (e > 7.294) & (e < 7.75)), (0, True))
0009
0010 In [23]: pw.subs(b, 1.)
0011 Out[23]: Piecewise((Max(0, 0.225957188630962*e - 1.06222481205998), (e > 7.294) & (e < 7.75)), (0, True))
0012
0013
0014
0015 """
0016 import numpy as np
0017 import sympy as sym
0018
0019 ri = np.array([
0020        [ 1.55 ,  1.478],
0021        [ 1.795,  1.48 ],
0022        [ 2.105,  1.484],
0023        [ 2.271,  1.486],
0024        [ 2.551,  1.492],
0025        [ 2.845,  1.496],
0026        [ 3.064,  1.499],
0027        [ 4.133,  1.526],
0028        [ 6.2  ,  1.619],
0029        [ 6.526,  1.618],
0030        [ 6.889,  1.527],
0031        [ 7.294,  1.554],
0032        [ 7.75 ,  1.793],
0033        [ 8.267,  1.783],
0034        [ 8.857,  1.664],
0035        [ 9.538,  1.554],
0036        [10.33 ,  1.454],
0037        [15.5  ,  1.454]
0038       ])
0039
0040
0041
0042 if __name__ == '__main__':
0043
0044
0045      e, b = sym.symbols("e b")
0046
0047      i = 11
0048
0049      e0, r0 = ri[i]
0050      e1, r1 = ri[i+1]
0051
0052      em = (e0 + e1)/2.
0053
0054
0055      v0 = ( 1 - b/r0 ) * ( 1 + b/r0 )
0056      v1 = ( 1 - b/r1 ) * ( 1 + b/r1 )
0057
0058      fr = (e-e0)/(e1-e0)
0059
0060      pt = ( sym.Max(v0*(1-fr) + v1*fr,0),  (e > e0) & (e < e1) )
0061      ot = (0, True )
0062      pw = sym.Piecewise( pt, ot )
0063
0064      v = pw.subs(b, 1.55).subs(e, em)
0065      print(v)
0066
0067