Changeset 8678a34 in sasmodels


Ignore:
Timestamp:
Apr 9, 2017 5:45:20 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:
a0ebc96
Parents:
85190c2
Message:

code cleanup on jitter visualization

File:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/jitter.py

    r85190c2 r8678a34  
    33dispersity and possible replacement algorithms. 
    44""" 
     5import sys 
     6 
    57import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
    68import matplotlib.pyplot as plt 
    79from matplotlib.widgets import Slider, CheckButtons 
    810from matplotlib import cm 
    9  
    1011import numpy as np 
    1112from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1213 
    13 def draw_beam(ax): 
     14def draw_beam(ax, view=(0, 0)): 
    1415    #ax.plot([0,0],[0,0],[1,-1]) 
    1516    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
     
    2223    x = r*np.outer(np.cos(u), np.ones_like(v)) 
    2324    y = r*np.outer(np.sin(u), np.ones_like(v)) 
    24     z = np.outer(np.ones_like(u), v) 
     25    z = 1.3*np.outer(np.ones_like(u), v) 
     26 
     27    theta, phi = view 
     28    shape = x.shape 
     29    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
     30    points = Rz(phi)*Ry(theta)*points 
     31    x, y, z = [v.reshape(shape) for v in points] 
    2532 
    2633    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    2734 
    28 def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi): 
    29     size=[0.1, 0.4, 1.0] 
    30     view=[theta, phi, psi] 
    31     shimmy=[0,0,0] 
    32     #draw_shape = draw_parallelepiped 
    33     draw_shape = draw_ellipsoid 
     35def draw_jitter(ax, view, jitter): 
     36    size = [0.1, 0.4, 1.0] 
     37    draw_shape = draw_parallelepiped 
     38    #draw_shape = draw_ellipsoid 
    3439 
    3540    #np.random.seed(10) 
     
    6469        [ 1,  1,  1], 
    6570    ] 
     71    dtheta, dphi, dpsi = jitter 
    6672    if dtheta == 0: 
    6773        cloud = [v for v in cloud if v[0] == 0] 
     
    7076    if dpsi == 0: 
    7177        cloud = [v for v in cloud if v[2] == 0] 
    72     draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8) 
     78    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    7379    for point in cloud: 
    74         shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]] 
    75         draw_shape(ax, size, view, shimmy, alpha=0.8) 
     80        delta = [dtheta*point[0], dphi*point[1], dpsi*point[2]] 
     81        draw_shape(ax, size, view, delta, alpha=0.8) 
    7682    for v in 'xyz': 
    7783        a, b, c = size 
     
    8086        getattr(ax, v+'axis').label.set_text(v) 
    8187 
    82 def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1): 
     88def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
    8389    a,b,c = size 
    84     theta, phi, psi = view 
    85     dtheta, dphi, dpsi = shimmy 
    86  
    8790    u = np.linspace(0, 2 * np.pi, steps) 
    8891    v = np.linspace(0, np.pi, steps) 
     
    9093    y = b*np.outer(np.sin(u), np.sin(v)) 
    9194    z = c*np.outer(np.ones_like(u), np.cos(v)) 
    92  
    93     shape = x.shape 
    94     points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
    95     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    96     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    97     x,y,z = [v.reshape(shape) for v in points] 
     95    x, y, z = transform_xyz(view, jitter, x, y, z) 
    9896 
    9997    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    10098 
    101 def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 
     99    draw_labels(ax, view, jitter, [ 
     100         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     101         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     102         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     103         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     104         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     105         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     106    ]) 
     107 
     108def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
    102109    a,b,c = size 
    103     theta, phi, psi = view 
    104     dtheta, dphi, dpsi = shimmy 
    105  
    106110    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    107111    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
     
    118122    ]) 
    119123 
    120     points = np.matrix([x,y,z]) 
    121     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    122     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    123  
    124     x,y,z = [np.array(v).flatten() for v in points] 
     124    x, y, z = transform_xyz(view, jitter, x, y, z) 
    125125    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    126126 
    127 def draw_sphere(ax, radius=10., steps=100): 
    128     u = np.linspace(0, 2 * np.pi, steps) 
    129     v = np.linspace(0, np.pi, steps) 
    130  
    131     x = radius * np.outer(np.cos(u), np.sin(v)) 
    132     y = radius * np.outer(np.sin(u), np.sin(v)) 
    133     z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
    134     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    135  
    136 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    137     theta_center = radians(theta) 
    138     phi_center = radians(phi) 
    139     flow_center = radians(flow) 
    140     dtheta = radians(dtheta) 
    141     dphi = radians(dphi) 
    142  
    143     # 10 point 3-sigma gaussian weights 
    144     t = np.linspace(-3., 3., 11) 
     127    draw_labels(ax, view, jitter, [ 
     128         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     129         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     130         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     131         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     132         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     133         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     134    ]) 
     135 
     136def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gauss'): 
     137    theta, phi, psi = view 
     138    dtheta, dphi, dpsi = jitter 
    145139    if dist == 'gauss': 
     140        t = np.linspace(-3, 3, n) 
    146141        weights = exp(-0.5*t**2) 
    147142    elif dist == 'rect': 
     143        t = np.linspace(0, 1, n) 
    148144        weights = np.ones_like(t) 
    149145    else: 
    150146        raise ValueError("expected dist to be 'gauss' or 'rect'") 
    151     theta = dtheta*t 
    152     phi = dphi*t 
    153  
    154     x = radius * np.outer(cos(phi), cos(theta)) 
    155     y = radius * np.outer(sin(phi), cos(theta)) 
    156     z = radius * np.outer(np.ones_like(phi), sin(theta)) 
    157     #w = np.outer(weights, weights*abs(cos(dtheta*t))) 
    158     w = np.outer(weights, weights*abs(cos(theta))) 
    159  
    160     x, y, z, w = [v.flatten() for v in (x,y,z,w)] 
    161     x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center) 
    162  
    163     ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
    164  
    165 def rotate(x, y, z, phi, theta, psi): 
    166     R = Rz(phi)*Ry(theta)*Rz(psi) 
    167     p = np.vstack([x,y,z]) 
    168     return R*p 
     147 
     148    # mesh in theta, phi formed by rotating z 
     149    z = np.matrix([[0], [0], [radius]]) 
     150    points = np.hstack([Rx(phi_i)*Ry(theta_i)*z 
     151                        for theta_i in dtheta*t 
     152                        for phi_i in dphi*t]) 
     153    # rotate relative to beam 
     154    points = orient_relative_to_beam(view, points) 
     155 
     156    w = np.outer(weights, weights) 
     157 
     158    x, y, z = [np.array(v).flatten() for v in points] 
     159    ax.scatter(x, y, z, c=w.flatten(), marker='o', vmin=0., vmax=1.) 
    169160 
    170161def Rx(angle): 
     
    188179         [0., 0., 1.]] 
    189180    return np.matrix(R) 
     181 
     182def transform_xyz(view, jitter, x, y, z): 
     183    x, y, z = [np.asarray(v) for v in (x, y, z)] 
     184    shape = x.shape 
     185    points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     186    points = apply_jitter(jitter, points) 
     187    points = orient_relative_to_beam(view, points) 
     188    x, y, z = [np.array(v).reshape(shape) for v in points] 
     189    return x, y, z 
     190 
     191def apply_jitter(jitter, points): 
     192    dtheta, dphi, dpsi = jitter 
     193    points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
     194    return points 
     195 
     196def orient_relative_to_beam(view, points): 
     197    theta, phi, psi = view 
     198    points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     199    return points 
     200 
     201def draw_labels(ax, view, jitter, text): 
     202    labels, locations, orientations = zip(*text) 
     203    px, py, pz = zip(*locations) 
     204    dx, dy, dz = zip(*orientations) 
     205 
     206    px, py, pz = transform_xyz(view, jitter, px, py, pz) 
     207    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz) 
     208 
     209    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 
     210        zdir = np.asarray(zdir).flatten() 
     211        ax.text(p[0], p[1], p[2], label, zdir=zdir) 
     212 
     213def draw_sphere(ax, radius=10., steps=100): 
     214    u = np.linspace(0, 2 * np.pi, steps) 
     215    v = np.linspace(0, np.pi, steps) 
     216 
     217    x = radius * np.outer(np.cos(u), np.sin(v)) 
     218    y = radius * np.outer(np.sin(u), np.sin(v)) 
     219    z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
     220    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    190221 
    191222def main(): 
     
    206237 
    207238    axcolor = 'lightgoldenrodyellow' 
     239 
    208240    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    209241    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
     
    212244    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    213245    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
     246 
    214247    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    215248    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
     
    217250    sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta) 
    218251    sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi) 
    219     sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi) 
     252    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dpsi) 
    220253 
    221254    def update(val, axis=None): 
    222         theta, phi, psi = stheta.val, sphi.val, spsi.val 
    223         dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val 
     255        view = stheta.val, sphi.val, spsi.val 
     256        jitter = sdtheta.val, sdphi.val, sdpsi.val 
    224257        ax.cla() 
    225         draw_beam(ax) 
    226         draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi) 
    227         #if not axis.startswith('d'): 
    228         #    ax.view_init(elev=theta, azim=phi) 
     258        draw_beam(ax, (0, 0)) 
     259        if 0: 
     260            draw_jitter(ax, view, jitter) 
     261        else: 
     262            draw_jitter(ax, view, (0,0,0)) 
     263            draw_mesh(ax, view, jitter) 
    229264        plt.gcf().canvas.draw() 
    230265 
Note: See TracChangeset for help on using the changeset viewer.