Changeset 58a1be9 in sasmodels


Ignore:
Timestamp:
Feb 23, 2018 9:48:03 AM (6 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:
199bd07
Parents:
3db96b0
Message:

show paracrystal models as points in jitter explorer

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/jitter.py

    r42158d2 r58a1be9  
    4848    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    4949 
    50 def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0)): 
     50def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
     51    """Draw an ellipsoid.""" 
     52    a,b,c = size 
     53    u = np.linspace(0, 2 * np.pi, steps) 
     54    v = np.linspace(0, np.pi, steps) 
     55    x = a*np.outer(np.cos(u), np.sin(v)) 
     56    y = b*np.outer(np.sin(u), np.sin(v)) 
     57    z = c*np.outer(np.ones_like(u), np.cos(v)) 
     58    x, y, z = transform_xyz(view, jitter, x, y, z) 
     59 
     60    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
     61 
     62    draw_labels(ax, view, jitter, [ 
     63         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     64         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     65         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     66         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     67         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     68         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     69    ]) 
     70 
     71def draw_sc(ax, size, view, jitter, steps=None, alpha=1): 
     72    atoms = _build_sc() 
     73    _draw_crystal(ax, size, view, jitter, atoms=atoms) 
     74 
     75def draw_fcc(ax, size, view, jitter, steps=None, alpha=1): 
     76    # Build the simple cubic crystal 
     77    atoms = _build_sc() 
     78    # Define the centers for each face 
     79    # x planes at -1, 0, 1 have four centers per plane, at +/- 0.5 in y and z 
     80    x, y, z = ( 
     81        [-1]*4 + [0]*4 + [1]*4, 
     82        ([-0.5]*2 + [0.5]*2)*3, 
     83        [-0.5, 0.5]*12, 
     84    ) 
     85    # y and z planes can be generated by substituting x for y and z respectively 
     86    atoms.extend(zip(x+y+z, y+z+x, z+x+y)) 
     87    _draw_crystal(ax, size, view, jitter, atoms=atoms) 
     88 
     89def draw_bcc(ax, size, view, jitter, steps=None, alpha=1): 
     90    # Build the simple cubic crystal 
     91    atoms = _build_sc() 
     92    # Define the centers for each octant 
     93    # x plane at +/- 0.5 have four centers per plane at +/- 0.5 in y and z 
     94    x, y, z = ( 
     95        [-0.5]*4 + [0.5]*4, 
     96        ([-0.5]*2 + [0.5]*2)*2, 
     97        [-0.5, 0.5]*8, 
     98    ) 
     99    atoms.extend(zip(x, y, z)) 
     100    _draw_crystal(ax, size, view, jitter, atoms=atoms) 
     101 
     102def _draw_crystal(ax, size, view, jitter, steps=None, alpha=1, atoms=None): 
     103    atoms, size = np.asarray(atoms, 'd').T, np.asarray(size, 'd') 
     104    x, y, z = atoms*size[:, None] 
     105    x, y, z = transform_xyz(view, jitter, x, y, z) 
     106    ax.scatter([x[0]], [y[0]], [z[0]], c='yellow', marker='o') 
     107    ax.scatter(x[1:], y[1:], z[1:], c='r', marker='o') 
     108 
     109def _build_sc(): 
     110    # three planes of 9 dots for x at -1, 0 and 1 
     111    x, y, z = ( 
     112        [-1]*9 + [0]*9 + [1]*9, 
     113        ([-1]*3 + [0]*3 + [1]*3)*3, 
     114        [-1, 0, 1]*9, 
     115    ) 
     116    atoms = list(zip(x, y, z)) 
     117    #print(list(enumerate(atoms))) 
     118    # Pull the dot at (1, 0, 0) to the front of the list 
     119    # It will be highlighted in the view 
     120    index = 22 
     121    highlight = atoms[index] 
     122    del atoms[index] 
     123    atoms.insert(0, highlight) 
     124    return atoms 
     125 
     126def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
     127    """Draw a parallelepiped.""" 
     128    a, b, c = size 
     129    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
     130    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
     131    z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1]) 
     132    tri = np.array([ 
     133        # counter clockwise triangles 
     134        # z: up/down, x: right/left, y: front/back 
     135        [0,1,2], [3,2,1], # top face 
     136        [6,5,4], [5,6,7], # bottom face 
     137        [0,2,6], [6,4,0], # right face 
     138        [1,5,7], [7,3,1], # left face 
     139        [2,3,6], [7,6,3], # front face 
     140        [4,1,0], [5,1,4], # back face 
     141    ]) 
     142 
     143    x, y, z = transform_xyz(view, jitter, x, y, z) 
     144    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
     145 
     146    # Draw pink face on box. 
     147    # Since I can't control face color, instead draw a thin box situated just 
     148    # in front of the "c+" face.  Use the c face so that rotations about psi 
     149    # rotate that face. 
     150    if 1: 
     151        x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
     152        y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
     153        z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1]) 
     154        x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001) 
     155        ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1,0.6,0.6], alpha=alpha) 
     156 
     157    draw_labels(ax, view, jitter, [ 
     158         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     159         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     160         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     161         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     162         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     163         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     164    ]) 
     165 
     166def draw_sphere(ax, radius=10., steps=100): 
     167    """Draw a sphere""" 
     168    u = np.linspace(0, 2 * np.pi, steps) 
     169    v = np.linspace(0, np.pi, steps) 
     170 
     171    x = radius * np.outer(np.cos(u), np.sin(v)) 
     172    y = radius * np.outer(np.sin(u), np.sin(v)) 
     173    z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
     174    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
     175 
     176def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), 
     177                draw_shape=draw_parallelepiped): 
    51178    """ 
    52179    Represent jitter as a set of shapes at different orientations. 
     
    55182    scale = 0.95/sqrt(sum(v**2 for v in size)) 
    56183    size = tuple(scale*v for v in size) 
    57     draw_shape = draw_parallelepiped 
    58     #draw_shape = draw_ellipsoid 
    59184 
    60185    #np.random.seed(10) 
     
    106231        getattr(ax, 'set_'+v+'lim')([-lim, lim]) 
    107232        getattr(ax, v+'axis').label.set_text(v) 
    108  
    109 def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
    110     """Draw an ellipsoid.""" 
    111     a,b,c = size 
    112     u = np.linspace(0, 2 * np.pi, steps) 
    113     v = np.linspace(0, np.pi, steps) 
    114     x = a*np.outer(np.cos(u), np.sin(v)) 
    115     y = b*np.outer(np.sin(u), np.sin(v)) 
    116     z = c*np.outer(np.ones_like(u), np.cos(v)) 
    117     x, y, z = transform_xyz(view, jitter, x, y, z) 
    118  
    119     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    120  
    121     draw_labels(ax, view, jitter, [ 
    122          ('c+', [ 0, 0, c], [ 1, 0, 0]), 
    123          ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
    124          ('a+', [ a, 0, 0], [ 0, 0, 1]), 
    125          ('a-', [-a, 0, 0], [ 0, 0,-1]), 
    126          ('b+', [ 0, b, 0], [-1, 0, 0]), 
    127          ('b-', [ 0,-b, 0], [-1, 0, 0]), 
    128     ]) 
    129  
    130 def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
    131     """Draw a parallelepiped.""" 
    132     a, b, c = size 
    133     x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    134     y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    135     z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1]) 
    136     tri = np.array([ 
    137         # counter clockwise triangles 
    138         # z: up/down, x: right/left, y: front/back 
    139         [0,1,2], [3,2,1], # top face 
    140         [6,5,4], [5,6,7], # bottom face 
    141         [0,2,6], [6,4,0], # right face 
    142         [1,5,7], [7,3,1], # left face 
    143         [2,3,6], [7,6,3], # front face 
    144         [4,1,0], [5,1,4], # back face 
    145     ]) 
    146  
    147     x, y, z = transform_xyz(view, jitter, x, y, z) 
    148     ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    149  
    150     # Draw pink face on box. 
    151     # Since I can't control face color, instead draw a thin box situated just 
    152     # in front of the "a+" face. 
    153     if 1: 
    154         x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    155         y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    156         z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1]) 
    157         x, y, z = transform_xyz(view, jitter, abs(x)*1.05, y, z) 
    158         ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1,0.6,0.6], alpha=alpha) 
    159  
    160     draw_labels(ax, view, jitter, [ 
    161          ('c+', [ 0, 0, c], [ 1, 0, 0]), 
    162          ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
    163          ('a+', [ a, 0, 0], [ 0, 0, 1]), 
    164          ('a-', [-a, 0, 0], [ 0, 0,-1]), 
    165          ('b+', [ 0, b, 0], [-1, 0, 0]), 
    166          ('b-', [ 0,-b, 0], [-1, 0, 0]), 
    167     ]) 
    168  
    169 def draw_sphere(ax, radius=10., steps=100): 
    170     """Draw a sphere""" 
    171     u = np.linspace(0, 2 * np.pi, steps) 
    172     v = np.linspace(0, np.pi, steps) 
    173  
    174     x = radius * np.outer(np.cos(u), np.sin(v)) 
    175     y = radius * np.outer(np.sin(u), np.sin(v)) 
    176     z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
    177     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    178233 
    179234PROJECTIONS = [ 
     
    577632    """ 
    578633    a, b, c = size 
     634    d_factor = 0.06  # for paracrystal models 
    579635    if model_name == 'sphere': 
    580636        calculator = build_model('sphere', n=n, radius=c) 
    581637        a = b = c 
     638    elif model_name == 'sc_paracrystal': 
     639        a = b = c 
     640        dnn = c 
     641        radius = 0.5*c 
     642        calculator = build_model('sc_paracrystal', n=n, dnn=dnn, 
     643                                  d_factor=d_factor, radius=(1-d_factor)*radius, 
     644                                  background=0) 
     645    elif model_name == 'fcc_paracrystal': 
     646        a = b = c 
     647        # nearest neigbour distance dnn should be 2 radius, but I think the 
     648        # model uses lattice spacing rather than dnn in its calculations 
     649        dnn = 0.5*c 
     650        radius = sqrt(2)/4 * c 
     651        calculator = build_model('fcc_paracrystal', n=n, dnn=dnn, 
     652                                  d_factor=d_factor, radius=(1-d_factor)*radius, 
     653                                  background=0) 
    582654    elif model_name == 'bcc_paracrystal': 
    583         calculator = build_model('bcc_paracrystal', n=n, dnn=c, 
    584                                   d_factor=0.06, radius=40) 
    585655        a = b = c 
     656        # nearest neigbour distance dnn should be 2 radius, but I think the 
     657        # model uses lattice spacing rather than dnn in its calculations 
     658        dnn = 0.5*c 
     659        radius = sqrt(3)/2 * c 
     660        calculator = build_model('bcc_paracrystal', n=n, dnn=dnn, 
     661                                  d_factor=d_factor, radius=(1-d_factor)*radius, 
     662                                  background=0) 
    586663    elif model_name == 'cylinder': 
    587664        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c) 
     
    604681 
    605682SHAPES = [ 
    606     'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal', 
    607     'cylinder', 'ellipsoid', 
    608     'sphere', 
     683    'parallelepiped', 
     684    'sphere', 'ellipsoid', 'triaxial_ellipsoid', 
     685    'cylinder', 
     686    'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal', 
    609687 ] 
     688 
     689DRAW_SHAPES = { 
     690    'fcc_paracrystal': draw_fcc, 
     691    'bcc_paracrystal': draw_bcc, 
     692    'sc_paracrystal': draw_sc, 
     693    'parallelepiped': draw_parallelepiped, 
     694} 
    610695 
    611696DISTRIBUTIONS = [ 
     
    624709    Show an interactive orientation and jitter demo. 
    625710 
    626     *model_name* is one of: sphere, cylinder, ellipsoid, parallelepiped, 
    627     triaxial_ellipsoid, or bcc_paracrystal. 
     711    *model_name* is one of: sphere, ellipsoid, triaxial_ellipsoid, 
     712    parallelepiped, cylinder, or sc/fcc/bcc_paracrystal 
    628713 
    629714    *size* gives the dimensions (a, b, c) of the shape. 
     
    646731    # set up calculator 
    647732    calculator, size = select_calculator(model_name, n=150, size=size) 
     733    draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped) 
    648734 
    649735    ## uncomment to set an independent the colour range for every view 
     
    693779    sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
    694780 
     781 
    695782    ## callback to draw the new view 
    696783    def update(val, axis=None): 
     
    703790        ax.cla() 
    704791        draw_beam(ax, (0, 0)) 
    705         draw_jitter(ax, view, jitter, dist=dist, size=size) 
     792        draw_jitter(ax, view, jitter, dist=dist, size=size, draw_shape=draw_shape) 
    706793        #draw_jitter(ax, view, (0,0,0)) 
    707794        draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection) 
Note: See TracChangeset for help on using the changeset viewer.