[4f611f1] | 1 | """ |
---|
| 2 | Asymmetric shape integration |
---|
| 3 | """ |
---|
| 4 | |
---|
| 5 | from __future__ import print_function, division |
---|
| 6 | |
---|
| 7 | import os, sys |
---|
| 8 | sys.path.insert(0, os.path.dirname(os.path.dirname(__file__))) |
---|
| 9 | |
---|
| 10 | import numpy as np |
---|
[31eea1f] | 11 | import mpmath as mp |
---|
[4f611f1] | 12 | from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10 |
---|
[31eea1f] | 13 | from numpy.polynomial.legendre import leggauss |
---|
[4f611f1] | 14 | from scipy.integrate import dblquad, simps, romb, romberg |
---|
| 15 | import pylab |
---|
| 16 | |
---|
| 17 | from sasmodels.special import square |
---|
| 18 | from sasmodels.special import Gauss20Wt, Gauss20Z |
---|
| 19 | from sasmodels.special import Gauss76Wt, Gauss76Z |
---|
| 20 | from sasmodels.special import Gauss150Wt, Gauss150Z |
---|
| 21 | from sasmodels.special import sas_2J1x_x, sas_sinx_x, sas_3j1x_x |
---|
| 22 | |
---|
[31eea1f] | 23 | def mp_3j1x_x(x): |
---|
| 24 | return 3*(mp.sin(x)/x - mp.cos(x))/(x*x) |
---|
| 25 | def mp_2J1x_x(x): |
---|
| 26 | return 2*mp.j1(x)/x |
---|
| 27 | def mp_sinx_x(x): |
---|
| 28 | return mp.sin(x)/x |
---|
| 29 | |
---|
| 30 | SLD = 3 |
---|
| 31 | SLD_SOLVENT = 6 |
---|
[4f611f1] | 32 | CONTRAST = SLD - SLD_SOLVENT |
---|
| 33 | |
---|
[31eea1f] | 34 | def make_parallelepiped(a, b, c): |
---|
| 35 | def Fq(qa, qb, qc): |
---|
| 36 | siA = sas_sinx_x(0.5*a*qa) |
---|
| 37 | siB = sas_sinx_x(0.5*b*qb) |
---|
| 38 | siC = sas_sinx_x(0.5*c*qc) |
---|
| 39 | return siA * siB * siC |
---|
| 40 | volume = a*b*c |
---|
| 41 | norm = volume*CONTRAST**2/10**4 |
---|
| 42 | return norm, Fq |
---|
| 43 | |
---|
| 44 | def make_parallelepiped_mp(a, b, c): |
---|
| 45 | a, b, c = mp.mpf(a), mp.mpf(b), mp.mpf(c) |
---|
| 46 | def Fq(qa, qb, qc): |
---|
| 47 | siA = mp_sinx_x(a*qa/2) |
---|
| 48 | siB = mp_sinx_x(b*qb/2) |
---|
| 49 | siC = mp_sinx_x(c*qc/2) |
---|
| 50 | return siA * siB * siC |
---|
| 51 | volume = a*b*c |
---|
| 52 | norm = (volume*CONTRAST**2)/10000 # mpf since volume=a*b*c is mpf |
---|
| 53 | return norm, Fq |
---|
| 54 | |
---|
| 55 | def make_triellip(a, b, c): |
---|
| 56 | def Fq(qa, qb, qc): |
---|
| 57 | qr = sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2) |
---|
| 58 | return sas_3j1x_x(qr) |
---|
| 59 | volume = 4*pi*a*b*c/3 |
---|
| 60 | norm = volume*CONTRAST**2/10**4 |
---|
| 61 | return norm, Fq |
---|
| 62 | |
---|
| 63 | def make_triellip_mp(a, b, c): |
---|
| 64 | a, b, c = mp.mpf(a), mp.mpf(b), mp.mpf(c) |
---|
| 65 | def Fq(qa, qb, qc): |
---|
| 66 | qr = mp.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2) |
---|
| 67 | return mp_3j1x_x(qr) |
---|
| 68 | volume = (4*mp.pi*a*b*c)/3 |
---|
| 69 | norm = (volume*CONTRAST**2)/10000 # mpf since mp.pi is mpf |
---|
| 70 | return norm, Fq |
---|
| 71 | |
---|
[4f611f1] | 72 | def make_cylinder(radius, length): |
---|
[31eea1f] | 73 | def Fq(qa, qb, qc): |
---|
[4f611f1] | 74 | qab = sqrt(qa**2 + qb**2) |
---|
| 75 | return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length) |
---|
| 76 | volume = pi*radius**2*length |
---|
[31eea1f] | 77 | norm = volume*CONTRAST**2/10**4 |
---|
| 78 | return norm, Fq |
---|
| 79 | |
---|
| 80 | def make_cylinder_mp(radius, length): |
---|
| 81 | radius, length = mp.mpf(radius), mp.mpf(length) |
---|
| 82 | def Fq(qa, qb, qc): |
---|
| 83 | qab = mp.sqrt(qa**2 + qb**2) |
---|
| 84 | return mp_2J1x_x(qab*radius) * mp_sinx_x((qc*length)/2) |
---|
| 85 | volume = mp.pi*radius**2*length |
---|
| 86 | norm = (volume*CONTRAST**2)/10000 # mpf since mp.pi is mpf |
---|
| 87 | return norm, Fq |
---|
[4f611f1] | 88 | |
---|
| 89 | def make_sphere(radius): |
---|
[31eea1f] | 90 | def Fq(qa, qb, qc): |
---|
| 91 | q = sqrt(qa**2 + qb**2 + qc**2) |
---|
[4f611f1] | 92 | return sas_3j1x_x(q*radius) |
---|
| 93 | volume = 4*pi*radius**3/3 |
---|
[31eea1f] | 94 | norm = volume*CONTRAST**2/10**4 |
---|
| 95 | return norm, Fq |
---|
| 96 | |
---|
| 97 | def make_sphere_mp(radius): |
---|
| 98 | radius = mp.mpf(radius) |
---|
| 99 | def Fq(qa, qb, qc): |
---|
| 100 | q = mp.sqrt(qa**2 + qb**2 + qc**2) |
---|
| 101 | return mp_3j1x_x(q*radius) |
---|
| 102 | volume = (4*mp.pi*radius**3)/3 |
---|
| 103 | norm = (volume*CONTRAST**2)/10000 # mpf since mp.pi is mpf |
---|
| 104 | return norm, Fq |
---|
| 105 | |
---|
| 106 | shape = 'parallelepiped' |
---|
| 107 | #shape = 'triellip' |
---|
| 108 | #shape = 'sphere' |
---|
| 109 | #shape = 'cylinder' |
---|
| 110 | if shape == 'cylinder': |
---|
| 111 | #RADIUS, LENGTH = 10, 100000 |
---|
| 112 | RADIUS, LENGTH = 10, 300 # integer for the sake of mpf |
---|
| 113 | NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) |
---|
| 114 | NORM_MP, KERNEL_MP = make_cylinder_mp(radius=RADIUS, length=LENGTH) |
---|
| 115 | elif shape == 'triellip': |
---|
| 116 | #A, B, C = 4450, 14000, 47 |
---|
| 117 | A, B, C = 445, 140, 47 # integer for the sake of mpf |
---|
| 118 | NORM, KERNEL = make_triellip(A, B, C) |
---|
| 119 | NORM_MP, KERNEL_MP = make_triellip_mp(A, B, C) |
---|
| 120 | elif shape == 'parallelepiped': |
---|
| 121 | #A, B, C = 4450, 14000, 47 |
---|
| 122 | A, B, C = 445, 140, 47 # integer for the sake of mpf |
---|
| 123 | NORM, KERNEL = make_parallelepiped(A, B, C) |
---|
| 124 | NORM_MP, KERNEL_MP = make_parallelepiped_mp(A, B, C) |
---|
| 125 | elif shape == 'sphere': |
---|
| 126 | RADIUS = 50 # integer for the sake of mpf |
---|
| 127 | NORM, KERNEL = make_sphere(radius=RADIUS) |
---|
| 128 | NORM_MP, KERNEL_MP = make_sphere_mp(radius=RADIUS) |
---|
| 129 | else: |
---|
| 130 | raise ValueError("Unknown shape %r"%shape) |
---|
[4f611f1] | 131 | |
---|
| 132 | THETA_LOW, THETA_HIGH = 0, pi |
---|
| 133 | PHI_LOW, PHI_HIGH = 0, 2*pi |
---|
| 134 | SCALE = 1 |
---|
| 135 | |
---|
[31eea1f] | 136 | # mathematica code for triaxial_ellipsoid (untested) |
---|
| 137 | _ = """ |
---|
| 138 | R[theta_, phi_, a_, b_, c_] := Sqrt[(a Sin[theta]Cos[phi])^2 + (b Sin[theta]Sin[phi])^2 + (c Cos[theta])^2] |
---|
| 139 | Sphere[q_, r_] := 3 SphericalBesselJ[q r]/(q r) |
---|
| 140 | V[a_, b_, c_] := 4/3 pi a b c |
---|
| 141 | Norm[sld_, solvent_, a_, b_, c_] := V[a, b, c] (solvent - sld)^2 |
---|
| 142 | F[q_, theta_, phi_, a_, b_, c_] := Sphere[q, R[theta, phi, a, b, c]] |
---|
| 143 | I[q_, sld_, solvent_, a_, b_, c_] := Norm[sld, solvent, a, b, c]/(4 pi) Integrate[F[q, theta, phi, a, b, c]^2 Sin[theta], {phi, 0, 2 pi}, {theta, 0, pi}] |
---|
| 144 | I[6/10^3, 63/10, 3, 445, 140, 47] |
---|
| 145 | """ |
---|
| 146 | |
---|
| 147 | |
---|
| 148 | def mp_quad(q, shape): |
---|
| 149 | evals = [0] |
---|
| 150 | def integrand(theta, phi): |
---|
| 151 | evals[0] += 1 |
---|
| 152 | qab = q*mp.sin(theta) |
---|
| 153 | qa = qab*mp.cos(phi) |
---|
| 154 | qb = qab*mp.sin(phi) |
---|
| 155 | qc = q*mp.cos(theta) |
---|
| 156 | Zq = KERNEL_MP(qa, qb, qc)**2 |
---|
| 157 | return Zq*mp.sin(theta) |
---|
| 158 | ans = mp.quad(integrand, (0, mp.pi), (0, 2*mp.pi)) |
---|
| 159 | Iq = NORM_MP*ans/(4*mp.pi) |
---|
| 160 | return evals[0], Iq |
---|
[4f611f1] | 161 | |
---|
| 162 | def kernel(q, theta, phi): |
---|
| 163 | """ |
---|
| 164 | S(q) kernel for paracrystal forms. |
---|
| 165 | """ |
---|
| 166 | qab = q*sin(theta) |
---|
| 167 | qa = qab*cos(phi) |
---|
| 168 | qb = qab*sin(phi) |
---|
| 169 | qc = q*cos(theta) |
---|
| 170 | return NORM*KERNEL(qa, qb, qc)**2 |
---|
| 171 | |
---|
| 172 | def scipy_dblquad(q): |
---|
| 173 | """ |
---|
| 174 | Compute the integral using scipy dblquad. This gets the correct answer |
---|
| 175 | eventually, but it is slow. |
---|
| 176 | """ |
---|
| 177 | evals = [0] |
---|
[31eea1f] | 178 | def integrand(phi, theta): |
---|
[4f611f1] | 179 | evals[0] += 1 |
---|
| 180 | Zq = kernel(q, theta=theta, phi=phi) |
---|
| 181 | return Zq*sin(theta) |
---|
[31eea1f] | 182 | ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0] |
---|
| 183 | return evals[0], ans*SCALE/(4*pi) |
---|
[4f611f1] | 184 | |
---|
| 185 | |
---|
| 186 | def scipy_romberg_2d(q): |
---|
| 187 | """ |
---|
| 188 | Compute the integral using romberg integration. This function does not |
---|
| 189 | complete in a reasonable time. No idea if it is accurate. |
---|
| 190 | """ |
---|
[31eea1f] | 191 | evals = [0] |
---|
[4f611f1] | 192 | def inner(phi, theta): |
---|
[31eea1f] | 193 | evals[0] += 1 |
---|
[4f611f1] | 194 | return kernel(q, theta=theta, phi=phi) |
---|
| 195 | def outer(theta): |
---|
[31eea1f] | 196 | Zq = romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,)) |
---|
| 197 | return Zq*sin(theta) |
---|
| 198 | ans = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100) |
---|
| 199 | return evals[0], ans*SCALE/(4*pi) |
---|
[4f611f1] | 200 | |
---|
| 201 | |
---|
| 202 | def semi_romberg(q, n=100): |
---|
| 203 | """ |
---|
| 204 | Use 1D romberg integration in phi and regular simpsons rule in theta. |
---|
| 205 | """ |
---|
| 206 | evals = [0] |
---|
| 207 | def inner(phi, theta): |
---|
| 208 | evals[0] += 1 |
---|
| 209 | return kernel(q, theta=theta, phi=phi) |
---|
| 210 | theta = np.linspace(THETA_LOW, THETA_HIGH, n) |
---|
[31eea1f] | 211 | Zq = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,)) for t in theta] |
---|
| 212 | ans = simps(np.array(Zq)*sin(theta), dx=theta[1]-theta[0]) |
---|
| 213 | return evals[0], ans*SCALE/(4*pi) |
---|
[4f611f1] | 214 | |
---|
| 215 | def gauss_quad(q, n=150): |
---|
| 216 | """ |
---|
| 217 | Compute the integral using gaussian quadrature for n = 20, 76 or 150. |
---|
| 218 | """ |
---|
| 219 | if n == 20: |
---|
| 220 | z, w = Gauss20Z, Gauss20Wt |
---|
| 221 | elif n == 76: |
---|
| 222 | z, w = Gauss76Z, Gauss76Wt |
---|
[31eea1f] | 223 | elif n == 150: |
---|
[4f611f1] | 224 | z, w = Gauss150Z, Gauss150Wt |
---|
[31eea1f] | 225 | else: |
---|
| 226 | z, w = leggauss(n) |
---|
[4f611f1] | 227 | theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW |
---|
| 228 | phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW |
---|
| 229 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
| 230 | Aw = w[None, :] * w[:, None] |
---|
[31eea1f] | 231 | sin_theta = abs(sin(Atheta)) |
---|
[4f611f1] | 232 | Zq = kernel(q=q, theta=Atheta, phi=Aphi) |
---|
[31eea1f] | 233 | # change from [-1,1] x [-1,1] range to [0, pi] x [0, 2 pi] range |
---|
| 234 | dxdy_stretch = (THETA_HIGH-THETA_LOW)/2 * (PHI_HIGH-PHI_LOW)/2 |
---|
| 235 | Iq = np.sum(Zq*Aw*sin_theta)*SCALE/(4*pi) * dxdy_stretch |
---|
| 236 | return n**2, Iq |
---|
[4f611f1] | 237 | |
---|
| 238 | |
---|
| 239 | def gridded_integrals(q, n=300): |
---|
| 240 | """ |
---|
| 241 | Compute the integral on a regular grid using rectangular, trapezoidal, |
---|
| 242 | simpsons, and romberg integration. Romberg integration requires that |
---|
| 243 | the grid be of size n = 2**k + 1. |
---|
| 244 | """ |
---|
| 245 | theta = np.linspace(THETA_LOW, THETA_HIGH, n) |
---|
| 246 | phi = np.linspace(PHI_LOW, PHI_HIGH, n) |
---|
| 247 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
| 248 | Zq = kernel(q=q, theta=Atheta, phi=Aphi) |
---|
| 249 | Zq *= abs(sin(Atheta)) |
---|
| 250 | dx, dy = theta[1]-theta[0], phi[1]-phi[0] |
---|
[31eea1f] | 251 | print("rect-%d"%n, n**2, np.sum(Zq)*dx*dy*SCALE/(4*pi)) |
---|
| 252 | print("trapz-%d"%n, n**2, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) |
---|
| 253 | print("simpson-%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) |
---|
| 254 | print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) |
---|
[4f611f1] | 255 | |
---|
| 256 | def plot(q, n=300): |
---|
| 257 | """ |
---|
| 258 | Plot the 2D surface that needs to be integrated in order to compute |
---|
| 259 | the BCC S(q) at a particular q, dnn and d_factor. *n* is the number |
---|
| 260 | of points in the grid. |
---|
| 261 | """ |
---|
| 262 | theta = np.linspace(THETA_LOW, THETA_HIGH, n) |
---|
| 263 | phi = np.linspace(PHI_LOW, PHI_HIGH, n) |
---|
| 264 | Atheta, Aphi = np.meshgrid(theta, phi) |
---|
| 265 | Zq = kernel(q=q, theta=Atheta, phi=Aphi) |
---|
| 266 | #Zq *= abs(sin(Atheta)) |
---|
| 267 | pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6))) |
---|
| 268 | pylab.axis('tight') |
---|
| 269 | pylab.title("%s Z(q) for q=%g" % (KERNEL.__name__, q)) |
---|
| 270 | pylab.xlabel("theta (degrees)") |
---|
| 271 | pylab.ylabel("phi (degrees)") |
---|
| 272 | cbar = pylab.colorbar() |
---|
| 273 | cbar.set_label('log10 S(q)') |
---|
| 274 | pylab.show() |
---|
| 275 | |
---|
| 276 | if __name__ == "__main__": |
---|
[31eea1f] | 277 | Qstr = '0.005' |
---|
| 278 | #Qstr = '0.8' |
---|
| 279 | #Qstr = '0.0003' |
---|
| 280 | Q = float(Qstr) |
---|
| 281 | if shape == 'sphere': |
---|
| 282 | print("exact", NORM*sas_3j1x_x(Q*RADIUS)**2) |
---|
| 283 | print("gauss-20", *gauss_quad(Q, n=20)) |
---|
| 284 | print("gauss-76", *gauss_quad(Q, n=76)) |
---|
| 285 | print("gauss-150", *gauss_quad(Q, n=150)) |
---|
| 286 | print("gauss-500", *gauss_quad(Q, n=500)) |
---|
| 287 | print("dblquad", *scipy_dblquad(Q)) |
---|
| 288 | print("semi-romberg-100", *semi_romberg(Q, n=100)) |
---|
| 289 | print("romberg", *scipy_romberg_2d(Q)) |
---|
[4f611f1] | 290 | #gridded_integrals(Q, n=2**8+1) |
---|
[31eea1f] | 291 | gridded_integrals(Q, n=2**10+1) |
---|
[4f611f1] | 292 | #gridded_integrals(Q, n=2**13+1) |
---|
[31eea1f] | 293 | #gridded_integrals(Q, n=2**15+1) |
---|
| 294 | with mp.workprec(100): |
---|
| 295 | print("mpmath", *mp_quad(mp.mpf(Qstr), shape)) |
---|
| 296 | #plot(Q, n=200) |
---|