source: sasmodels/explore/cyl.py @ 7c2935c

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 7c2935c was 0a6da3c, checked in by Paul Kienzle <pkienzle@…>, 9 years ago

explore numerical precision of cylinder form

  • Property mode set to 100644
File size: 4.3 KB
Line 
1r"""
2Show numerical precision of cylinder form.
3
4Using::
5
6    qr = q r sin(t)
7    qh = q h/2 cos(t)
8    F = 2 J_1(qr)/qr sin(qh)/qh
9"""
10
11import numpy as np
12from sympy.mpmath import mp
13#import matplotlib; matplotlib.use('TkAgg')
14import pylab
15
16SHOW_DIFF = True  # True if show diff rather than function value
17LINEAR_X = False  # True if q is linearly spaced instead of log spaced
18
19RADIUS = 20
20LENGTH = 300
21CONTRAST = 5
22THETA = 45
23
24def mp_form(vec, bits=500):
25    """
26    Direct calculation using sympy multiprecision library.
27    """
28    with mp.workprec(bits):
29        return [_mp_f(mp.mpf(x)) for x in vec]
30
31def _mp_f(x):
32    """
33    Helper function for mp_j1c
34    """
35    f = mp.mpf
36    theta = f(THETA)*mp.pi/f(180)
37    qr = x * f(RADIUS)*mp.sin(theta)
38    qh = x * f(LENGTH)/f(2)*mp.cos(theta)
39    return (f(2)*mp.j1(qr)/qr * mp.sin(qh)/qh)**f(2)
40
41def np_form(x, dtype):
42    """
43    Direct calculation using scipy.
44    """
45    from scipy.special import j1
46    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
47    x = np.asarray(x, dtype)
48    theta = f(THETA)*f(np.pi)/f(180)
49    qr = x * f(RADIUS)*np.sin(theta)
50    qh = x * f(LENGTH)/f(2)*np.cos(theta)
51    return (f(2)*j1(qr)/qr*np.sin(qh)/qh)**f(2)
52
53def sasmodels_form(x, dtype):
54    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
55    x = np.asarray(x, dtype)
56    theta = f(THETA)*f(np.pi)/f(180)
57    qr = x * f(RADIUS)*np.sin(theta)
58    qh = x * f(LENGTH)/f(2)*np.cos(theta)
59    return (J1c(qr, dtype)*np.sin(qh)/qh)**f(2)
60
61def J1c(x, dtype):
62    x = np.asarray(x, dtype)
63    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32
64    return np.asarray([_J1c(xi, f) for xi in x], dtype)
65
66def _J1c(x, f):
67    ax = abs(x)
68    if ax < f(8.0):
69        y = x*x
70        ans1 = f(2)*(f(72362614232.0)
71                  + y*(f(-7895059235.0)
72                  + y*(f(242396853.1)
73                  + y*(f(-2972611.439)
74                  + y*(f(15704.48260)
75                  + y*(f(-30.16036606)))))))
76        ans2 = (f(144725228442.0)
77                  + y*(f(2300535178.0)
78                  + y*(f(18583304.74)
79                  + y*(f(99447.43394)
80                  + y*(f(376.9991397)
81                  + y)))))
82        return ans1/ans2
83    else:
84        y = f(64.0)/(ax*ax)
85        xx = ax - f(2.356194491)
86        ans1 = (f(1.0)
87                  + y*(f(0.183105e-2)
88                  + y*(f(-0.3516396496e-4)
89                  + y*(f(0.2457520174e-5)
90                  + y*f(-0.240337019e-6)))))
91        ans2 = (f(0.04687499995)
92                  + y*(f(-0.2002690873e-3)
93                  + y*(f(0.8449199096e-5)
94                  + y*(f(-0.88228987e-6)
95                  + y*f(0.105787412e-6)))))
96        sn, cn = np.sin(xx), np.cos(xx)
97        ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) * f(2)/x
98        return -ans if (x < f(0.0)) else ans
99
100def plotdiff(x, target, actual, label):
101    """
102    Plot the computed value.
103
104    Use relative error if SHOW_DIFF, otherwise just plot the value directly.
105    """
106    if SHOW_DIFF:
107        err = np.clip(abs((target-actual)/target), 0, 1)
108        pylab.loglog(x, err, '-', label=label)
109    else:
110        limits = np.min(target), np.max(target)
111        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label)
112
113def compare(x, precision):
114    r"""
115    Compare the different computation methods using the given precision.
116    """
117    target = np.asarray(mp_form(x), 'double')
118    plotdiff(x, target, mp_form(x, bits=11), '11-bit')
119    plotdiff(x, target, np_form(x, precision), 'direct '+precision)
120    plotdiff(x, target, sasmodels_form(x, precision), 'sasmodels '+precision)
121    pylab.xlabel("qr (1/Ang)")
122    if SHOW_DIFF:
123        pylab.ylabel("relative error")
124    else:
125        pylab.ylabel("2 J1(x)/x")
126        pylab.semilogx(x, target,  '-', label="true value")
127    if LINEAR_X:
128        pylab.xscale('linear')
129
130def main():
131    r"""
132    Compare accuracy of different methods for computing $3 j_1(x)/x$.
133    :return:
134    """
135    if LINEAR_X:
136        qr = np.linspace(1,1000,2000)
137    else:
138        qr = np.logspace(-3,5,400)
139    pylab.subplot(121)
140    compare(qr, 'single')
141    pylab.legend(loc='best')
142    pylab.subplot(122)
143    compare(qr, 'double')
144    pylab.legend(loc='best')
145    pylab.suptitle('2 J1(x)/x')
146
147if __name__ == "__main__":
148    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3]))
149    main()
150    pylab.show()
Note: See TracBrowser for help on using the repository browser.