Changes in / [bea7acb:0d0aee1] in sasmodels
- Files:
-
- 1 deleted
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
explore/J1c.py
r0a6da3c re6f1410 8 8 import pylab 9 9 10 mp.dps = 150 # number of digits to use in estimating true value 10 11 11 12 SHOW_DIFF = True # True if show diff rather than function value 12 13 LINEAR_X = False # True if q is linearly spaced instead of log spaced 13 14 14 def mp_J1c(vec , bits=500):15 def mp_J1c(vec): 15 16 """ 16 17 Direct calculation using sympy multiprecision library. 17 18 """ 18 with mp.workprec(bits): 19 return [_mp_J1c(mp.mpf(x)) for x in vec] 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 as J131 from scipy.special import 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.float3240 39 x = np.asarray(x, dtype) 41 40 ans = np.empty_like(x) … … 43 42 44 43 # Branch a 45 a_idx = ax < f(8.0)44 a_idx = ax < 8.0 46 45 a_xsq = x[a_idx]**2 47 46 a_coeff1 = list(reversed((72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606))) 48 47 a_coeff2 = list(reversed((144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0))) 49 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_ans248 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_ans2 52 51 53 52 # Branch b … … 56 55 b_x = x[b_idx] 57 56 58 b_y = f(64.0)/(b_ax**2)59 b_xx = b_ax - f(2.356194491)57 b_y = 64.0/(b_ax**2) 58 b_xx = b_ax - 2.356194491 60 59 b_coeff1 = list(reversed((1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6))) 61 60 b_coeff2 = list(reversed((0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6))) 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)61 b_ans1 = np.polyval(b_coeff1[n:], b_y) 62 b_ans2 = np.polyval(b_coeff2[n:], b_y) 64 63 b_sn, b_cn = np.sin(b_xx), np.cos(b_xx) 65 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_x64 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_x 66 65 67 66 return ans 68 69 def div_J1c(x, dtype):70 f = np.float64 if np.dtype(dtype) == np.float64 else np.float3271 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*x78 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/ans291 else:92 y = f(64.0)/(ax*ax)93 xx = ax - f(2.356194491)94 ans1 = (f(1.0)95 + y*(f(0.183105e-2)96 + y*(f(-0.3516396496e-4)97 + y*(f(0.2457520174e-5)98 + y*f(-0.240337019e-6)))))99 ans2 = (f(0.04687499995)100 + y*(f(-0.2002690873e-3)101 + y*(f(0.8449199096e-5)102 + y*(f(-0.88228987e-6)103 + y*f(0.105787412e-6)))))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 ans107 67 108 68 def plotdiff(x, target, actual, label): … … 124 84 """ 125 85 target = np.asarray(mp_J1c(x), 'double') 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) 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) 131 92 pylab.xlabel("qr (1/Ang)") 132 93 if SHOW_DIFF: … … 156 117 157 118 if __name__ == "__main__": 158 #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))119 print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3])) 159 120 main() 160 121 pylab.show() -
sasmodels/models/power_law.py
rb15849c r841753c 1 1 #power_law model 2 2 #conversion of PowerLawAbsModel.py 3 #converted by Steve King, Dec 20154 3 5 4 r"""
Note: See TracChangeset
for help on using the changeset viewer.