Changeset c462169 in sasmodels


Ignore:
Timestamp:
Apr 2, 2018 12:48:10 PM (4 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:
29afc50
Parents:
b3703f5
git-author:
Paul Kienzle <pkienzle@…> (04/02/18 12:15:58)
git-committer:
Paul Kienzle <pkienzle@…> (04/02/18 12:48:10)
Message:

show u-sub integration as used in models in explore/asymint.py

File:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/asymint.py

    ra1c32c2 rc462169  
    3030import numpy as np 
    3131import mpmath as mp 
    32 from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10 
     32from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10, arccos 
    3333from numpy.polynomial.legendre import leggauss 
    3434from scipy.integrate import dblquad, simps, romb, romberg 
     
    4141shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1] 
    4242Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2] 
     43 
     44DTYPE = 'd' 
    4345 
    4446class MPenv: 
     
    7274    sas_sinx_x = staticmethod(sp.sas_sinx_x) 
    7375    pi = np.pi 
    74     mpf = staticmethod(float) 
     76    #mpf = staticmethod(float) 
     77    mpf = staticmethod(lambda x: np.array(x, DTYPE)) 
    7578 
    7679SLD = 3 
     
    220223    #A, B, C = 4450, 14000, 47 
    221224    #A, B, C = 445, 140, 47  # integer for the sake of mpf 
    222     A, B, C = 6800, 114, 1380 
    223     DA, DB, DC = 2300, 21, 58 
     225    A, B, C = 114, 1380, 6800 
     226    DA, DB, DC = 21, 58, 2300 
    224227    SLDA, SLDB, SLDC = "5", "-0.3", "11.5" 
     228    ## default parameters from sasmodels 
     229    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 
     230    ## swap A-B-C to C-B-A 
     231    #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 
    225232    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
    226233    #SLD_SOLVENT,CONTRAST = 0, 4 
     
    233240        DA, DC = DC, DA 
    234241        SLDA, SLDC = SLDC, SLDA 
     242    #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    235243    NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    236244    NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
     
    348356    return n**2, Iq 
    349357 
     358def gauss_quad_usub(q, n=150, dtype=DTYPE): 
     359    """ 
     360    Compute the integral using gaussian quadrature for n = 20, 76 or 150. 
     361 
     362    Use *u = sin theta* substitution, and restrict integration over a single 
     363    quadrant for shapes that are mirror symmetric about AB, AC and BC planes. 
     364 
     365    Note that this doesn't work for fcc/bcc paracrystals, which instead step 
     366    over the entire 4 pi surface uniformly in theta-phi. 
     367    """ 
     368    z, w = leggauss(n) 
     369    cos_theta = 0.5 * (z + 1) 
     370    theta = arccos(cos_theta) 
     371    phi = pi/2*(0.5 * (z + 1)) 
     372    Atheta, Aphi = np.meshgrid(theta, phi) 
     373    Aw = w[None, :] * w[:, None] 
     374    q, Atheta, Aphi, Aw = [np.asarray(v, dtype=dtype) for v in (q, Atheta, Aphi, Aw)] 
     375    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi) 
     376    Iq = np.sum(Zq*Aw)*0.25 
     377    return n**2, Iq 
     378 
    350379def gridded_2d(q, n=300): 
    351380    """ 
     
    395424    print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 
    396425    print("gauss-2049", *gauss_quad_2d(Q, n=2049)) 
     426    print("gauss-20 usub", *gauss_quad_usub(Q, n=20)) 
     427    print("gauss-76 usub", *gauss_quad_usub(Q, n=76)) 
     428    print("gauss-150 usub", *gauss_quad_usub(Q, n=150)) 
    397429    #gridded_2d(Q, n=2**8+1) 
    398430    gridded_2d(Q, n=2**10+1) 
Note: See TracChangeset for help on using the changeset viewer.