| 1 | r""" |
|---|
| 2 | Show numerical precision of cylinder form. |
|---|
| 3 | |
|---|
| 4 | Using:: |
|---|
| 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 | |
|---|
| 11 | import numpy as np |
|---|
| 12 | from sympy.mpmath import mp |
|---|
| 13 | #import matplotlib; matplotlib.use('TkAgg') |
|---|
| 14 | import pylab |
|---|
| 15 | |
|---|
| 16 | SHOW_DIFF = True # True if show diff rather than function value |
|---|
| 17 | LINEAR_X = False # True if q is linearly spaced instead of log spaced |
|---|
| 18 | |
|---|
| 19 | RADIUS = 20 |
|---|
| 20 | LENGTH = 300 |
|---|
| 21 | CONTRAST = 5 |
|---|
| 22 | THETA = 45 |
|---|
| 23 | |
|---|
| 24 | def 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 | |
|---|
| 31 | def _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 | |
|---|
| 41 | def 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 | |
|---|
| 53 | def 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 | |
|---|
| 61 | def 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 | |
|---|
| 66 | def _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 | |
|---|
| 100 | def 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 | |
|---|
| 113 | def 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 | |
|---|
| 130 | def 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 | |
|---|
| 147 | if __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() |
|---|
