# source:sasmodels/explore/J1c.py@e6f1410

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since e6f1410 was e6f1410, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

examine numerical precision of J1c, sinc, j1c and log gamma

• Property mode set to `100644`
File size: 3.5 KB
Line
1r"""
2Show numerical precision of \$2 J_1(x)/x\$.
3"""
4
5import numpy as np
6from sympy.mpmath import mp
7#import matplotlib; matplotlib.use('TkAgg')
8import pylab
9
10mp.dps = 150 # number of digits to use in estimating true value
11
12SHOW_DIFF = True  # True if show diff rather than function value
13LINEAR_X = False  # True if q is linearly spaced instead of log spaced
14
15def mp_J1c(vec):
16    """
17    Direct calculation using sympy multiprecision library.
18    """
19    return [_mp_J1c(mp.mpf(x)) for x in vec]
20
21def _mp_J1c(x):
22    """
23    Helper funciton for mp_j1c
24    """
25    return mp.mpf(2)*mp.j1(x)/x
26
27def np_j1c(x, dtype):
28    """
29    Direct calculation using scipy.
30    """
31    from scipy.special import j1
32    x = np.asarray(x, dtype)
33    return np.asarray(2, dtype)*j1(x)/x
34
35def cephes_j1c(x, dtype, n):
36    """
38    """
39    x = np.asarray(x, dtype)
40    ans = np.empty_like(x)
41    ax = abs(x)
42
43    # Branch a
44    a_idx = ax < 8.0
45    a_xsq = x[a_idx]**2
46    a_coeff1 = list(reversed((72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606)))
47    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_ans2
51
52    # Branch b
53    b_idx = ~a_idx
54    b_ax = ax[b_idx]
55    b_x = x[b_idx]
56
57    b_y = 64.0/(b_ax**2)
58    b_xx = b_ax - 2.356194491
59    b_coeff1 = list(reversed((1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6)))
60    b_coeff2 = list(reversed((0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6)))
61    b_ans1 = np.polyval(b_coeff1[n:], b_y)
62    b_ans2 = np.polyval(b_coeff2[n:], b_y)
63    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_x
65
66    return ans
67
68def plotdiff(x, target, actual, label):
69    """
70    Plot the computed value.
71
72    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
73    """
74    if SHOW_DIFF:
75        err = np.clip(abs((target-actual)/target), 0, 1)
76        pylab.loglog(x, err, '-', label=label)
77    else:
78        limits = np.min(target), np.max(target)
79        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
80
81def compare(x, precision):
82    r"""
83    Compare the different computation methods using the given precision.
84    """
85    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)
92    pylab.xlabel("qr (1/Ang)")
93    if SHOW_DIFF:
94        pylab.ylabel("relative error")
95    else:
96        pylab.ylabel("2 J1(x)/x")
97        pylab.semilogx(x, target,  '-', label="true value")
98    if LINEAR_X:
99        pylab.xscale('linear')
100
101def main():
102    r"""
103    Compare accuracy of different methods for computing \$3 j_1(x)/x\$.
104    :return:
105    """
106    if LINEAR_X:
107        qr = np.linspace(1,1000,2000)
108    else:
109        qr = np.logspace(-3,5,400)
110    pylab.subplot(121)
111    compare(qr, 'single')
112    pylab.legend(loc='best')
113    pylab.subplot(122)
114    compare(qr, 'double')
115    pylab.legend(loc='best')
116    pylab.suptitle('2 J1(x)/x')
117
118if __name__ == "__main__":
119    print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
120    main()
121    pylab.show()
Note: See TracBrowser for help on using the repository browser.