Changeset 0a6da3c in sasmodels
 Timestamp:
 Mar 8, 2016 5:00:20 PM (7 years ago)
 Branches:
 master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 b15849c
 Parents:
 5a483877
 Location:
 explore
 Files:

 1 added
 1 edited
Legend:
 Unmodified
 Added
 Removed

explore/J1c.py
re6f1410 r0a6da3c 8 8 import pylab 9 9 10 mp.dps = 150 # number of digits to use in estimating true value11 10 12 11 SHOW_DIFF = True # True if show diff rather than function value 13 12 LINEAR_X = False # True if q is linearly spaced instead of log spaced 14 13 15 def mp_J1c(vec ):14 def mp_J1c(vec, bits=500): 16 15 """ 17 16 Direct calculation using sympy multiprecision library. 18 17 """ 19 return [_mp_J1c(mp.mpf(x)) for x in vec] 18 with mp.workprec(bits): 19 return [_mp_J1c(mp.mpf(x)) for x in vec] 20 20 21 21 def _mp_J1c(x): … … 25 25 return mp.mpf(2)*mp.j1(x)/x 26 26 27 def np_ j1c(x, dtype):27 def np_J1c(x, dtype): 28 28 """ 29 29 Direct calculation using scipy. 30 30 """ 31 from scipy.special import j1 31 from scipy.special import j1 as J1 32 32 x = np.asarray(x, dtype) 33 return np.asarray(2, dtype)* j1(x)/x33 return np.asarray(2, dtype)*J1(x)/x 34 34 35 def cephes_ j1c(x, dtype, n):35 def cephes_J1c(x, dtype, n): 36 36 """ 37 37 Calculation using pade approximant. 38 38 """ 39 f = np.float64 if np.dtype(dtype) == np.float64 else np.float32 39 40 x = np.asarray(x, dtype) 40 41 ans = np.empty_like(x) … … 42 43 43 44 # Branch a 44 a_idx = ax < 8.045 a_idx = ax < f(8.0) 45 46 a_xsq = x[a_idx]**2 46 47 a_coeff1 = list(reversed((72362614232.0, 7895059235.0, 242396853.1, 2972611.439, 15704.48260, 30.16036606))) 47 48 a_coeff2 = list(reversed((144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0))) 48 a_ans1 = np.polyval( a_coeff1[n:], a_xsq)49 a_ans2 = np.polyval( a_coeff2[n:], a_xsq)50 ans[a_idx] = 2*a_ans1/a_ans249 a_ans1 = np.polyval(np.asarray(a_coeff1[n:], dtype), a_xsq) 50 a_ans2 = np.polyval(np.asarray(a_coeff2[n:], dtype), a_xsq) 51 ans[a_idx] = f(2.0)*a_ans1/a_ans2 51 52 52 53 # Branch b … … 55 56 b_x = x[b_idx] 56 57 57 b_y = 64.0/(b_ax**2)58 b_xx = b_ax  2.35619449158 b_y = f(64.0)/(b_ax**2) 59 b_xx = b_ax  f(2.356194491) 59 60 b_coeff1 = list(reversed((1.0, 0.183105e2, 0.3516396496e4, 0.2457520174e5, 0.240337019e6))) 60 61 b_coeff2 = list(reversed((0.04687499995, 0.2002690873e3, 0.8449199096e5, 0.88228987e6, 0.105787412e6))) 61 b_ans1 = np.polyval( b_coeff1[n:],b_y)62 b_ans2 = np.polyval( b_coeff2[n:], b_y)62 b_ans1 = np.polyval(np.asarray(b_coeff1[n:], dtype),b_y) 63 b_ans2 = np.polyval(np.asarray(b_coeff2[n:], dtype), b_y) 63 64 b_sn, b_cn = np.sin(b_xx), np.cos(b_xx) 64 ans[b_idx] = np.sign(b_x)*np.sqrt( 0.636619772/b_ax) * (b_cn*b_ans1  (8.0/b_ax)*b_sn*b_ans2)*2.0/b_x65 ans[b_idx] = np.sign(b_x)*np.sqrt(f(0.636619772)/b_ax) * (b_cn*b_ans1  (f(8.0)/b_ax)*b_sn*b_ans2)*f(2.0)/b_x 65 66 66 67 return ans 68 69 def div_J1c(x, dtype): 70 f = np.float64 if np.dtype(dtype) == np.float64 else np.float32 71 x = np.asarray(x, dtype) 72 return f(2.0)*np.asarray([_J1(xi, f)/xi for xi in x], dtype) 73 74 def _J1(x, f): 75 ax = abs(x) 76 if ax < f(8.0): 77 y = x*x 78 ans1 = x*(f(72362614232.0) 79 + y*(f(7895059235.0) 80 + y*(f(242396853.1) 81 + y*(f(2972611.439) 82 + y*(f(15704.48260) 83 + y*(f(30.16036606))))))) 84 ans2 = (f(144725228442.0) 85 + y*(f(2300535178.0) 86 + y*(f(18583304.74) 87 + y*(f(99447.43394) 88 + y*(f(376.9991397) 89 + y))))) 90 return ans1/ans2 91 else: 92 y = f(64.0)/(ax*ax) 93 xx = ax  f(2.356194491) 94 ans1 = (f(1.0) 95 + y*(f(0.183105e2) 96 + y*(f(0.3516396496e4) 97 + y*(f(0.2457520174e5) 98 + y*f(0.240337019e6))))) 99 ans2 = (f(0.04687499995) 100 + y*(f(0.2002690873e3) 101 + y*(f(0.8449199096e5) 102 + y*(f(0.88228987e6) 103 + y*f(0.105787412e6))))) 104 sn, cn = np.sin(xx), np.cos(xx) 105 ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1  (f(8.0)/ax)*sn*ans2) 106 return ans if (x < f(0.0)) else ans 67 107 68 108 def plotdiff(x, target, actual, label): … … 84 124 """ 85 125 target = np.asarray(mp_J1c(x), 'double') 86 direct = np_j1c(x, precision) 87 approx0 = cephes_j1c(x, precision, 0) 88 approx1 = cephes_j1c(x, precision, 1) 89 plotdiff(x, target, direct, 'direct '+precision) 90 plotdiff(x, target, approx0, 'cephes '+precision) 91 #plotdiff(x, target, approx1, 'reduced '+precision) 126 #plotdiff(x, target, mp_J1c(x, 11), 'mp 11 bits') 127 plotdiff(x, target, np_J1c(x, precision), 'direct '+precision) 128 plotdiff(x, target, cephes_J1c(x, precision, 0), 'cephes '+precision) 129 #plotdiff(x, target, cephes_J1c(x, precision, 1), 'cephes '+precision) 130 #plotdiff(x, target, div_J1c(x, precision), 'cephes 2 J1(x)/x '+precision) 92 131 pylab.xlabel("qr (1/Ang)") 93 132 if SHOW_DIFF: … … 117 156 118 157 if __name__ == "__main__": 119 print "\n".join(str(x) for x in mp_J1c([1e6,1e5,1e4,1e3]))158 #print "\n".join(str(x) for x in mp_J1c([1e6,1e5,1e4,1e3])) 120 159 main() 121 160 pylab.show()
Note: See TracChangeset
for help on using the changeset viewer.