Changeset 58a1be9 in sasmodels

Timestamp:

Feb 23, 2018 11:48:03 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:

199bd07

Parents:

3db96b0

Message:

show paracrystal models as points in jitter explorer

File:

• sasmodels/jitter.py (modified) (9 diffs)

sasmodels/jitter.py

@@ -48 +48 @@
     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5)
 
-def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0)):
+def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1):
+    """Draw an ellipsoid."""
+    a,b,c = size
+    u = np.linspace(0, 2 * np.pi, steps)
+    v = np.linspace(0, np.pi, steps)
+    x = a*np.outer(np.cos(u), np.sin(v))
+    y = b*np.outer(np.sin(u), np.sin(v))
+    z = c*np.outer(np.ones_like(u), np.cos(v))
+    x, y, z = transform_xyz(view, jitter, x, y, z)
+
+    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
+
+    draw_labels(ax, view, jitter, [
+         ('c+', [ 0, 0, c], [ 1, 0, 0]),
+         ('c-', [ 0, 0,-c], [ 0, 0,-1]),
+         ('a+', [ a, 0, 0], [ 0, 0, 1]),
+         ('a-', [-a, 0, 0], [ 0, 0,-1]),
+         ('b+', [ 0, b, 0], [-1, 0, 0]),
+         ('b-', [ 0,-b, 0], [-1, 0, 0]),
+    ])
+
+def draw_sc(ax, size, view, jitter, steps=None, alpha=1):
+    atoms = _build_sc()
+    _draw_crystal(ax, size, view, jitter, atoms=atoms)
+
+def draw_fcc(ax, size, view, jitter, steps=None, alpha=1):
+    # Build the simple cubic crystal
+    atoms = _build_sc()
+    # Define the centers for each face
+    # x planes at -1, 0, 1 have four centers per plane, at +/- 0.5 in y and z
+    x, y, z = (
+        [-1]*4 + [0]*4 + [1]*4,
+        ([-0.5]*2 + [0.5]*2)*3,
+        [-0.5, 0.5]*12,
+    )
+    # y and z planes can be generated by substituting x for y and z respectively
+    atoms.extend(zip(x+y+z, y+z+x, z+x+y))
+    _draw_crystal(ax, size, view, jitter, atoms=atoms)
+
+def draw_bcc(ax, size, view, jitter, steps=None, alpha=1):
+    # Build the simple cubic crystal
+    atoms = _build_sc()
+    # Define the centers for each octant
+    # x plane at +/- 0.5 have four centers per plane at +/- 0.5 in y and z
+    x, y, z = (
+        [-0.5]*4 + [0.5]*4,
+        ([-0.5]*2 + [0.5]*2)*2,
+        [-0.5, 0.5]*8,
+    )
+    atoms.extend(zip(x, y, z))
+    _draw_crystal(ax, size, view, jitter, atoms=atoms)
+
+def _draw_crystal(ax, size, view, jitter, steps=None, alpha=1, atoms=None):
+    atoms, size = np.asarray(atoms, 'd').T, np.asarray(size, 'd')
+    x, y, z = atoms*size[:, None]
+    x, y, z = transform_xyz(view, jitter, x, y, z)
+    ax.scatter([x[0]], [y[0]], [z[0]], c='yellow', marker='o')
+    ax.scatter(x[1:], y[1:], z[1:], c='r', marker='o')
+
+def _build_sc():
+    # three planes of 9 dots for x at -1, 0 and 1
+    x, y, z = (
+        [-1]*9 + [0]*9 + [1]*9,
+        ([-1]*3 + [0]*3 + [1]*3)*3,
+        [-1, 0, 1]*9,
+    )
+    atoms = list(zip(x, y, z))
+    #print(list(enumerate(atoms)))
+    # Pull the dot at (1, 0, 0) to the front of the list
+    # It will be highlighted in the view
+    index = 22
+    highlight = atoms[index]
+    del atoms[index]
+    atoms.insert(0, highlight)
+    return atoms
+
+def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1):
+    """Draw a parallelepiped."""
+    a, b, c = size
+    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
+    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
+    z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1])
+    tri = np.array([
+        # counter clockwise triangles
+        # z: up/down, x: right/left, y: front/back
+        [0,1,2], [3,2,1], # top face
+        [6,5,4], [5,6,7], # bottom face
+        [0,2,6], [6,4,0], # right face
+        [1,5,7], [7,3,1], # left face
+        [2,3,6], [7,6,3], # front face
+        [4,1,0], [5,1,4], # back face
+    ])
+
+    x, y, z = transform_xyz(view, jitter, x, y, z)
+    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
+
+    # Draw pink face on box.
+    # Since I can't control face color, instead draw a thin box situated just
+    # in front of the "c+" face.  Use the c face so that rotations about psi
+    # rotate that face.
+    if 1:
+        x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
+        y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
+        z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1])
+        x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001)
+        ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1,0.6,0.6], alpha=alpha)
+
+    draw_labels(ax, view, jitter, [
+         ('c+', [ 0, 0, c], [ 1, 0, 0]),
+         ('c-', [ 0, 0,-c], [ 0, 0,-1]),
+         ('a+', [ a, 0, 0], [ 0, 0, 1]),
+         ('a-', [-a, 0, 0], [ 0, 0,-1]),
+         ('b+', [ 0, b, 0], [-1, 0, 0]),
+         ('b-', [ 0,-b, 0], [-1, 0, 0]),
+    ])
+
+def draw_sphere(ax, radius=10., steps=100):
+    """Draw a sphere"""
+    u = np.linspace(0, 2 * np.pi, steps)
+    v = np.linspace(0, np.pi, steps)
+
+    x = radius * np.outer(np.cos(u), np.sin(v))
+    y = radius * np.outer(np.sin(u), np.sin(v))
+    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))
+    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
+
+def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0),
+                draw_shape=draw_parallelepiped):
     """
     Represent jitter as a set of shapes at different orientations.
@@ -55 +182 @@
     scale = 0.95/sqrt(sum(v**2 for v in size))
     size = tuple(scale*v for v in size)
-    draw_shape = draw_parallelepiped
-    #draw_shape = draw_ellipsoid
 
     #np.random.seed(10)
@@ -106 +231 @@
         getattr(ax, 'set_'+v+'lim')([-lim, lim])
         getattr(ax, v+'axis').label.set_text(v)
-
-def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1):
-    """Draw an ellipsoid."""
-    a,b,c = size
-    u = np.linspace(0, 2 * np.pi, steps)
-    v = np.linspace(0, np.pi, steps)
-    x = a*np.outer(np.cos(u), np.sin(v))
-    y = b*np.outer(np.sin(u), np.sin(v))
-    z = c*np.outer(np.ones_like(u), np.cos(v))
-    x, y, z = transform_xyz(view, jitter, x, y, z)
-
-    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
-
-    draw_labels(ax, view, jitter, [
-         ('c+', [ 0, 0, c], [ 1, 0, 0]),
-         ('c-', [ 0, 0,-c], [ 0, 0,-1]),
-         ('a+', [ a, 0, 0], [ 0, 0, 1]),
-         ('a-', [-a, 0, 0], [ 0, 0,-1]),
-         ('b+', [ 0, b, 0], [-1, 0, 0]),
-         ('b-', [ 0,-b, 0], [-1, 0, 0]),
-    ])
-
-def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1):
-    """Draw a parallelepiped."""
-    a, b, c = size
-    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
-    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
-    z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1])
-    tri = np.array([
-        # counter clockwise triangles
-        # z: up/down, x: right/left, y: front/back
-        [0,1,2], [3,2,1], # top face
-        [6,5,4], [5,6,7], # bottom face
-        [0,2,6], [6,4,0], # right face
-        [1,5,7], [7,3,1], # left face
-        [2,3,6], [7,6,3], # front face
-        [4,1,0], [5,1,4], # back face
-    ])
-
-    x, y, z = transform_xyz(view, jitter, x, y, z)
-    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
-
-    # Draw pink face on box.
-    # Since I can't control face color, instead draw a thin box situated just
-    # in front of the "a+" face.
-    if 1:
-        x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
-        y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
-        z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1])
-        x, y, z = transform_xyz(view, jitter, abs(x)*1.05, y, z)
-        ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1,0.6,0.6], alpha=alpha)
-
-    draw_labels(ax, view, jitter, [
-         ('c+', [ 0, 0, c], [ 1, 0, 0]),
-         ('c-', [ 0, 0,-c], [ 0, 0,-1]),
-         ('a+', [ a, 0, 0], [ 0, 0, 1]),
-         ('a-', [-a, 0, 0], [ 0, 0,-1]),
-         ('b+', [ 0, b, 0], [-1, 0, 0]),
-         ('b-', [ 0,-b, 0], [-1, 0, 0]),
-    ])
-
-def draw_sphere(ax, radius=10., steps=100):
-    """Draw a sphere"""
-    u = np.linspace(0, 2 * np.pi, steps)
-    v = np.linspace(0, np.pi, steps)
-
-    x = radius * np.outer(np.cos(u), np.sin(v))
-    y = radius * np.outer(np.sin(u), np.sin(v))
-    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))
-    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
 
 PROJECTIONS = [
@@ -577 +632 @@
     """
     a, b, c = size
+    d_factor = 0.06  # for paracrystal models
     if model_name == 'sphere':
         calculator = build_model('sphere', n=n, radius=c)
         a = b = c
+    elif model_name == 'sc_paracrystal':
+        a = b = c
+        dnn = c
+        radius = 0.5*c
+        calculator = build_model('sc_paracrystal', n=n, dnn=dnn,
+                                  d_factor=d_factor, radius=(1-d_factor)*radius,
+                                  background=0)
+    elif model_name == 'fcc_paracrystal':
+        a = b = c
+        # nearest neigbour distance dnn should be 2 radius, but I think the
+        # model uses lattice spacing rather than dnn in its calculations
+        dnn = 0.5*c
+        radius = sqrt(2)/4 * c
+        calculator = build_model('fcc_paracrystal', n=n, dnn=dnn,
+                                  d_factor=d_factor, radius=(1-d_factor)*radius,
+                                  background=0)
     elif model_name == 'bcc_paracrystal':
-        calculator = build_model('bcc_paracrystal', n=n, dnn=c,
-                                  d_factor=0.06, radius=40)
         a = b = c
+        # nearest neigbour distance dnn should be 2 radius, but I think the
+        # model uses lattice spacing rather than dnn in its calculations
+        dnn = 0.5*c
+        radius = sqrt(3)/2 * c
+        calculator = build_model('bcc_paracrystal', n=n, dnn=dnn,
+                                  d_factor=d_factor, radius=(1-d_factor)*radius,
+                                  background=0)
     elif model_name == 'cylinder':
         calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c)
@@ -604 +681 @@
 
 SHAPES = [
-    'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal',
-    'cylinder', 'ellipsoid',
-    'sphere',
+    'parallelepiped',
+    'sphere', 'ellipsoid', 'triaxial_ellipsoid',
+    'cylinder',
+    'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal',
  ]
+
+DRAW_SHAPES = {
+    'fcc_paracrystal': draw_fcc,
+    'bcc_paracrystal': draw_bcc,
+    'sc_paracrystal': draw_sc,
+    'parallelepiped': draw_parallelepiped,
+}
 
 DISTRIBUTIONS = [
@@ -624 +709 @@
     Show an interactive orientation and jitter demo.
 
-    *model_name* is one of: sphere, cylinder, ellipsoid, parallelepiped,
-    triaxial_ellipsoid, or bcc_paracrystal.
+    *model_name* is one of: sphere, ellipsoid, triaxial_ellipsoid,
+    parallelepiped, cylinder, or sc/fcc/bcc_paracrystal
 
     *size* gives the dimensions (a, b, c) of the shape.
@@ -646 +731 @@
     # set up calculator
     calculator, size = select_calculator(model_name, n=150, size=size)
+    draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped)
 
     ## uncomment to set an independent the colour range for every view
@@ -693 +779 @@
     sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi)
 
+
     ## callback to draw the new view
     def update(val, axis=None):
@@ -703 +790 @@
         ax.cla()
         draw_beam(ax, (0, 0))
-        draw_jitter(ax, view, jitter, dist=dist, size=size)
+        draw_jitter(ax, view, jitter, dist=dist, size=size, draw_shape=draw_shape)
         #draw_jitter(ax, view, (0,0,0))
         draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection)