source: sasmodels/explore/asymint.py @ 66dbbfb

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 66dbbfb was 66dbbfb, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

prettier output for explore/asymint.py

  • Property mode set to 100755
File size: 20.1 KB
RevLine 
[20fe0cd]1#!/usr/bin/env python
[4f611f1]2"""
3Asymmetric shape integration
[5110e16]4
5Usage:
6
7    explore/asymint.py [MODEL] [q-value]
8
9Computes the numerical integral over theta and phi of the given model at a
10single point q using different algorithms or the same algorithm with different
11precision.  It also displays a 2-D image of the theta-phi surface that is
12being integrated.
13
14The available models are:
15
16    triaxial_ellipsoid, parallelpiped, paracrystal, cylinder, sphere
17
18Cylinder and sphere are included as simple checks on the integration
19algorithms. Cylinder is better investigated using 1-D integration methods in
20explore/symint.py.  Sphere has an easily computed analytic value which is
21identical for all theta-phi for a given q, so it is useful for checking
22that the normalization constants are correct for the different algorithms.
[4f611f1]23"""
24
25from __future__ import print_function, division
26
27import os, sys
[20fe0cd]28sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
[242b361]29import warnings
[4f611f1]30
31import numpy as np
[31eea1f]32import mpmath as mp
[c462169]33from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10, arccos
[31eea1f]34from numpy.polynomial.legendre import leggauss
[4f611f1]35from scipy.integrate import dblquad, simps, romb, romberg
36import pylab
37
[20fe0cd]38import sasmodels.special as sp
[4f611f1]39
[c462169]40DTYPE = 'd'
41
[20fe0cd]42class MPenv:
[6e604f8]43    sqrt = staticmethod(mp.sqrt)
44    exp = staticmethod(mp.exp)
45    expm1 = staticmethod(mp.expm1)
46    cos = staticmethod(mp.cos)
47    sin = staticmethod(mp.sin)
48    tan = staticmethod(mp.tan)
49    @staticmethod
[20fe0cd]50    def sas_3j1x_x(x):
51        return 3*(mp.sin(x)/x - mp.cos(x))/(x*x)
[6e604f8]52    @staticmethod
[20fe0cd]53    def sas_2J1x_x(x):
54        return 2*mp.j1(x)/x
[6e604f8]55    @staticmethod
[20fe0cd]56    def sas_sinx_x(x):
57        return mp.sin(x)/x
58    pi = mp.pi
[6e604f8]59    mpf = staticmethod(mp.mpf)
[20fe0cd]60
61class NPenv:
[6e604f8]62    sqrt = staticmethod(np.sqrt)
63    exp = staticmethod(np.exp)
64    expm1 = staticmethod(np.expm1)
65    cos = staticmethod(np.cos)
66    sin = staticmethod(np.sin)
67    tan = staticmethod(np.tan)
68    sas_3j1x_x = staticmethod(sp.sas_3j1x_x)
69    sas_2J1x_x = staticmethod(sp.sas_2J1x_x)
70    sas_sinx_x = staticmethod(sp.sas_sinx_x)
[20fe0cd]71    pi = np.pi
[c462169]72    #mpf = staticmethod(float)
73    mpf = staticmethod(lambda x: np.array(x, DTYPE))
[31eea1f]74
75SLD = 3
76SLD_SOLVENT = 6
[4f611f1]77CONTRAST = SLD - SLD_SOLVENT
78
[20fe0cd]79# Carefully code models so that mpmath will use full precision.  That means:
80#    * wrap inputs in env.mpf
81#    * don't use floating point constants, only integers
82#    * for division, make sure the numerator or denominator is env.mpf
83#    * use env.pi, env.sas_sinx_x, etc. for functions
84def make_parallelepiped(a, b, c, env=NPenv):
85    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
[31eea1f]86    def Fq(qa, qb, qc):
[49eb251]87        siA = env.sas_sinx_x(a*qa/2)
88        siB = env.sas_sinx_x(b*qb/2)
89        siC = env.sas_sinx_x(c*qc/2)
[31eea1f]90        return siA * siB * siC
[20fe0cd]91    Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c)
[31eea1f]92    volume = a*b*c
[20fe0cd]93    norm = CONTRAST**2*volume/10000
[31eea1f]94    return norm, Fq
95
[49eb251]96def make_core_shell_parallelepiped(a, b, c, da, db, dc, slda, sldb, sldc, env=NPenv):
[a189283]97    overlapping = False
[49eb251]98    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
99    da, db, dc = env.mpf(da), env.mpf(db), env.mpf(dc)
100    slda, sldb, sldc = env.mpf(slda), env.mpf(sldb), env.mpf(sldc)
[a189283]101    dr0 = CONTRAST
102    drA, drB, drC = slda-SLD_SOLVENT, sldb-SLD_SOLVENT, sldc-SLD_SOLVENT
103    tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc
[49eb251]104    def Fq(qa, qb, qc):
[a189283]105        siA = a*env.sas_sinx_x(a*qa/2)
106        siB = b*env.sas_sinx_x(b*qb/2)
107        siC = c*env.sas_sinx_x(c*qc/2)
108        siAt = tA*env.sas_sinx_x(tA*qa/2)
109        siBt = tB*env.sas_sinx_x(tB*qb/2)
110        siCt = tC*env.sas_sinx_x(tC*qc/2)
111        if overlapping:
112            return (dr0*siA*siB*siC
113                    + drA*(siAt-siA)*siB*siC
114                    + drB*siAt*(siBt-siB)*siC
115                    + drC*siAt*siBt*(siCt-siC))
116        else:
117            return (dr0*siA*siB*siC
118                    + drA*(siAt-siA)*siB*siC
119                    + drB*siA*(siBt-siB)*siC
120                    + drC*siA*siB*(siCt-siC))
[49eb251]121    Fq.__doc__ = "core-shell parallelepiped a=%g, b=%g c=%g"%(a, b, c)
[a189283]122    if overlapping:
123        volume = a*b*c + 2*da*b*c + 2*tA*db*c + 2*tA*tB*dc
124    else:
125        volume = a*b*c + 2*da*b*c + 2*a*db*c + 2*a*b*dc
[49eb251]126    norm = 1/(volume*10000)
127    return norm, Fq
128
[20fe0cd]129def make_triaxial_ellipsoid(a, b, c, env=NPenv):
130    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c)
[31eea1f]131    def Fq(qa, qb, qc):
[20fe0cd]132        qr = env.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2)
133        return env.sas_3j1x_x(qr)
[1820208]134    Fq.__doc__ = "triaxial ellipsoid minor=%g, major=%g polar=%g"%(a, b, c)
[20fe0cd]135    volume = 4*env.pi*a*b*c/3
136    norm = CONTRAST**2*volume/10000
[31eea1f]137    return norm, Fq
138
[20fe0cd]139def make_cylinder(radius, length, env=NPenv):
140    radius, length = env.mpf(radius), env.mpf(length)
[31eea1f]141    def Fq(qa, qb, qc):
[20fe0cd]142        qab = env.sqrt(qa**2 + qb**2)
143        return env.sas_2J1x_x(qab*radius) * env.sas_sinx_x((qc*length)/2)
144    Fq.__doc__ = "cylinder radius=%g, length=%g"%(radius, length)
145    volume = env.pi*radius**2*length
146    norm = CONTRAST**2*volume/10000
[31eea1f]147    return norm, Fq
148
[20fe0cd]149def make_sphere(radius, env=NPenv):
150    radius = env.mpf(radius)
[31eea1f]151    def Fq(qa, qb, qc):
[20fe0cd]152        q = env.sqrt(qa**2 + qb**2 + qc**2)
153        return env.sas_3j1x_x(q*radius)
154    Fq.__doc__ = "sphere radius=%g"%(radius, )
155    volume = 4*pi*radius**3
156    norm = CONTRAST**2*volume/10000
[31eea1f]157    return norm, Fq
158
[20fe0cd]159def make_paracrystal(radius, dnn, d_factor, lattice='bcc', env=NPenv):
160    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor)
161    def sc(qa, qb, qc):
162        return qa, qb, qc
163    def bcc(qa, qb, qc):
164        a1 = (+qa + qb + qc)/2
165        a2 = (-qa - qb + qc)/2
166        a3 = (-qa + qb - qc)/2
167        return a1, a2, a3
168    def fcc(qa, qb, qc):
[5110e16]169        a1 = ( 0 + qb + qc)/2
170        a2 = (-qa + 0 + qc)/2
171        a3 = (-qa + qb + 0)/2
[20fe0cd]172        return a1, a2, a3
173    lattice_fn = {'sc': sc, 'bcc': bcc, 'fcc': fcc}[lattice]
174    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor)
[31eea1f]175    def Fq(qa, qb, qc):
[20fe0cd]176        a1, a2, a3 = lattice_fn(qa, qb, qc)
177        # Note: paper says that different directions can have different
178        # distoration factors.  Easy enough to add to the code.
179        arg = -(dnn*d_factor)**2*(a1**2 + a2**2 + a3**2)/2
180        exp_arg = env.exp(arg)
181        den = [((exp_arg - 2*env.cos(dnn*a))*exp_arg + 1) for a in (a1, a2, a3)]
182        Sq = -env.expm1(2*arg)**3/(den[0]*den[1]*den[2])
[31eea1f]183
[20fe0cd]184        q = env.sqrt(qa**2 + qb**2 + qc**2)
185        Fq = env.sas_3j1x_x(q*radius)
[5110e16]186        # the caller computes F(q)**2, but we need it to compute S(q)*F(q)**2
[20fe0cd]187        return env.sqrt(Sq)*Fq
188    Fq.__doc__ = "%s paracrystal a=%g da=%g r=%g"%(lattice, dnn, d_factor, radius)
189    def sphere_volume(r): return 4*env.pi*r**3/3
190    Vf = {
191        'sc': sphere_volume(radius/dnn),
192        'bcc': 2*sphere_volume(env.sqrt(3)/2*radius/dnn),
[5110e16]193        'fcc': 4*sphere_volume(1/env.sqrt(2)*radius/dnn),
[20fe0cd]194    }[lattice]
195    volume = sphere_volume(radius)
[5110e16]196    norm = CONTRAST**2*volume/10000*Vf
[31eea1f]197    return norm, Fq
[4f611f1]198
[242b361]199NORM = 1.0  # type: float
200KERNEL = None  # type: CALLABLE[[ndarray, ndarray, ndarray], ndarray]
201NORM_MP = 1  # type: mpf
202KERNEL = None  # type: CALLABLE[[mpf, mpf, mpf], mpf]
203
204SHAPES = [
205    'sphere',
206    'cylinder',
207    'triaxial_ellipsoid',
208    'parallelepiped',
209    'core_shell_parallelepiped',
210    'fcc_paracrystal',
211    'bcc_paracrystal',
212    'sc_paracrystal',
213]
214def build_shape(shape, **pars):
215    global NORM, KERNEL
216    global NORM_MP, KERNEL_MP
217
218    # Note: using integer or string defaults for the sake of mpf
219    if shape == 'sphere':
220        RADIUS = pars.get('radius', 50)
221        NORM, KERNEL = make_sphere(radius=RADIUS)
222        NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv)
223    elif shape == 'cylinder':
224        #RADIUS, LENGTH = 10, 100000
225        RADIUS = pars.get('radius', 10)
226        LENGTH = pars.get('radius', 300)
227        NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH)
228        NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv)
229    elif shape == 'triaxial_ellipsoid':
230        #A, B, C = 4450, 14000, 47
231        A = pars.get('a', 445)
232        B = pars.get('b', 140)
233        C = pars.get('c', 47)
234        NORM, KERNEL = make_triaxial_ellipsoid(A, B, C)
235        NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv)
236    elif shape == 'parallelepiped':
237        #A, B, C = 4450, 14000, 47
238        A = pars.get('a', 445)
239        B = pars.get('b', 140)
240        C = pars.get('c', 47)
241        NORM, KERNEL = make_parallelepiped(A, B, C)
242        NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv)
243    elif shape == 'core_shell_parallelepiped':
244        #A, B, C = 4450, 14000, 47
245        #A, B, C = 445, 140, 47  # integer for the sake of mpf
246        A = pars.get('a', 114)
247        B = pars.get('b', 1380)
248        C = pars.get('c', 6800)
249        DA = pars.get('da', 21)
250        DB = pars.get('db', 58)
251        DC = pars.get('dc', 2300)
252        SLDA = pars.get('slda', "5")
253        SLDB = pars.get('sldb', "-0.3")
254        SLDC = pars.get('sldc', "11.5")
255        ## default parameters from sasmodels
256        #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2
257        ## swap A-B-C to C-B-A
258        #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA
259        #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3
260        #SLD_SOLVENT,CONTRAST = 0, 4
261        if 1: # C shortest
262            B, C = C, B
263            DB, DC = DC, DB
264            SLDB, SLDC = SLDC, SLDB
265        elif 0: # C longest
266            A, C = C, A
267            DA, DC = DC, DA
268            SLDA, SLDC = SLDC, SLDA
269        #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
270        NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC)
271        NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv)
272    elif shape.endswith('paracrystal'):
273        LATTICE, _ = shape.split('_')
274        DNN = pars.get('dnn', 220)
275        D_FACTOR = pars.get('d_factor', '0.06')
276        RADIUS = pars.get('radius', 40)
277        NORM, KERNEL = make_paracrystal(
278            radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE)
279        NORM_MP, KERNEL_MP = make_paracrystal(
280            radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv)
281    else:
282        raise ValueError("Unknown shape %r"%shape)
[4f611f1]283
[20fe0cd]284# Note: hardcoded in mp_quad
[4f611f1]285THETA_LOW, THETA_HIGH = 0, pi
286PHI_LOW, PHI_HIGH = 0, 2*pi
287SCALE = 1
288
[31eea1f]289# mathematica code for triaxial_ellipsoid (untested)
290_ = """
291R[theta_, phi_, a_, b_, c_] := Sqrt[(a Sin[theta]Cos[phi])^2 + (b Sin[theta]Sin[phi])^2 + (c Cos[theta])^2]
292Sphere[q_, r_] := 3 SphericalBesselJ[q r]/(q r)
293V[a_, b_, c_] := 4/3 pi a b c
294Norm[sld_, solvent_, a_, b_, c_] := V[a, b, c] (solvent - sld)^2
295F[q_, theta_, phi_, a_, b_, c_] := Sphere[q, R[theta, phi, a, b, c]]
296I[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}]
297I[6/10^3, 63/10, 3, 445, 140, 47]
298"""
299
[20fe0cd]300# 2D integration functions
[242b361]301def mp_quad_2d(q):
[31eea1f]302    evals = [0]
303    def integrand(theta, phi):
304        evals[0] += 1
305        qab = q*mp.sin(theta)
306        qa = qab*mp.cos(phi)
307        qb = qab*mp.sin(phi)
308        qc = q*mp.cos(theta)
309        Zq = KERNEL_MP(qa, qb, qc)**2
310        return Zq*mp.sin(theta)
311    ans = mp.quad(integrand, (0, mp.pi), (0, 2*mp.pi))
312    Iq = NORM_MP*ans/(4*mp.pi)
313    return evals[0], Iq
[4f611f1]314
[20fe0cd]315def kernel_2d(q, theta, phi):
[4f611f1]316    """
317    S(q) kernel for paracrystal forms.
318    """
319    qab = q*sin(theta)
320    qa = qab*cos(phi)
321    qb = qab*sin(phi)
322    qc = q*cos(theta)
323    return NORM*KERNEL(qa, qb, qc)**2
324
[20fe0cd]325def scipy_dblquad_2d(q):
[4f611f1]326    """
327    Compute the integral using scipy dblquad.  This gets the correct answer
328    eventually, but it is slow.
329    """
330    evals = [0]
[31eea1f]331    def integrand(phi, theta):
[4f611f1]332        evals[0] += 1
[20fe0cd]333        Zq = kernel_2d(q, theta=theta, phi=phi)
[4f611f1]334        return Zq*sin(theta)
[31eea1f]335    ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0]
336    return evals[0], ans*SCALE/(4*pi)
[4f611f1]337
338
339def scipy_romberg_2d(q):
340    """
341    Compute the integral using romberg integration.  This function does not
342    complete in a reasonable time.  No idea if it is accurate.
343    """
[31eea1f]344    evals = [0]
[4f611f1]345    def inner(phi, theta):
[31eea1f]346        evals[0] += 1
[20fe0cd]347        return kernel_2d(q, theta=theta, phi=phi)
[4f611f1]348    def outer(theta):
[31eea1f]349        Zq = romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,))
350        return Zq*sin(theta)
351    ans = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)
352    return evals[0], ans*SCALE/(4*pi)
[4f611f1]353
354
[20fe0cd]355def semi_romberg_2d(q, n=100):
[4f611f1]356    """
357    Use 1D romberg integration in phi and regular simpsons rule in theta.
358    """
359    evals = [0]
360    def inner(phi, theta):
361        evals[0] += 1
[20fe0cd]362        return kernel_2d(q, theta=theta, phi=phi)
[4f611f1]363    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
[31eea1f]364    Zq = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,)) for t in theta]
365    ans = simps(np.array(Zq)*sin(theta), dx=theta[1]-theta[0])
366    return evals[0], ans*SCALE/(4*pi)
[4f611f1]367
[20fe0cd]368def gauss_quad_2d(q, n=150):
[4f611f1]369    """
370    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
371    """
[20fe0cd]372    z, w = leggauss(n)
[4f611f1]373    theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW
374    phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW
375    Atheta, Aphi = np.meshgrid(theta, phi)
376    Aw = w[None, :] * w[:, None]
[31eea1f]377    sin_theta = abs(sin(Atheta))
[20fe0cd]378    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
[31eea1f]379    # change from [-1,1] x [-1,1] range to [0, pi] x [0, 2 pi] range
380    dxdy_stretch = (THETA_HIGH-THETA_LOW)/2 * (PHI_HIGH-PHI_LOW)/2
381    Iq = np.sum(Zq*Aw*sin_theta)*SCALE/(4*pi) * dxdy_stretch
382    return n**2, Iq
[4f611f1]383
[c462169]384def gauss_quad_usub(q, n=150, dtype=DTYPE):
385    """
386    Compute the integral using gaussian quadrature for n = 20, 76 or 150.
387
388    Use *u = sin theta* substitution, and restrict integration over a single
389    quadrant for shapes that are mirror symmetric about AB, AC and BC planes.
390
391    Note that this doesn't work for fcc/bcc paracrystals, which instead step
392    over the entire 4 pi surface uniformly in theta-phi.
393    """
394    z, w = leggauss(n)
395    cos_theta = 0.5 * (z + 1)
396    theta = arccos(cos_theta)
397    phi = pi/2*(0.5 * (z + 1))
398    Atheta, Aphi = np.meshgrid(theta, phi)
399    Aw = w[None, :] * w[:, None]
400    q, Atheta, Aphi, Aw = [np.asarray(v, dtype=dtype) for v in (q, Atheta, Aphi, Aw)]
401    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
402    Iq = np.sum(Zq*Aw)*0.25
403    return n**2, Iq
404
[20fe0cd]405def gridded_2d(q, n=300):
[4f611f1]406    """
407    Compute the integral on a regular grid using rectangular, trapezoidal,
408    simpsons, and romberg integration.  Romberg integration requires that
409    the grid be of size n = 2**k + 1.
410    """
411    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
412    phi = np.linspace(PHI_LOW, PHI_HIGH, n)
413    Atheta, Aphi = np.meshgrid(theta, phi)
[20fe0cd]414    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
[4f611f1]415    Zq *= abs(sin(Atheta))
416    dx, dy = theta[1]-theta[0], phi[1]-phi[0]
[31eea1f]417    print("rect-%d"%n, n**2, np.sum(Zq)*dx*dy*SCALE/(4*pi))
418    print("trapz-%d"%n, n**2, np.trapz(np.trapz(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
419    print("simpson-%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
420    print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi))
[4f611f1]421
[242b361]422def quadpy_method(q, rule):
423    """
424    Use *rule*="name:index" where name and index are chosen from below.
425
426    Available rule names and the corresponding indices::
427
428        AlbrechtCollatz: [1-5]
429        BazantOh: 9, 11, 13
430        HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3],
431            29, 31, 33, 35, 37, 39-[1-2]
432        FliegeMaier: 4, 9, 16, 25
433        Lebedev: 3[a-c], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35,
434            41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131
435        McLaren: [1-10]
436        Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3],
437            U3 11-[1-3], U3 14-1
438    """
439    try:
440        import quadpy
441    except ImportError:
442        warnings.warn("use 'pip install quadpy' to enable quadpy.sphere tests")
443        return
444
445    from quadpy.sphere import (AlbrechtCollatz, BazantOh, HeoXu,
446        FliegeMaier, Lebedev, McLaren, Stroud, integrate_spherical)
447    RULES = {
448        'AlbrechtCollatz': AlbrechtCollatz,
449        'BazantOh': BazantOh,
450        'HeoXu': HeoXu,
451        'FliegeMaier': FliegeMaier,
452        'Lebedev': Lebedev,
453        'McLaren': McLaren,
454        'Stroud': Stroud,
455    }
456    int_index = 'AlbrechtCollatz', 'McLaren'
457
458    rule_name, rule_index = rule.split(':')
459    index = int(rule_index) if rule_name in int_index else rule_index
460    rule_obj = RULES[rule_name](index)
461    fn = lambda azimuthal, polar: kernel_2d(q=q, theta=polar, phi=azimuthal)
462    Iq = integrate_spherical(fn, rule=rule_obj)/(4*pi)
[66dbbfb]463    print("%s degree=%d points=%s => %.15g"
[242b361]464          % (rule, rule_obj.degree, len(rule_obj.points), Iq))
465
[20fe0cd]466def plot_2d(q, n=300):
[4f611f1]467    """
468    Plot the 2D surface that needs to be integrated in order to compute
469    the BCC S(q) at a particular q, dnn and d_factor.  *n* is the number
470    of points in the grid.
471    """
472    theta = np.linspace(THETA_LOW, THETA_HIGH, n)
473    phi = np.linspace(PHI_LOW, PHI_HIGH, n)
474    Atheta, Aphi = np.meshgrid(theta, phi)
[20fe0cd]475    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi)
[4f611f1]476    #Zq *= abs(sin(Atheta))
477    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6)))
478    pylab.axis('tight')
[20fe0cd]479    pylab.title("%s I(q,t) sin(t) for q=%g" % (KERNEL.__doc__, q))
[4f611f1]480    pylab.xlabel("theta (degrees)")
481    pylab.ylabel("phi (degrees)")
482    cbar = pylab.colorbar()
483    cbar.set_label('log10 S(q)')
484    pylab.show()
485
[242b361]486def main():
487    import argparse
488
489    parser = argparse.ArgumentParser(
490        description="asymmetric integration explorer",
491        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
492        )
493    parser.add_argument('-s', '--shape', choices=SHAPES,
494                        default='parallelepiped',
495                        help='oriented shape')
496    parser.add_argument('-q', '--q_value', type=str, default='0.005',
497                        help='Q value to evaluate')
498    parser.add_argument('pars', type=str, nargs='*', default=[],
499                        help='p=val for p in shape parameters')
500    opts = parser.parse_args()
501    pars = {k: v for par in opts.pars for k, v in [par.split('=')]}
502    build_shape(opts.shape, **pars)
503
504    Q = float(opts.q_value)
505    if opts.shape == 'sphere':
[1820208]506        print("exact", NORM*sp.sas_3j1x_x(Q*RADIUS)**2)
[242b361]507
508    # Methods from quadpy, if quadpy is available
509    #  AlbrechtCollatz: [1-5]
510    #  BazantOh: 9, 11, 13
511    #  HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3],
512    #     29, 31, 33, 35, 37, 39-[1-2]
513    #  FliegeMaier: 4, 9, 16, 25
514    #  Lebedev: 3[a-c], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35,
515    #     41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131
516    #  McLaren: [1-10]
517    #  Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3],
518    #     U3 11-[1-3], U3 14-1
519    quadpy_method(Q, "AlbrechtCollatz:5")
520    quadpy_method(Q, "HeoXu:39-2")
521    quadpy_method(Q, "FliegeMaier:25")
522    quadpy_method(Q, "Lebedev:19")
523    quadpy_method(Q, "Lebedev:131")
524    quadpy_method(Q, "McLaren:10")
525    quadpy_method(Q, "Stroud:U3 14-1")
526
[66dbbfb]527    print("gauss-20 points=%d => %.15g" % gauss_quad_2d(Q, n=20))
528    print("gauss-76 points=%d => %.15g" % gauss_quad_2d(Q, n=76))
529    print("gauss-150 points=%d => %.15g" % gauss_quad_2d(Q, n=150))
530    print("gauss-500 points=%d => %.15g" % gauss_quad_2d(Q, n=500))
531    print("gauss-1025 points=%d => %.15g" % gauss_quad_2d(Q, n=1025))
532    print("gauss-2049 points=%d => %.15g" % gauss_quad_2d(Q, n=2049))
533    print("gauss-20 usub points=%d => %.15g" % gauss_quad_usub(Q, n=20))
534    print("gauss-76 usub points=%d => %.15g" % gauss_quad_usub(Q, n=76))
535    print("gauss-150 usub points=%d => %.15g" % gauss_quad_usub(Q, n=150))
[242b361]536
[20fe0cd]537    #gridded_2d(Q, n=2**8+1)
538    gridded_2d(Q, n=2**10+1)
[49eb251]539    #gridded_2d(Q, n=2**12+1)
[20fe0cd]540    #gridded_2d(Q, n=2**15+1)
[242b361]541    # adaptive forms on models for which the calculations are fast enough
542    SLOW_SHAPES = {
543        'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal',
544        'core_shell_parallelepiped',
545    }
546    if opts.shape not in SLOW_SHAPES:
[20fe0cd]547        print("dblquad", *scipy_dblquad_2d(Q))
548        print("semi-romberg-100", *semi_romberg_2d(Q, n=100))
549        print("romberg", *scipy_romberg_2d(Q))
550        with mp.workprec(100):
[242b361]551            print("mpmath", *mp_quad_2d(mp.mpf(opts.q_value)))
[20fe0cd]552    plot_2d(Q, n=200)
553
554if __name__ == "__main__":
[242b361]555    main()
Note: See TracBrowser for help on using the repository browser.