Changeset 20fe0cd in sasmodels


Ignore:
Timestamp:
Oct 23, 2017 11:17:40 PM (5 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:
6e604f8
Parents:
36b3154
Message:

move paracrystal integration tests with the rest of the non-rotationally symmetric tests

Files:
1 added
1 deleted
4 edited

Legend:

Unmodified
Added
Removed
  • explore/angles.py

    • Property mode changed from 100644 to 100755
    rea60e08 r20fe0cd  
    229229 
    230230if 0: 
    231     comment("==== asymmetric (old) ====") 
     231    comment("==== asymmetric (3.x) ====") 
    232232    view_inv = Rz(-psi)*Rx(theta)*Ry(-(pi/2 - phi)) 
    233233    vprint(q_abc, view_inv*q_xy, "reverse view") 
     234    print("""  existing code 
     235    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; 
     236    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; 
     237    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; 
     238    """) 
  • explore/asymint.py

    • Property mode changed from 100644 to 100755
    r31eea1f r20fe0cd  
     1#!/usr/bin/env python 
    12""" 
    23Asymmetric shape integration 
     
    67 
    78import os, sys 
    8 sys.path.insert(0, os.path.dirname(os.path.dirname(__file__))) 
     9sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
    910 
    1011import numpy as np 
     
    1516import pylab 
    1617 
    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  
    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 
     18import sasmodels.special as sp 
     19 
     20# Need to parse shape early since it determines the kernel function 
     21# that will be used for the various integrators 
     22shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1] 
     23Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2] 
     24 
     25class MPenv: 
     26    sqrt = mp.sqrt 
     27    exp = mp.exp 
     28    expm1 = mp.expm1 
     29    cos = mp.cos 
     30    sin = mp.sin 
     31    tan = mp.tan 
     32    def sas_3j1x_x(x): 
     33        return 3*(mp.sin(x)/x - mp.cos(x))/(x*x) 
     34    def sas_2J1x_x(x): 
     35        return 2*mp.j1(x)/x 
     36    def sas_sinx_x(x): 
     37        return mp.sin(x)/x 
     38    pi = mp.pi 
     39    mpf = mp.mpf 
     40 
     41class NPenv: 
     42    sqrt = np.sqrt 
     43    exp = np.exp 
     44    expm1 = np.expm1 
     45    cos = np.cos 
     46    sin = np.sin 
     47    tan = np.tan 
     48    sas_3j1x_x = sp.sas_3j1x_x 
     49    sas_2J1x_x = sp.sas_2J1x_x 
     50    sas_sinx_x = sp.sas_sinx_x 
     51    pi = np.pi 
     52    mpf = float 
    2953 
    3054SLD = 3 
     
    3256CONTRAST = SLD - SLD_SOLVENT 
    3357 
    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) 
     58# Carefully code models so that mpmath will use full precision.  That means: 
     59#    * wrap inputs in env.mpf 
     60#    * don't use floating point constants, only integers 
     61#    * for division, make sure the numerator or denominator is env.mpf 
     62#    * use env.pi, env.sas_sinx_x, etc. for functions 
     63def make_parallelepiped(a, b, c, env=NPenv): 
     64    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c) 
     65    def Fq(qa, qb, qc): 
     66        siA = env.sas_sinx_x(0.5*a*qa/2) 
     67        siB = env.sas_sinx_x(0.5*b*qb/2) 
     68        siC = env.sas_sinx_x(0.5*c*qc/2) 
    3969        return siA * siB * siC 
     70    Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c) 
    4071    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  
    72 def make_cylinder(radius, length): 
    73     def Fq(qa, qb, qc): 
    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 
    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 
    88  
    89 def make_sphere(radius): 
    90     def Fq(qa, qb, qc): 
    91         q = sqrt(qa**2 + qb**2 + qc**2) 
    92         return sas_3j1x_x(q*radius) 
    93     volume = 4*pi*radius**3/3 
    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': 
     72    norm = CONTRAST**2*volume/10000 
     73    return norm, Fq 
     74 
     75def make_triaxial_ellipsoid(a, b, c, env=NPenv): 
     76    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c) 
     77    def Fq(qa, qb, qc): 
     78        qr = env.sqrt((a*qa)**2 + (b*qb)**2 + (c*qc)**2) 
     79        return env.sas_3j1x_x(qr) 
     80    Fq.__doc__ = "triaxial ellipse minor=%g, major=%g polar=%g"%(a, b, c) 
     81    volume = 4*env.pi*a*b*c/3 
     82    norm = CONTRAST**2*volume/10000 
     83    return norm, Fq 
     84 
     85def make_cylinder(radius, length, env=NPenv): 
     86    radius, length = env.mpf(radius), env.mpf(length) 
     87    def Fq(qa, qb, qc): 
     88        qab = env.sqrt(qa**2 + qb**2) 
     89        return env.sas_2J1x_x(qab*radius) * env.sas_sinx_x((qc*length)/2) 
     90    Fq.__doc__ = "cylinder radius=%g, length=%g"%(radius, length) 
     91    volume = env.pi*radius**2*length 
     92    norm = CONTRAST**2*volume/10000 
     93    return norm, Fq 
     94 
     95def make_sphere(radius, env=NPenv): 
     96    radius = env.mpf(radius) 
     97    def Fq(qa, qb, qc): 
     98        q = env.sqrt(qa**2 + qb**2 + qc**2) 
     99        return env.sas_3j1x_x(q*radius) 
     100    Fq.__doc__ = "sphere radius=%g"%(radius, ) 
     101    volume = 4*pi*radius**3 
     102    norm = CONTRAST**2*volume/10000 
     103    return norm, Fq 
     104 
     105def make_paracrystal(radius, dnn, d_factor, lattice='bcc', env=NPenv): 
     106    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor) 
     107    def sc(qa, qb, qc): 
     108        return qa, qb, qc 
     109    def bcc(qa, qb, qc): 
     110        a1 = (+qa + qb + qc)/2 
     111        a2 = (-qa - qb + qc)/2 
     112        a3 = (-qa + qb - qc)/2 
     113        return a1, a2, a3 
     114    def fcc(qa, qb, qc): 
     115        a1 = ( 0. + qb + qc)/2 
     116        a2 = (-qa + 0. + qc)/2 
     117        a3 = (-qa + qb + 0.)/2 
     118        return a1, a2, a3 
     119    lattice_fn = {'sc': sc, 'bcc': bcc, 'fcc': fcc}[lattice] 
     120    radius, dnn, d_factor = env.mpf(radius), env.mpf(dnn), env.mpf(d_factor) 
     121    def Fq(qa, qb, qc): 
     122        a1, a2, a3 = lattice_fn(qa, qb, qc) 
     123        # Note: paper says that different directions can have different 
     124        # distoration factors.  Easy enough to add to the code. 
     125        # This would definitely break 8-fold symmetry. 
     126        arg = -(dnn*d_factor)**2*(a1**2 + a2**2 + a3**2)/2 
     127        exp_arg = env.exp(arg) 
     128        den = [((exp_arg - 2*env.cos(dnn*a))*exp_arg + 1) for a in (a1, a2, a3)] 
     129        Sq = -env.expm1(2*arg)**3/(den[0]*den[1]*den[2]) 
     130 
     131        q = env.sqrt(qa**2 + qb**2 + qc**2) 
     132        Fq = env.sas_3j1x_x(q*radius) 
     133        # the kernel computes F(q)**2, but we need S(q)*F(q)**2 
     134        return env.sqrt(Sq)*Fq 
     135    Fq.__doc__ = "%s paracrystal a=%g da=%g r=%g"%(lattice, dnn, d_factor, radius) 
     136    def sphere_volume(r): return 4*env.pi*r**3/3 
     137    Vf = { 
     138        'sc': sphere_volume(radius/dnn), 
     139        'bcc': 2*sphere_volume(env.sqrt(3)/2*radius/dnn), 
     140        'fcc': 4*sphere_volume(radius/dnn/env.sqrt(2)), 
     141    }[lattice] 
     142    volume = sphere_volume(radius) 
     143    norm = CONTRAST**2*volume*Vf/10000 
     144    return norm, Fq 
     145 
     146if shape == 'sphere': 
     147    RADIUS = 50  # integer for the sake of mpf 
     148    NORM, KERNEL = make_sphere(radius=RADIUS) 
     149    NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv) 
     150elif shape == 'cylinder': 
    111151    #RADIUS, LENGTH = 10, 100000 
    112152    RADIUS, LENGTH = 10, 300  # integer for the sake of mpf 
    113153    NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) 
    114     NORM_MP, KERNEL_MP = make_cylinder_mp(radius=RADIUS, length=LENGTH) 
    115 elif shape == 'triellip': 
     154    NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv) 
     155elif shape == 'triaxial_ellipsoid': 
    116156    #A, B, C = 4450, 14000, 47 
    117157    A, B, C = 445, 140, 47  # integer for the sake of mpf 
    118158    NORM, KERNEL = make_triellip(A, B, C) 
    119     NORM_MP, KERNEL_MP = make_triellip_mp(A, B, C) 
     159    NORM_MP, KERNEL_MP = make_triellip(A, B, C, env=MPenv) 
    120160elif shape == 'parallelepiped': 
    121161    #A, B, C = 4450, 14000, 47 
    122162    A, B, C = 445, 140, 47  # integer for the sake of mpf 
    123163    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) 
     164    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
     165elif shape == 'paracrystal': 
     166    LATTICE = 'bcc' 
     167    #LATTICE = 'fcc' 
     168    #LATTICE = 'sc' 
     169    DNN, D_FACTOR = 220, '0.06'  # mpmath needs to initialize floats from string 
     170    RADIUS = 40  # integer for the sake of mpf 
     171    NORM, KERNEL = make_paracrystal( 
     172        radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE) 
     173    NORM_MP, KERNEL_MP = make_paracrystal( 
     174        radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv) 
    129175else: 
    130176    raise ValueError("Unknown shape %r"%shape) 
    131177 
     178# Note: hardcoded in mp_quad 
    132179THETA_LOW, THETA_HIGH = 0, pi 
    133180PHI_LOW, PHI_HIGH = 0, 2*pi 
     
    145192""" 
    146193 
    147  
    148 def mp_quad(q, shape): 
     194# 2D integration functions 
     195def mp_quad_2d(q, shape): 
    149196    evals = [0] 
    150197    def integrand(theta, phi): 
     
    160207    return evals[0], Iq 
    161208 
    162 def kernel(q, theta, phi): 
     209def kernel_2d(q, theta, phi): 
    163210    """ 
    164211    S(q) kernel for paracrystal forms. 
     
    170217    return NORM*KERNEL(qa, qb, qc)**2 
    171218 
    172 def scipy_dblquad(q): 
     219def scipy_dblquad_2d(q): 
    173220    """ 
    174221    Compute the integral using scipy dblquad.  This gets the correct answer 
     
    178225    def integrand(phi, theta): 
    179226        evals[0] += 1 
    180         Zq = kernel(q, theta=theta, phi=phi) 
     227        Zq = kernel_2d(q, theta=theta, phi=phi) 
    181228        return Zq*sin(theta) 
    182229    ans = dblquad(integrand, THETA_LOW, THETA_HIGH, lambda x: PHI_LOW, lambda x: PHI_HIGH)[0] 
     
    192239    def inner(phi, theta): 
    193240        evals[0] += 1 
    194         return kernel(q, theta=theta, phi=phi) 
     241        return kernel_2d(q, theta=theta, phi=phi) 
    195242    def outer(theta): 
    196243        Zq = romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(theta,)) 
     
    200247 
    201248 
    202 def semi_romberg(q, n=100): 
     249def semi_romberg_2d(q, n=100): 
    203250    """ 
    204251    Use 1D romberg integration in phi and regular simpsons rule in theta. 
     
    207254    def inner(phi, theta): 
    208255        evals[0] += 1 
    209         return kernel(q, theta=theta, phi=phi) 
     256        return kernel_2d(q, theta=theta, phi=phi) 
    210257    theta = np.linspace(THETA_LOW, THETA_HIGH, n) 
    211258    Zq = [romberg(inner, PHI_LOW, PHI_HIGH, divmax=100, args=(t,)) for t in theta] 
     
    213260    return evals[0], ans*SCALE/(4*pi) 
    214261 
    215 def gauss_quad(q, n=150): 
     262def gauss_quad_2d(q, n=150): 
    216263    """ 
    217264    Compute the integral using gaussian quadrature for n = 20, 76 or 150. 
    218265    """ 
    219     if n == 20: 
    220         z, w = Gauss20Z, Gauss20Wt 
    221     elif n == 76: 
    222         z, w = Gauss76Z, Gauss76Wt 
    223     elif n == 150: 
    224         z, w = Gauss150Z, Gauss150Wt 
    225     else: 
    226         z, w = leggauss(n) 
     266    z, w = leggauss(n) 
    227267    theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW 
    228268    phi = (PHI_HIGH-PHI_LOW)*(z + 1)/2 + PHI_LOW 
     
    230270    Aw = w[None, :] * w[:, None] 
    231271    sin_theta = abs(sin(Atheta)) 
    232     Zq = kernel(q=q, theta=Atheta, phi=Aphi) 
     272    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi) 
    233273    # change from [-1,1] x [-1,1] range to [0, pi] x [0, 2 pi] range 
    234274    dxdy_stretch = (THETA_HIGH-THETA_LOW)/2 * (PHI_HIGH-PHI_LOW)/2 
     
    236276    return n**2, Iq 
    237277 
    238  
    239 def gridded_integrals(q, n=300): 
     278def gridded_2d(q, n=300): 
    240279    """ 
    241280    Compute the integral on a regular grid using rectangular, trapezoidal, 
     
    246285    phi = np.linspace(PHI_LOW, PHI_HIGH, n) 
    247286    Atheta, Aphi = np.meshgrid(theta, phi) 
    248     Zq = kernel(q=q, theta=Atheta, phi=Aphi) 
     287    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi) 
    249288    Zq *= abs(sin(Atheta)) 
    250289    dx, dy = theta[1]-theta[0], phi[1]-phi[0] 
     
    254293    print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 
    255294 
    256 def plot(q, n=300): 
     295def plot_2d(q, n=300): 
    257296    """ 
    258297    Plot the 2D surface that needs to be integrated in order to compute 
     
    263302    phi = np.linspace(PHI_LOW, PHI_HIGH, n) 
    264303    Atheta, Aphi = np.meshgrid(theta, phi) 
    265     Zq = kernel(q=q, theta=Atheta, phi=Aphi) 
     304    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi) 
    266305    #Zq *= abs(sin(Atheta)) 
    267306    pylab.pcolor(degrees(theta), degrees(phi), log10(np.fmax(Zq, 1.e-6))) 
    268307    pylab.axis('tight') 
    269     pylab.title("%s Z(q) for q=%g" % (KERNEL.__name__, q)) 
     308    pylab.title("%s I(q,t) sin(t) for q=%g" % (KERNEL.__doc__, q)) 
    270309    pylab.xlabel("theta (degrees)") 
    271310    pylab.ylabel("phi (degrees)") 
     
    274313    pylab.show() 
    275314 
    276 if __name__ == "__main__": 
    277     Qstr = '0.005' 
    278     #Qstr = '0.8' 
    279     #Qstr = '0.0003' 
     315def main(Qstr): 
    280316    Q = float(Qstr) 
    281317    if shape == 'sphere': 
    282318        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)) 
    290     #gridded_integrals(Q, n=2**8+1) 
    291     gridded_integrals(Q, n=2**10+1) 
    292     #gridded_integrals(Q, n=2**13+1) 
    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) 
     319    print("gauss-20", *gauss_quad_2d(Q, n=20)) 
     320    print("gauss-76", *gauss_quad_2d(Q, n=76)) 
     321    print("gauss-150", *gauss_quad_2d(Q, n=150)) 
     322    print("gauss-500", *gauss_quad_2d(Q, n=500)) 
     323    #gridded_2d(Q, n=2**8+1) 
     324    gridded_2d(Q, n=2**10+1) 
     325    #gridded_2d(Q, n=2**13+1) 
     326    #gridded_2d(Q, n=2**15+1) 
     327    if shape != 'paracrystal':  # adaptive forms are too slow! 
     328        print("dblquad", *scipy_dblquad_2d(Q)) 
     329        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 
     330        print("romberg", *scipy_romberg_2d(Q)) 
     331        with mp.workprec(100): 
     332            print("mpmath", *mp_quad_2d(mp.mpf(Qstr), shape)) 
     333    plot_2d(Q, n=200) 
     334 
     335if __name__ == "__main__": 
     336    main(Qstr) 
  • explore/symint.py

    re47b06b r20fe0cd  
    1010import numpy as np 
    1111from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10 
     12from numpy.polynomial.legendre import leggauss 
    1213from scipy.integrate import dblquad, simps, romb, romberg 
    13 from scipy.special.orthogonal import p_roots 
    1414import pylab 
    1515 
     
    2121 
    2222SLD = 3.0 
    23 SLD_SOLVENT = 6.3 
     23SLD_SOLVENT = 6 
    2424CONTRAST = SLD - SLD_SOLVENT 
    2525 
     
    2929    cylinder.__doc__ = "cylinder radius=%g, length=%g"%(radius, length) 
    3030    volume = pi*radius**2*length 
    31     norm = 1e-4*volume*CONTRAST**2 
     31    norm = CONTRAST**2*volume/10000 
    3232    return norm, cylinder 
    3333 
     
    3737    long_cylinder.__doc__ = "long cylinder radius=%g, length=%g"%(radius, length) 
    3838    volume = pi*radius**2*length 
    39     norm = 1e-4*volume*CONTRAST**2*pi/length 
     39    norm = CONTRAST**2*volume/10000*pi/length 
    4040    return long_cylinder 
    4141 
     
    4646    sphere.__doc__ = "sphere radius=%g"%(radius,) 
    4747    volume = 4*pi*radius**3/3 
    48     norm = 1e-4*volume*CONTRAST**2 
     48    norm = CONTRAST**2*volume/10000 
    4949    return norm, sphere 
    5050 
     
    5353 
    5454 
    55 def kernel(q, theta): 
     55def kernel_1d(q, theta): 
    5656    """ 
    5757    S(q) kernel for paracrystal forms. 
     
    6161    return NORM*KERNEL(qab, qc)**2 
    6262 
    63  
    64 def gauss_quad(q, n=150): 
     63def gauss_quad_1d(q, n=150): 
    6564    """ 
    6665    Compute the integral using gaussian quadrature for n = 20, 76 or 150. 
    6766    """ 
    68     z, w = p_roots(n) 
     67    z, w = leggauss(n) 
    6968    theta = (THETA_HIGH-THETA_LOW)*(z + 1)/2 + THETA_LOW 
    7069    sin_theta = abs(sin(theta)) 
    71     Zq = kernel(q=q, theta=theta) 
     70    Zq = kernel_1d(q=q, theta=theta) 
    7271    return np.sum(Zq*w*sin_theta)*(THETA_HIGH-THETA_LOW)/2 
    7372 
    74  
    75 def gridded_integrals(q, n=300): 
     73def gridded_1d(q, n=300): 
    7674    """ 
    7775    Compute the integral on a regular grid using rectangular, trapezoidal, 
     
    8078    """ 
    8179    theta = np.linspace(THETA_LOW, THETA_HIGH, n) 
    82     Zq = kernel(q=q, theta=theta) 
     80    Zq = kernel_1d(q=q, theta=theta) 
    8381    Zq *= abs(sin(theta)) 
    8482    dx = theta[1]-theta[0] 
    85     print("rect", n, np.sum(Zq)*dx*SCALE) 
    86     print("trapz", n, np.trapz(Zq, dx=dx)*SCALE) 
    87     print("simpson", n, simps(Zq, dx=dx)*SCALE) 
    88     print("romb", n, romb(Zq, dx=dx)*SCALE) 
     83    print("rect-%d"%n, np.sum(Zq)*dx*SCALE) 
     84    print("trapz-%d"%n, np.trapz(Zq, dx=dx)*SCALE) 
     85    print("simpson-%d"%n, simps(Zq, dx=dx)*SCALE) 
     86    print("romb-%d"%n, romb(Zq, dx=dx)*SCALE) 
    8987 
    90 def scipy_romberg(q): 
     88def scipy_romberg_1d(q): 
    9189    """ 
    9290    Compute the integral using romberg integration.  This function does not 
     
    9694    def outer(theta): 
    9795        evals[0] += 1 
    98         return kernel(q, theta=theta)*abs(sin(theta)) 
     96        return kernel_1d(q, theta=theta)*abs(sin(theta)) 
    9997    result = romberg(outer, THETA_LOW, THETA_HIGH, divmax=100)*SCALE 
    10098    print("scipy romberg", evals[0], result) 
    10199 
    102 def plot(q, n=300): 
     100def plot_1d(q, n=300): 
    103101    """ 
    104     Plot the 2D surface that needs to be integrated in order to compute 
    105     the BCC S(q) at a particular q, dnn and d_factor.  *n* is the number 
    106     of points in the grid. 
     102    Plot the function that needs to be integrated in order to compute 
     103    the I(q) at a particular q.  *n* is the number of points in the grid. 
    107104    """ 
    108105    theta = np.linspace(THETA_LOW, THETA_HIGH, n) 
    109     Zq = kernel(q=q, theta=theta) 
     106    Zq = kernel_1d(q=q, theta=theta) 
    110107    Zq *= abs(sin(theta)) 
    111108    pylab.semilogy(degrees(theta), np.fmax(Zq, 1.e-6), label="Q=%g"%q) 
     
    116113def Iq_trapz(q, n): 
    117114    theta = np.linspace(THETA_LOW, THETA_HIGH, n) 
    118     Zq = kernel(q=q, theta=theta) 
     115    Zq = kernel_1d(q=q, theta=theta) 
    119116    Zq *= abs(sin(theta)) 
    120117    dx = theta[1]-theta[0] 
     
    142139    Q = 0.386 
    143140    for n in (20, 76, 150, 300, 1000): #, 10000, 30000): 
    144         print("gauss", n, gauss_quad(Q, n=n)) 
     141        print("gauss-%d"%n, gauss_quad(Q, n=n)) 
    145142    for k in (8, 10, 13, 16, 19): 
    146         gridded_integrals(Q, n=2**k+1) 
     143        gridded_1d(Q, n=2**k+1) 
    147144    #print("inf cyl", 0, long_cyl(Q)) 
    148145    #scipy_romberg(Q) 
    149146 
    150     plot(0.386, n=2000) 
    151     plot(0.5, n=2000) 
    152     plot(0.8, n=2000) 
     147    plot_1d(0.386, n=2000) 
     148    plot_1d(0.5, n=2000) 
     149    plot_1d(0.8, n=2000) 
    153150    pylab.legend() 
    154151    pylab.figure() 
  • sasmodels/model_test.py

    r8db25bf r20fe0cd  
    202202                ] 
    203203            tests = smoke_tests 
     204            #tests = [] 
    204205            if self.info.tests is not None: 
    205206                tests += self.info.tests 
Note: See TracChangeset for help on using the changeset viewer.