[cb038a2] | 1 | """ |
---|
| 2 | The current 1D calculations for BCC paracrystal are very wrong at low q, orders |
---|
| 3 | of magnitude wrong. The integration fails to capture a very narrow, |
---|
| 4 | very steep ridge. |
---|
| 5 | |
---|
| 6 | Uncomment the plot() line at the bottom of the code to show an image of the |
---|
| 7 | set of S(q, theta, phi) values that must be summed to form the 1D S(q=0.001) |
---|
| 8 | for the BCC paracrystal form. Note particularly that this is a log scale |
---|
| 9 | image spanning 10 orders of magnitude. This pattern repeats itself 8 times |
---|
| 10 | over the entire 4 pi surface integral. |
---|
| 11 | |
---|
| 12 | You can explore various integration options by uncommenting more lines. |
---|
| 13 | Adaptive integration using scpy.integrate.dbsquad is very slow. Romberg |
---|
| 14 | didn't even complete in the time I gave it. Accurate brute force calculation |
---|
| 15 | requires a 4000x4000 grid to get enough precision. |
---|
| 16 | |
---|
| 17 | We may need a specialized integrator for low q which can identify and integrate |
---|
| 18 | the ridges properly. |
---|
| 19 | |
---|
| 20 | This program need sasmodels on the path so it is inserted automatically, |
---|
| 21 | assuming that the explore/bccpy.py is beside sasmodels/special.py in the |
---|
| 22 | source tree. Run from the sasmodels directory using: |
---|
| 23 | |
---|
| 24 | python explore/bccpy.py |
---|
| 25 | """ |
---|
| 26 | |
---|
[859567e] | 27 | from __future__ import print_function, division |
---|
| 28 | |
---|
[cb038a2] | 29 | import os, sys |
---|
| 30 | sys.path.insert(0, os.path.dirname(os.path.dirname(__file__))) |
---|
| 31 | |
---|
[859567e] | 32 | import numpy as np |
---|
| 33 | from numpy import pi, sin, cos, exp, expm1, degrees, log10 |
---|
| 34 | from scipy.integrate import dblquad, simps, romb, romberg |
---|
| 35 | import pylab |
---|
| 36 | |
---|
| 37 | from sasmodels.special import square |
---|
| 38 | from sasmodels.special import Gauss20Wt, Gauss20Z |
---|
| 39 | from sasmodels.special import Gauss76Wt, Gauss76Z |
---|
| 40 | from sasmodels.special import Gauss150Wt, Gauss150Z |
---|
| 41 | |
---|
| 42 | Q = 0.001 |
---|
| 43 | DNN = 220. |
---|
| 44 | D_FACTOR = 0.06 |
---|
| 45 | RADIUS = 40.0 |
---|
| 46 | SLD = 3.0 |
---|
| 47 | SLD_SOLVENT = 6.3 |
---|
| 48 | |
---|
| 49 | def kernel(q, dnn, d_factor, theta, phi): |
---|
[cb038a2] | 50 | """ |
---|
| 51 | S(q) kernel for paracrystal forms. |
---|
| 52 | """ |
---|
[859567e] | 53 | qab = q*sin(theta) |
---|
| 54 | qa = qab*cos(phi) |
---|
| 55 | qb = qab*sin(phi) |
---|
| 56 | qc = q*cos(theta) |
---|
| 57 | |
---|
| 58 | if 0: # sc |
---|
| 59 | a1, a2, a3 = qa, qb, qc |
---|
| 60 | dcos = dnn |
---|
| 61 | if 1: # bcc |
---|
| 62 | a1 = +qa - qc + qb |
---|
| 63 | a2 = +qa + qc - qb |
---|
| 64 | a3 = -qa + qc + qb |
---|
| 65 | dcos = dnn/2 |
---|
[cb038a2] | 66 | if 0: # fcc |
---|
| 67 | a1 = qb + qa |
---|
| 68 | a2 = qa + qc |
---|
| 69 | a3 = qb + qc |
---|
| 70 | dcos = dnn/2 |
---|
[859567e] | 71 | |
---|
| 72 | arg = 0.5*square(dnn*d_factor)*(a1**2 + a2**2 + a3**2) |
---|
| 73 | exp_arg = exp(-arg) |
---|
| 74 | den = [((exp_arg - 2*cos(dcos*a))*exp_arg + 1.0) for a in (a1, a2, a3)] |
---|
| 75 | Sq = -expm1(-2*arg)**3/np.prod(den, axis=0) |
---|
| 76 | return Sq |
---|
| 77 | |
---|
| 78 | |
---|
| 79 | def scipy_dblquad(q=Q, dnn=DNN, d_factor=D_FACTOR): |
---|
[cb038a2] | 80 | """ |
---|
| 81 | Compute the integral using scipy dblquad. This gets the correct answer |
---|
| 82 | eventually, but it is slow. |
---|
| 83 | """ |
---|
| 84 | evals = [0] |
---|
[859567e] | 85 | def integrand(theta, phi): |
---|
[cb038a2] | 86 | evals[0] += 1 |
---|
[859567e] | 87 | Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi) |
---|
| 88 | return Sq*sin(theta) |
---|
[cb038a2] | 89 | ans = dblquad(integrand, 0, pi/2, lambda x: 0, lambda x: pi/2)[0]*8/(4*pi) |
---|
| 90 | print("dblquad evals =", evals[0]) |
---|
| 91 | return ans |
---|
[859567e] | 92 | |
---|
| 93 | |
---|
[cb038a2] | 94 | def scipy_romberg_2d(q=Q, dnn=DNN, d_factor=D_FACTOR): |
---|
| 95 | """ |
---|
| 96 | Compute the integral using romberg integration. This function does not |
---|
| 97 | complete in a reasonable time. No idea if it is accurate. |
---|
| 98 | """ |
---|
[859567e] | 99 | def inner(phi, theta): |
---|
| 100 | return kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi) |
---|
| 101 | def outer(theta): |
---|
| 102 | return romberg(inner, 0, pi/2, divmax=100, args=(theta,))*sin(theta) |
---|
| 103 | return romberg(outer, 0, pi/2, divmax=100)*8/(4*pi) |
---|
| 104 | |
---|
| 105 | |
---|
[cb038a2] | 106 | def semi_romberg(q=Q, dnn=DNN, d_factor=D_FACTOR, n=100): |
---|
| 107 | """ |
---|
| 108 | Use 1D romberg integration in phi and regular simpsons rule in theta. |
---|
| 109 | """ |
---|
| 110 | evals = [0] |
---|
| 111 | def inner(phi, theta): |
---|
| 112 | evals[0] += 1 |
---|
| 113 | return kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi) |
---|
| 114 | theta = np.linspace(0, pi/2, n) |
---|
| 115 | f_phi = [romberg(inner, 0, pi/2, divmax=100, args=(t,)) |
---|
| 116 | for t in theta] |
---|
| 117 | ans = simps(sin(theta)*np.array(f_phi), dx=theta[1]-theta[0]) |
---|
| 118 | print("semi romberg evals =", evals[0]) |
---|
| 119 | return ans*8/(4*pi) |
---|
| 120 | |
---|
[859567e] | 121 | def gauss_quad(q=Q, dnn=DNN, d_factor=D_FACTOR, n=150): |
---|
[cb038a2] | 122 | """ |
---|
| 123 | Compute the integral using gaussian quadrature for n = 20, 76 or 150. |
---|
| 124 | """ |
---|
[859567e] | 125 | if n == 20: |
---|
| 126 | z, w = Gauss20Z, Gauss20Wt |
---|
| 127 | elif n == 76: |
---|
| 128 | z, w = Gauss76Z, Gauss76Wt |
---|
| 129 | else: |
---|
| 130 | z, w = Gauss150Z, Gauss150Wt |
---|
| 131 | theta = pi/4*(z + 1) |
---|
| 132 | phi = pi/4*(z + 1) |
---|
| 133 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
| 134 | Aw = w[None, :] * w[:, None] |
---|
| 135 | sin_theta = np.fmax(abs(sin(Atheta)), 1e-6) |
---|
| 136 | Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) |
---|
[cb038a2] | 137 | print("gauss %d evals ="%n, n**2) |
---|
[859567e] | 138 | return np.sum(Sq*Aw*sin_theta)*8/(4*pi) |
---|
| 139 | |
---|
| 140 | |
---|
| 141 | def gridded_integrals(q=0.001, dnn=DNN, d_factor=D_FACTOR, n=300): |
---|
[cb038a2] | 142 | """ |
---|
| 143 | Compute the integral on a regular grid using rectangular, trapezoidal, |
---|
| 144 | simpsons, and romberg integration. Romberg integration requires that |
---|
| 145 | the grid be of size n = 2**k + 1. |
---|
| 146 | """ |
---|
[859567e] | 147 | theta = np.linspace(0, pi/2, n) |
---|
| 148 | phi = np.linspace(0, pi/2, n) |
---|
| 149 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
| 150 | Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) |
---|
| 151 | Sq *= abs(sin(Atheta)) |
---|
| 152 | dx, dy = theta[1]-theta[0], phi[1]-phi[0] |
---|
| 153 | print("rect", n, np.sum(Sq)*dx*dy*8/(4*pi)) |
---|
| 154 | print("trapz", n, np.trapz(np.trapz(Sq, dx=dx), dx=dy)*8/(4*pi)) |
---|
| 155 | print("simpson", n, simps(simps(Sq, dx=dx), dx=dy)*8/(4*pi)) |
---|
| 156 | print("romb", n, romb(romb(Sq, dx=dx), dx=dy)*8/(4*pi)) |
---|
[cb038a2] | 157 | print("gridded %d evals ="%n, n**2) |
---|
[859567e] | 158 | |
---|
| 159 | def plot(q=0.001, dnn=DNN, d_factor=D_FACTOR, n=300): |
---|
[cb038a2] | 160 | """ |
---|
| 161 | Plot the 2D surface that needs to be integrated in order to compute |
---|
| 162 | the BCC S(q) at a particular q, dnn and d_factor. *n* is the number |
---|
| 163 | of points in the grid. |
---|
| 164 | """ |
---|
| 165 | theta = np.linspace(0, pi, n) |
---|
| 166 | phi = np.linspace(0, 2*pi, n) |
---|
| 167 | #theta = np.linspace(0, pi/2, n) |
---|
| 168 | #phi = np.linspace(0, pi/2, n) |
---|
[859567e] | 169 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
| 170 | Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) |
---|
| 171 | Sq *= abs(sin(Atheta)) |
---|
| 172 | pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Sq, 1.e-6))) |
---|
[cb038a2] | 173 | pylab.axis('tight') |
---|
| 174 | pylab.title("BCC S(q) for q=%g, dnn=%g d_factor=%g" % (q, dnn, d_factor)) |
---|
[859567e] | 175 | pylab.xlabel("theta (degrees)") |
---|
| 176 | pylab.ylabel("phi (degrees)") |
---|
| 177 | cbar = pylab.colorbar() |
---|
[cb038a2] | 178 | cbar.set_label('log10 S(q)') |
---|
[859567e] | 179 | pylab.show() |
---|
| 180 | |
---|
| 181 | if __name__ == "__main__": |
---|
[cb038a2] | 182 | #print("gauss", 20, gauss_quad(n=20)) |
---|
| 183 | #print("gauss", 76, gauss_quad(n=76)) |
---|
| 184 | #print("gauss", 150, gauss_quad(n=150)) |
---|
| 185 | #print("dblquad", scipy_dblquad()) |
---|
| 186 | #print("semi romberg", semi_romberg()) |
---|
| 187 | #gridded_integrals(n=2**8+1) |
---|
| 188 | #gridded_integrals(n=2**10+1) |
---|
| 189 | #gridded_integrals(n=2**13+1) |
---|
[859567e] | 190 | #print("romberg", scipy_romberg()) |
---|
[cb038a2] | 191 | plot(n=400) |
---|