Changeset cb038a2 in sasmodels for explore/bccpy.py


Ignore:
Timestamp:
Apr 17, 2017 12:54:31 PM (8 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
7e0b281
Parents:
859567e
Message:

add some docs to the bcc/fcc/sc explorer

File:
1 edited

Legend:

Unmodified
Added
Removed

  • explore/bccpy.py

    r859567e rcb038a2  
     1""" 
     2The current 1D calculations for BCC paracrystal are very wrong at low q, orders 
     3of magnitude wrong.  The integration fails to capture a very narrow, 
     4very steep ridge. 
     5 
     6Uncomment the plot() line at the bottom of the code to show an image of the 
     7set of S(q, theta, phi) values that must be summed to form the 1D S(q=0.001) 
     8for the BCC paracrystal form.  Note particularly that this is a log scale 
     9image spanning 10 orders of magnitude.  This pattern repeats itself 8 times 
     10over the entire 4 pi surface integral. 
     11 
     12You can explore various integration options by uncommenting more lines. 
     13Adaptive integration using scpy.integrate.dbsquad is very slow.  Romberg 
     14didn't even complete in the time I gave it. Accurate brute force calculation 
     15requires a 4000x4000 grid to get enough precision. 
     16 
     17We may need a specialized integrator for low q which can identify and integrate 
     18the ridges properly. 
     19 
     20This program need sasmodels on the path so it is inserted automatically, 
     21assuming that the explore/bccpy.py is beside sasmodels/special.py in the 
     22source tree.  Run from the sasmodels directory using: 
     23 
     24    python explore/bccpy.py 
     25""" 
     26 
    127from __future__ import print_function, division 
     28 
     29import os, sys 
     30sys.path.insert(0, os.path.dirname(os.path.dirname(__file__))) 
    231 
    332import numpy as np 
    … …  
    1140from sasmodels.special import Gauss150Wt, Gauss150Z 
    1241 
    13 def square(x): 
    14     return x**2 
    15  
    1642Q = 0.001 
    1743DNN = 220. 
    … …  
    2248 
    2349def kernel(q, dnn, d_factor, theta, phi): 
     50    """ 
     51    S(q) kernel for paracrystal forms. 
     52    """ 
    2453    qab = q*sin(theta) 
    2554    qa = qab*cos(phi) 
    … …  
    3564        a3 = -qa + qc + qb 
    3665        dcos = dnn/2 
     66    if 0: # fcc 
     67        a1 = qb + qa 
     68        a2 = qa + qc 
     69        a3 = qb + qc 
     70        dcos = dnn/2 
    3771 
    3872    arg = 0.5*square(dnn*d_factor)*(a1**2 + a2**2 + a3**2) 
    … …  
    4478 
    4579def scipy_dblquad(q=Q, dnn=DNN, d_factor=D_FACTOR): 
     80    """ 
     81    Compute the integral using scipy dblquad.  This gets the correct answer 
     82    eventually, but it is slow. 
     83    """ 
     84    evals = [0] 
    4685    def integrand(theta, phi): 
     86        evals[0] += 1 
    4787        Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi) 
    4888        return Sq*sin(theta) 
    49     return dblquad(integrand, 0, pi/2, lambda x: 0, lambda x: pi/2)[0]*8/(4*pi) 
     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 
    5092 
    5193 
    52 def scipy_romberg(q=Q, dnn=DNN, d_factor=D_FACTOR): 
     94def 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    """ 
    5399    def inner(phi, theta): 
    54100        return kernel(q=q, dnn=dnn, d_factor=d_factor, theta=theta, phi=phi) 
    … …  
    58104 
    59105 
     106def 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 
    60121def gauss_quad(q=Q, dnn=DNN, d_factor=D_FACTOR, n=150): 
     122    """ 
     123    Compute the integral using gaussian quadrature for n = 20, 76 or 150. 
     124    """ 
    61125    if n == 20: 
    62126        z, w = Gauss20Z, Gauss20Wt 
    … …  
    71135    sin_theta = np.fmax(abs(sin(Atheta)), 1e-6) 
    72136    Sq = kernel(q=q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) 
     137    print("gauss %d evals ="%n, n**2) 
    73138    return np.sum(Sq*Aw*sin_theta)*8/(4*pi) 
    74139 
    75140 
    76141def gridded_integrals(q=0.001, dnn=DNN, d_factor=D_FACTOR, n=300): 
     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    """ 
    77147    theta = np.linspace(0, pi/2, n) 
    78148    phi = np.linspace(0, pi/2, n) 
    … …  
    85155    print("simpson", n, simps(simps(Sq, dx=dx), dx=dy)*8/(4*pi)) 
    86156    print("romb", n, romb(romb(Sq, dx=dx), dx=dy)*8/(4*pi)) 
     157    print("gridded %d evals ="%n, n**2) 
    87158 
    88159def plot(q=0.001, dnn=DNN, d_factor=D_FACTOR, n=300): 
    89     #theta = np.linspace(0, pi, n) 
    90     #phi = np.linspace(0, 2*pi, n) 
    91     theta = np.linspace(0, pi/2, n) 
    92     phi = np.linspace(0, pi/2, n) 
     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) 
    93169    Atheta, Aphi = np.meshgrid(theta, phi) 
    94170    Sq = kernel(q=Q, dnn=dnn, d_factor=d_factor, theta=Atheta, phi=Aphi) 
    95171    Sq *= abs(sin(Atheta)) 
    96172    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Sq, 1.e-6))) 
    97     pylab.title("I(q) for q=%g"%q) 
     173    pylab.axis('tight') 
     174    pylab.title("BCC S(q) for q=%g, dnn=%g d_factor=%g" % (q, dnn, d_factor)) 
    98175    pylab.xlabel("theta (degrees)") 
    99176    pylab.ylabel("phi (degrees)") 
    100177    cbar = pylab.colorbar() 
    101     cbar.set_label('log10 I(q)') 
     178    cbar.set_label('log10 S(q)') 
    102179    pylab.show() 
    103180 
    104181if __name__ == "__main__": 
    105     print("gauss", 20, gauss_quad(n=20)) 
    106     print("gauss", 76, gauss_quad(n=76)) 
    107     print("gauss", 150, gauss_quad(n=150)) 
    108     print("dblquad", scipy_dblquad()) 
    109     gridded_integrals(n=2**8+1) 
    110     gridded_integrals(n=2**10+1) 
    111     gridded_integrals(n=2**13+1) 
     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) 
    112190    #print("romberg", scipy_romberg()) 
    113     #plot(n=300) 
     191    plot(n=400) 
Note: See TracChangeset for help on using the changeset viewer.