Changeset 2ab331f in sasmodels for explore


Ignore:
Timestamp:
Jan 11, 2018 11:00:36 AM (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:
2a7e20e
Parents:
226473d
Message:

play with non-dilute cylinder models

File:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/realspace.py

    r226473d r2ab331f  
    393393        pass 
    394394    n = len(points) 
    395     index = np.random.choice(n, size=1000) if n > 1000 else slice(None, None) 
     395    #print("len points", n) 
     396    index = np.random.choice(n, size=500) if n > 500 else slice(None, None) 
    396397    ax.scatter(points[index, 0], points[index, 1], points[index, 2], c=rho[index]) 
    397398    #low, high = points.min(axis=0), points.max(axis=0) 
     
    476477    q = np.logspace(-3, 0, 200) 
    477478    r = shape.r_bins(q, r_step=0.01) 
    478     sampling_density = 15000 / shape.volume() 
     479    sampling_density = 50000 / shape.volume() 
    479480    rho, points = shape.sample(sampling_density) 
    480481    Pr = calc_Pr(r, rho-rho_solvent, points) 
     
    489490def check_shape_2d(shape, fn=None, view=(0, 0, 0)): 
    490491    rho_solvent = 0 
    491     qx = qy = np.linspace(-1, 1, 100) 
     492    nq, qmax = 100, 0.1 
     493    nq, qmax = 100, 0.03 
     494    qx = np.linspace(0.0, qmax, nq) 
     495    qy = np.linspace(0.0, qmax, nq) 
    492496    Qx, Qy = np.meshgrid(qx, qy) 
    493     sampling_density = 50000 / shape.volume() 
     497    sampling_density = 15000 / shape.volume() 
     498    #t0 = time.time() 
    494499    rho, points = shape.sample(sampling_density) 
     500    #print("sample time", time.time() - t0) 
     501    t0 = time.time() 
    495502    Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 
     503    print("calc time", time.time() - t0) 
     504    theory = fn(Qx, Qy) if fn is not None else None 
    496505    Iqxy += 0.001 * Iqxy.max() 
    497     theory = fn(Qx, Qy)+0.001 if fn is not None else None 
     506    if theory is not None: 
     507        theory += 0.001 * theory.max() 
    498508 
    499509    import pylab 
     510    #plot_points(rho, points); pylab.figure() 
    500511    plot_calc_2d(qx, qy, Iqxy, theory=theory) 
    501512    pylab.show() 
     
    508519def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)): 
    509520    shape = EllipticalCylinder(radius, radius, length, rho) 
     521    fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
     522    check_shape_2d(shape, fn, view=view) 
     523 
     524def check_cylinder_2d_lattice(radius=25, length=125, rho=2., 
     525                              view=(0, 0, 0)): 
     526    nx, dx = 1, 2*radius 
     527    ny, dy = 30, 2*radius 
     528    nz, dz = 30, length 
     529    dx, dy, dz = 2*dx, 2*dy, 2*dz 
     530    def center(*args): 
     531        sigma = 0.333 
     532        space = 2 
     533        return [(space*n+np.random.randn()*sigma)*x for n, x in args] 
     534    shapes = [EllipticalCylinder(radius, radius, length, rho, 
     535                                 #center=(ix*dx, iy*dy, iz*dz) 
     536                                 orientation=np.random.randn(3)*0, 
     537                                 center=center((ix, dx), (iy, dy), (iz, dz)) 
     538                                ) 
     539              for ix in range(nx) 
     540              for iy in range(ny) 
     541              for iz in range(nz)] 
     542    shape = Composite(shapes) 
    510543    fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
    511544    check_shape_2d(shape, fn, view=view) 
     
    531564    check_cylinder(radius=10, length=40) 
    532565    #check_cylinder_2d(radius=10, length=40, view=(90,30,0)) 
     566    #check_cylinder_2d_lattice(radius=10, length=50, view=(90,30,0)) 
    533567    #check_sphere() 
    534568    #check_csbox() 
Note: See TracChangeset for help on using the changeset viewer.