Changeset 1dd7a45f in sasmodels for explore


Ignore:
Timestamp:
Jul 5, 2017 6:30:55 PM (7 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:
e47b06b
Parents:
fac6657
Message:

cylinder: add equation for very long cylinder to numerical integration explorer

File:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/symint.py

    rfac23b3 r1dd7a45f  
    3131    norm = 1e-4*volume*CONTRAST**2 
    3232    return norm, cylinder 
     33 
     34def make_inf_cylinder(radius, length): 
     35    def inf_cylinder(q): 
     36        return norm/q * sas_2J1x_x(q*radius)**2 
     37    inf_cylinder.__doc__ = "inf cylinder radius=%g, length=%g"%(radius, length) 
     38    volume = pi*radius**2*length 
     39    norm = 1e-4*volume*CONTRAST**2*pi/length 
     40    return inf_cylinder 
    3341 
    3442def make_sphere(radius): 
     
    115123def plot_Iq(q, n, form="trapz"): 
    116124    if form == "trapz": 
    117         I = np.array([Iq_trapz(qk, n) for qk in q]) 
     125        Iq = np.array([Iq_trapz(qk, n) for qk in q]) 
    118126    elif form == "gauss": 
    119         I = np.array([gauss_quad(qk, n) for qk in q]) 
    120     pylab.loglog(q, I, label="%s, n=%d"%(form, n)) 
     127        Iq = np.array([gauss_quad(qk, n) for qk in q]) 
     128    pylab.loglog(q, Iq, label="%s, n=%d"%(form, n)) 
    121129    pylab.xlabel("q (1/A)") 
    122130    pylab.ylabel("Iq (1/cm)") 
    123131    pylab.title(KERNEL.__doc__ + " I(q) circular average") 
    124     return q, I 
     132    return Iq 
    125133 
    126 NORM, KERNEL = make_cylinder(radius=10., length=100000.) 
    127 #NORM, KERNEL = make_cylinder(radius=10., length=50000.) 
    128 #NORM, KERNEL = make_cylinder(radius=10., length=20000.) 
    129 #NORM, KERNEL = make_cylinder(radius=10., length=10000.) 
    130 #NORM, KERNEL = make_cylinder(radius=10., length=5000.) 
    131 #NORM, KERNEL = make_cylinder(radius=10., length=1000.) 
    132 #NORM, KERNEL = make_cylinder(radius=10., length=500.) 
    133 #NORM, KERNEL = make_cylinder(radius=10., length=100.) 
    134 #NORM, KERNEL = make_cylinder(radius=10., length=30.) 
     134radius = 10. 
     135length = 1e5 
     136NORM, KERNEL = make_cylinder(radius=radius, length=length) 
     137inf_cyl = make_inf_cylinder(radius=radius, length=length) 
    135138#NORM, KERNEL = make_sphere(radius=50.) 
    136139 
     140 
    137141if __name__ == "__main__": 
    138     Q = 0.8 
     142    Q = 0.386 
    139143    for n in (20, 76, 150, 300, 1000): #, 10000, 30000): 
    140144        print("gauss", n, gauss_quad(Q, n=n)) 
    141145    for k in (8, 10, 13, 16, 19): 
    142146        gridded_integrals(Q, n=2**k+1) 
     147    #print("inf cyl", 0, inf_cyl(Q)) 
    143148    #scipy_romberg(Q) 
     149 
     150    plot(0.386, n=2000) 
    144151    plot(0.5, n=2000) 
    145     plot(0.6, n=2000) 
    146152    plot(0.8, n=2000) 
    147153    pylab.legend() 
    148154    pylab.figure() 
    149     #plot_Iq(np.logspace(-3, 0, 400), n=2**19+1, form="trapz") 
    150     q1, I1 = plot_Iq(np.logspace(-3, 0, 400), n=2**16+1, form="trapz") 
    151     #plot_Iq(np.logspace(-3, 0, 400), n=2**10+1, form="trapz") 
    152     q2, I2 = plot_Iq(np.logspace(-3, 0, 400), n=1024, form="gauss") 
    153     #plot_Iq(np.logspace(-3, 0, 400), n=300, form="gauss") 
    154     plot_Iq(np.logspace(-3, 0, 400), n=150, form="gauss") 
    155     plot_Iq(np.logspace(-3, 0, 400), n=76, form="gauss") 
     155 
     156    q = np.logspace(-3, 0, 400) 
     157    I1 = inf_cyl(q) 
     158    I2 = plot_Iq(q, n=2**19+1, form="trapz") 
     159    #plot_Iq(q, n=2**16+1, form="trapz") 
     160    #plot_Iq(q, n=2**10+1, form="trapz") 
     161    plot_Iq(q, n=1024, form="gauss") 
     162    #plot_Iq(q, n=300, form="gauss") 
     163    #plot_Iq(q, n=150, form="gauss") 
     164    #plot_Iq(q, n=76, form="gauss") 
     165    pylab.loglog(q, inf_cyl(q), label="limit") 
    156166    pylab.legend() 
    157     #pylab.figure() 
    158     #pylab.semilogx(q1, (I2 - I1)/I1) 
     167 
     168    pylab.figure() 
     169    pylab.semilogx(q, (I2 - I1)/I1) 
     170 
    159171    pylab.show() 
Note: See TracChangeset for help on using the changeset viewer.