# source:sasmodels/sasmodels/jitter.py@0d5a655

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 0d5a655 was 0d5a655, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

make jitter explorer available from sasview app

• Property mode set to 100755
File size: 28.5 KB
Line
1#!/usr/bin/env python
2"""
3Jitter Explorer
4===============
5
6Application to explore orientation angle and angular dispersity.
7"""
8from __future__ import division, print_function
9
10import sys, os
11sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__))))
12
13import argparse
14
15import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot
16import matplotlib.pyplot as plt
17from matplotlib.widgets import Slider, CheckButtons
18from matplotlib import cm
19import numpy as np
20from numpy import pi, cos, sin, sqrt, exp, degrees, radians
21
22def draw_beam(ax, view=(0, 0)):
23    """
24    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1)
25    """
26    #ax.plot([0,0],[0,0],[1,-1])
27    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8)
28
29    steps = 25
30    u = np.linspace(0, 2 * np.pi, steps)
31    v = np.linspace(-1, 1, steps)
32
33    r = 0.02
34    x = r*np.outer(np.cos(u), np.ones_like(v))
35    y = r*np.outer(np.sin(u), np.ones_like(v))
36    z = 1.3*np.outer(np.ones_like(u), v)
37
38    theta, phi = view
39    shape = x.shape
40    points = np.matrix([x.flatten(), y.flatten(), z.flatten()])
41    points = Rz(phi)*Ry(theta)*points
42    x, y, z = [v.reshape(shape) for v in points]
43
44    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5)
45
46def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0)):
47    """
48    Represent jitter as a set of shapes at different orientations.
49    """
50    # set max diagonal to 0.95
51    scale = 0.95/sqrt(sum(v**2 for v in size))
52    size = tuple(scale*v for v in size)
53    draw_shape = draw_parallelepiped
54    #draw_shape = draw_ellipsoid
55
56    #np.random.seed(10)
57    #cloud = np.random.randn(10,3)
58    cloud = [
59        [-1, -1, -1],
60        [-1, -1,  0],
61        [-1, -1,  1],
62        [-1,  0, -1],
63        [-1,  0,  0],
64        [-1,  0,  1],
65        [-1,  1, -1],
66        [-1,  1,  0],
67        [-1,  1,  1],
68        [ 0, -1, -1],
69        [ 0, -1,  0],
70        [ 0, -1,  1],
71        [ 0,  0, -1],
72        [ 0,  0,  0],
73        [ 0,  0,  1],
74        [ 0,  1, -1],
75        [ 0,  1,  0],
76        [ 0,  1,  1],
77        [ 1, -1, -1],
78        [ 1, -1,  0],
79        [ 1, -1,  1],
80        [ 1,  0, -1],
81        [ 1,  0,  0],
82        [ 1,  0,  1],
83        [ 1,  1, -1],
84        [ 1,  1,  0],
85        [ 1,  1,  1],
86    ]
87    dtheta, dphi, dpsi = jitter
88    if dtheta == 0:
89        cloud = [v for v in cloud if v[0] == 0]
90    if dphi == 0:
91        cloud = [v for v in cloud if v[1] == 0]
92    if dpsi == 0:
93        cloud = [v for v in cloud if v[2] == 0]
94    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8)
95    scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist]
96    for point in cloud:
97        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]]
98        draw_shape(ax, size, view, delta, alpha=0.8)
99    for v in 'xyz':
100        a, b, c = size
101        lim = np.sqrt(a**2+b**2+c**2)
102        getattr(ax, 'set_'+v+'lim')([-lim, lim])
103        getattr(ax, v+'axis').label.set_text(v)
104
105def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1):
106    """Draw an ellipsoid."""
107    a,b,c = size
108    u = np.linspace(0, 2 * np.pi, steps)
109    v = np.linspace(0, np.pi, steps)
110    x = a*np.outer(np.cos(u), np.sin(v))
111    y = b*np.outer(np.sin(u), np.sin(v))
112    z = c*np.outer(np.ones_like(u), np.cos(v))
113    x, y, z = transform_xyz(view, jitter, x, y, z)
114
115    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
116
117    draw_labels(ax, view, jitter, [
118         ('c+', [ 0, 0, c], [ 1, 0, 0]),
119         ('c-', [ 0, 0,-c], [ 0, 0,-1]),
120         ('a+', [ a, 0, 0], [ 0, 0, 1]),
121         ('a-', [-a, 0, 0], [ 0, 0,-1]),
122         ('b+', [ 0, b, 0], [-1, 0, 0]),
123         ('b-', [ 0,-b, 0], [-1, 0, 0]),
124    ])
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 "a+" face.
149    if 1:
150        x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1])
151        y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1])
152        z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1])
153        x, y, z = transform_xyz(view, jitter, abs(x)*1.05, y, z)
154        ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1,0.6,0.6], alpha=alpha)
155
156    draw_labels(ax, view, jitter, [
157         ('c+', [ 0, 0, c], [ 1, 0, 0]),
158         ('c-', [ 0, 0,-c], [ 0, 0,-1]),
159         ('a+', [ a, 0, 0], [ 0, 0, 1]),
160         ('a-', [-a, 0, 0], [ 0, 0,-1]),
161         ('b+', [ 0, b, 0], [-1, 0, 0]),
162         ('b-', [ 0,-b, 0], [-1, 0, 0]),
163    ])
164
166    """Draw a sphere"""
167    u = np.linspace(0, 2 * np.pi, steps)
168    v = np.linspace(0, np.pi, steps)
169
170    x = radius * np.outer(np.cos(u), np.sin(v))
171    y = radius * np.outer(np.sin(u), np.sin(v))
172    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))
173    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
174
175PROJECTIONS = [
176    # in order of PROJECTION number; do not change without updating the
177    # constants in kernel_iq.c
178    'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance',
179    'azimuthal_equal_area',
180]
181def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian',
182              projection='equirectangular'):
183    """
184    Draw the dispersion mesh showing the theta-phi orientations at which
185    the model will be evaluated.
186
187    jitter projections
188    <https://en.wikipedia.org/wiki/List_of_map_projections>
189
190    equirectangular (standard latitude-longitude mesh)
191        <https://en.wikipedia.org/wiki/Equirectangular_projection>
192        Allows free movement in phi (around the equator), but theta is
193        limited to +/- 90, and points are cos-weighted. Jitter in phi is
194        uniform in weight along a line of latitude.  With small theta and
195        phi ranging over +/- 180 this forms a wobbling disk.  With small
196        phi and theta ranging over +/- 90 this forms a wedge like a slice
197        of an orange.
198    azimuthal_equidistance (Postel)
199        <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection>
200        Preserves distance from center, and so is an excellent map for
201        representing a bivariate gaussian on the surface.  Theta and phi
202        operate identically, cutting wegdes from the antipode of the viewing
203        angle.  This unfortunately does not allow free movement in either
204        theta or phi since the orthogonal wobble decreases to 0 as the body
205        rotates through 180 degrees.
206    sinusoidal (Sanson-Flamsteed, Mercator equal-area)
207        <https://en.wikipedia.org/wiki/Sinusoidal_projection>
208        Preserves arc length with latitude, giving bad behaviour at
209        theta near +/- 90.  Theta and phi operate somewhat differently,
210        so a system with a-b-c dtheta-dphi-dpsi will not give the same
211        value as one with b-a-c dphi-dtheta-dpsi, as would be the case
212        for azimuthal equidistance.  Free movement using theta or phi
213        uniform over +/- 180 will work, but not as well as equirectangular
214        phi, with theta being slightly worse.  Computationally it is much
215        cheaper for wide theta-phi meshes since it excludes points which
216        lie outside the sinusoid near theta +/- 90 rather than packing
217        them close together as in equirectangle.  Note that the poles
218        will be slightly overweighted for theta > 90 with the circle
219        from theta at 90+dt winding backwards around the pole, overlapping
220        the circle from theta at 90-dt.
221    Guyou (hemisphere-in-a-square) **not weighted**
222        <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection>
223        With tiling, allows rotation in phi or theta through +/- 180, with
224        uniform spacing.  Both theta and phi allow free rotation, with wobble
225        in the orthogonal direction reasonably well behaved (though not as
226        good as equirectangular phi). The forward/reverse transformations
227        relies on elliptic integrals that are somewhat expensive, so the
228        behaviour has to be very good to justify the cost and complexity.
229        The weighting function for each point has not yet been computed.
230        Note: run the module *guyou.py* directly and it will show the forward
231        and reverse mappings.
232    azimuthal_equal_area  **incomplete**
233        <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection>
234        Preserves the relative density of the surface patches.  Not that
235        useful and not completely implemented
236    Gauss-Kreuger **not implemented**
237        <https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Ellipsoidal_transverse_Mercator>
238        Should allow free movement in theta, but phi is distorted.
239    """
240    theta, phi, psi = view
241    dtheta, dphi, dpsi = jitter
242
243    t = np.linspace(-1, 1, n)
244    weights = np.ones_like(t)
245    if dist == 'gaussian':
246        t *= 3
247        weights = exp(-0.5*t**2)
248    elif dist == 'rectangle':
249        # Note: uses sasmodels ridiculous definition of rectangle width
250        t *= sqrt(3)
251    elif dist == 'uniform':
252        pass
253    else:
254        raise ValueError("expected dist to be gaussian, rectangle or uniform")
255
256    if projection == 'equirectangular':  #define PROJECTION 1
257        def rotate(theta_i, phi_j):
258            return Rx(phi_j)*Ry(theta_i)
259        def weight(theta_i, phi_j, wi, wj):
261    elif projection == 'sinusoidal':  #define PROJECTION 2
262        def rotate(theta_i, phi_j):
263            latitude = theta_i
265            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0
266            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
267            return Rx(longitude)*Ry(latitude)
268        def weight(theta_i, phi_j, wi, wj):
269            latitude = theta_i
271            w = 1 if abs(phi_j) < abs(scale)*180 else 0
272            return w*wi*wj
273    elif projection == 'guyou':  #define PROJECTION 3  (eventually?)
274        def rotate(theta_i, phi_j):
275            from guyou import guyou_invert
276            #latitude, longitude = guyou_invert([theta_i], [phi_j])
277            longitude, latitude = guyou_invert([phi_j], [theta_i])
278            return Rx(longitude[0])*Ry(latitude[0])
279        def weight(theta_i, phi_j, wi, wj):
280            return wi*wj
281    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry
282        def rotate(theta_i, phi_j):
283            latitude = sqrt(theta_i**2 + phi_j**2)
284            longitude = degrees(np.arctan2(phi_j, theta_i))
285            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
286            return Rz(longitude)*Ry(latitude)
287        def weight(theta_i, phi_j, wi, wj):
288            # Weighting for each point comes from the integral:
289            #     \int\int I(q, lat, log) sin(lat) dlat dlog
290            # We are doing a conformal mapping from disk to sphere, so we need
291            # a change of variables g(theta, phi) -> (lat, long):
292            #     lat, long = sqrt(theta^2 + phi^2), arctan(phi/theta)
293            # giving:
294            #     dtheta dphi = det(J) dlat dlong
295            # where J is the jacobian from the partials of g. Using
296            #     R = sqrt(theta^2 + phi^2),
297            # then
298            #     J = [[x/R, Y/R], -y/R^2, x/R^2]]
299            # and
300            #     det(J) = 1/R
301            # with the final integral being:
302            #    \int\int I(q, theta, phi) sin(R)/R dtheta dphi
303            #
304            # This does approximately the right thing, decreasing the weight
305            # of each point as you go farther out on the disk, but it hasn't
306            # yet been checked against the 1D integral results. Prior
307            # to declaring this "good enough" and checking that integrals
308            # work in practice, we will examine alternative mappings.
309            #
310            # The issue is that the mapping does not support the case of free
311            # rotation about a single axis correctly, with a small deviation
312            # in the orthogonal axis independent of the first axis.  Like the
313            # usual polar coordiates integration, the integrated sections
314            # form wedges, though at least in this case the wedge cuts through
315            # the entire sphere, and treats theta and phi identically.
316            latitude = sqrt(theta_i**2 + phi_j**2)
317            w = sin(radians(latitude))/latitude if latitude != 0 else 1
318            return w*wi*wj if latitude < 180 else 0
319    elif projection == 'azimuthal_equal_area':
320        def rotate(theta_i, phi_j):
321            R = min(1, sqrt(theta_i**2 + phi_j**2)/180)
322            latitude = 180-degrees(2*np.arccos(R))
323            longitude = degrees(np.arctan2(phi_j, theta_i))
324            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
325            return Rz(longitude)*Ry(latitude)
326        def weight(theta_i, phi_j, wi, wj):
327            latitude = sqrt(theta_i**2 + phi_j**2)
328            w = sin(radians(latitude))/latitude if latitude != 0 else 1
329            return w*wi*wj if latitude < 180 else 0
330    else:
331        raise ValueError("unknown projection %r"%projection)
332
333    # mesh in theta, phi formed by rotating z
334    z = np.matrix([[0], [0], [radius]])
335    points = np.hstack([rotate(theta_i, phi_j)*z
336                        for theta_i in dtheta*t
337                        for phi_j in dphi*t])
338    # select just the active points (i.e., those with phi < 180
339    w = np.array([weight(theta_i, phi_j, wi, wj)
340                  for wi, theta_i in zip(weights, dtheta*t)
341                  for wj, phi_j in zip(weights, dphi*t)])
342    #print(max(w), min(w), min(w[w>0]))
343    points = points[:, w>0]
344    w = w[w>0]
345    w /= max(w)
346
347    if 0: # Kent distribution
348        points = np.hstack([Rx(phi_j)*Ry(theta_i)*z for theta_i in 30*t for phi_j in 60*t])
349        xp, yp, zp = [np.array(v).flatten() for v in points]
352        w = exp(kappa*zp) #+ beta*(xp**2 + yp**2)
353        print(kappa, dtheta, radians(dtheta), min(w), max(w), sum(w))
355        #w /= sum(w)
356
357    # rotate relative to beam
358    points = orient_relative_to_beam(view, points)
359
360    x, y, z = [np.array(v).flatten() for v in points]
361    #plt.figure(2); plt.clf(); plt.hist(z, bins=np.linspace(-1, 1, 51))
362    ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.)
363
364def draw_labels(ax, view, jitter, text):
365    """
366    Draw text at a particular location.
367    """
368    labels, locations, orientations = zip(*text)
369    px, py, pz = zip(*locations)
370    dx, dy, dz = zip(*orientations)
371
372    px, py, pz = transform_xyz(view, jitter, px, py, pz)
373    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz)
374
375    # TODO: zdir for labels is broken, and labels aren't appearing.
376    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)):
377        zdir = np.asarray(zdir).flatten()
378        ax.text(p[0], p[1], p[2], label, zdir=zdir)
379
380# Definition of rotation matrices comes from wikipedia:
381#    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
382def Rx(angle):
383    """Construct a matrix to rotate points about *x* by *angle* degrees."""
385    R = [[1, 0, 0],
386         [0, +cos(a), -sin(a)],
387         [0, +sin(a), +cos(a)]]
388    return np.matrix(R)
389
390def Ry(angle):
391    """Construct a matrix to rotate points about *y* by *angle* degrees."""
393    R = [[+cos(a), 0, +sin(a)],
394         [0, 1, 0],
395         [-sin(a), 0, +cos(a)]]
396    return np.matrix(R)
397
398def Rz(angle):
399    """Construct a matrix to rotate points about *z* by *angle* degrees."""
401    R = [[+cos(a), -sin(a), 0],
402         [+sin(a), +cos(a), 0],
403         [0, 0, 1]]
404    return np.matrix(R)
405
406def transform_xyz(view, jitter, x, y, z):
407    """
408    Send a set of (x,y,z) points through the jitter and view transforms.
409    """
410    x, y, z = [np.asarray(v) for v in (x, y, z)]
411    shape = x.shape
412    points = np.matrix([x.flatten(),y.flatten(),z.flatten()])
413    points = apply_jitter(jitter, points)
414    points = orient_relative_to_beam(view, points)
415    x, y, z = [np.array(v).reshape(shape) for v in points]
416    return x, y, z
417
418def apply_jitter(jitter, points):
419    """
420    Apply the jitter transform to a set of points.
421
422    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
423    """
424    dtheta, dphi, dpsi = jitter
425    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
426    return points
427
428def orient_relative_to_beam(view, points):
429    """
430    Apply the view transform to a set of points.
431
432    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
433    """
434    theta, phi, psi = view
435    points = Rz(phi)*Ry(theta)*Rz(psi)*points
436    return points
437
438# translate between number of dimension of dispersity and the number of
439# points along each dimension.
440PD_N_TABLE = {
441    (0, 0, 0): (0, 0, 0),     # 0
442    (1, 0, 0): (100, 0, 0),   # 100
443    (0, 1, 0): (0, 100, 0),
444    (0, 0, 1): (0, 0, 100),
445    (1, 1, 0): (30, 30, 0),   # 900
446    (1, 0, 1): (30, 0, 30),
447    (0, 1, 1): (0, 30, 30),
448    (1, 1, 1): (15, 15, 15),  # 3375
449}
450
451def clipped_range(data, portion=1.0, mode='central'):
452    """
453    Determine range from data.
454
455    If *portion* is 1, use full range, otherwise use the center of the range
456    or the top of the range, depending on whether *mode* is 'central' or 'top'.
457    """
458    if portion == 1.0:
459        return data.min(), data.max()
460    elif mode == 'central':
461        data = np.sort(data.flatten())
462        offset = int(portion*len(data)/2 + 0.5)
463        return data[offset], data[-offset]
464    elif mode == 'top':
465        data = np.sort(data.flatten())
466        offset = int(portion*len(data) + 0.5)
467        return data[offset], data[-1]
468
469def draw_scattering(calculator, ax, view, jitter, dist='gaussian'):
470    """
471    Plot the scattering for the particular view.
472
473    *calculator* is returned from :func:build_model.  *ax* are the 3D axes
474    on which the data will be plotted.  *view* and *jitter* are the current
475    orientation and orientation dispersity.  *dist* is one of the sasmodels
476    weight distributions.
477    """
478    if dist == 'uniform':  # uniform is not yet in this branch
479        dist, scale = 'rectangle', 1/sqrt(3)
480    else:
481        scale = 1
482
483    # add the orientation parameters to the model parameters
484    theta, phi, psi = view
485    theta_pd, phi_pd, psi_pd = [scale*v for v in jitter]
486    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd>0, phi_pd>0, psi_pd>0)]
487    ## increase pd_n for testing jitter integration rather than simple viz
488    #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)]
489
490    pars = dict(
491        theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n,
492        phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n,
493        psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n,
494    )
495    pars.update(calculator.pars)
496
497    # compute the pattern
498    qx, qy = calculator._data.x_bins, calculator._data.y_bins
499    Iqxy = calculator(**pars).reshape(len(qx), len(qy))
500
501    # scale it and draw it
502    Iqxy = np.log(Iqxy)
503    if calculator.limits:
504        # use limits from orientation (0,0,0)
505        vmin, vmax = calculator.limits
506    else:
507        vmax = Iqxy.max()
508        vmin = vmax*10**-7
509        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top')
510    #print("range",(vmin,vmax))
511    #qx, qy = np.meshgrid(qx, qy)
512    if 0:
513        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i')
514        level[level<0] = 0
515        colors = plt.get_cmap()(level)
516        ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors)
517    elif 1:
518        ax.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1,
519                    levels=np.linspace(vmin, vmax, 24))
520    else:
521        ax.pcolormesh(qx, qy, Iqxy)
522
523def build_model(model_name, n=150, qmax=0.5, **pars):
524    """
525    Build a calculator for the given shape.
526
527    *model_name* is any sasmodels model.  *n* and *qmax* define an n x n mesh
528    on which to evaluate the model.  The remaining parameters are stored in
529    the returned calculator as *calculator.pars*.  They are used by
530    :func:draw_scattering to set the non-orientation parameters in the
531    calculation.
532
533    Returns a *calculator* function which takes a dictionary or parameters and
534    produces Iqxy.  The Iqxy value needs to be reshaped to an n x n matrix
535    for plotting.  See the :class:sasmodels.direct_model.DirectModel class
536    for details.
537    """
538    from sasmodels.core import load_model_info, build_model
539    from sasmodels.data import empty_data2D
540    from sasmodels.direct_model import DirectModel
541
543    model = build_model(model_info) #, dtype='double!')
544    q = np.linspace(-qmax, qmax, n)
545    data = empty_data2D(q, q)
546    calculator = DirectModel(data, model)
547
548    # stuff the values for non-orientation parameters into the calculator
549    calculator.pars = pars.copy()
550    calculator.pars.setdefault('backgound', 1e-3)
551
552    # fix the data limits so that we can see if the pattern fades
553    # under rotation or angular dispersion
554    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars)
555    Iqxy = np.log(Iqxy)
556    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top')
557    calculator.limits = vmin, vmax+1
558
559    return calculator
560
561def select_calculator(model_name, n=150, size=(10,40,100)):
562    """
563    Create a model calculator for the given shape.
564
565    *model_name* is one of sphere, cylinder, ellipsoid, triaxial_ellipsoid,
566    parallelepiped or bcc_paracrystal. *n* is the number of points to use
567    in the q range.  *qmax* is chosen based on model parameters for the
568    given model to show something intersting.
569
570    Returns *calculator* and tuple *size* (a,b,c) giving minor and major
571    equitorial axes and polar axis respectively.  See :func:build_model
572    for details on the returned calculator.
573    """
574    a, b, c = size
575    if model_name == 'sphere':
576        calculator = build_model('sphere', n=n, radius=c)
577        a = b = c
578    elif model_name == 'bcc_paracrystal':
579        calculator = build_model('bcc_paracrystal', n=n, dnn=c,
581        a = b = c
582    elif model_name == 'cylinder':
583        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c)
584        a = b
585    elif model_name == 'ellipsoid':
586        calculator = build_model('ellipsoid', n=n, qmax=1.0,
588        a = b
589    elif model_name == 'triaxial_ellipsoid':
590        calculator = build_model('triaxial_ellipsoid', n=n, qmax=0.5,
594    elif model_name == 'parallelepiped':
595        calculator = build_model('parallelepiped', n=n, a=a, b=b, c=c)
596    else:
597        raise ValueError("unknown model %s"%model_name)
598
599    return calculator, (a, b, c)
600
601SHAPES = [
602    'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal',
603    'cylinder', 'ellipsoid',
604    'sphere',
605 ]
606
607DISTRIBUTIONS = [
608    'gaussian', 'rectangle', 'uniform',
609]
610DIST_LIMITS = {
611    'gaussian': 30,
612    'rectangle': 90/sqrt(3),
613    'uniform': 90,
614}
615
616def run(model_name='parallelepiped', size=(10, 40, 100),
617        dist='gaussian', mesh=30,
618        projection='equirectangular'):
619    """
620    Show an interactive orientation and jitter demo.
621
622    *model_name* is one of: sphere, cylinder, ellipsoid, parallelepiped,
623    triaxial_ellipsoid, or bcc_paracrystal.
624
625    *size* gives the dimensions (a, b, c) of the shape.
626
627    *dist* is the type of dispersition: gaussian, rectangle, or uniform.
628
629    *mesh* is the number of points in the dispersion mesh.
630
631    *projection* is the map projection to use for the mesh: equirectangular,
632    sinusoidal, guyou, azimuthal_equidistance, or azimuthal_equal_area.
633    """
634    # projection number according to 1-order position in list, but
635    # only 1 and 2 are implemented so far.
636    from sasmodels import generate
637    generate.PROJECTION = PROJECTIONS.index(projection) + 1
638    if generate.PROJECTION > 2:
639        print("*** PROJECTION %s not implemented in scattering function ***"%projection)
640        generate.PROJECTION = 2
641
642    # set up calculator
643    calculator, size = select_calculator(model_name, n=150, size=size)
644
645    ## uncomment to set an independent the colour range for every view
646    ## If left commented, the colour range is fixed for all views
647    calculator.limits = None
648
649    ## initial view
650    #theta, dtheta = 70., 10.
651    #phi, dphi = -45., 3.
652    #psi, dpsi = -45., 3.
653    theta, phi, psi = 0, 0, 0
654    dtheta, dphi, dpsi = 0, 0, 0
655
656    ## create the plot window
657    #plt.hold(True)
658    plt.subplots(num=None, figsize=(5.5, 5.5))
659    plt.set_cmap('gist_earth')
660    plt.clf()
661    plt.gcf().canvas.set_window_title(projection)
662    #gs = gridspec.GridSpec(2,1,height_ratios=[4,1])
663    #ax = plt.subplot(gs[0], projection='3d')
664    ax = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d')
665    try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection
666        ax.axis('square')
667    except Exception:
668        pass
669
670    axcolor = 'lightgoldenrodyellow'
671
672    ## add control widgets to plot
673    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor)
674    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor)
675    axpsi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor)
676    stheta = Slider(axtheta, 'Theta', -90, 90, valinit=theta)
677    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi)
678    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi)
679
680    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor)
681    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor)
682    axdpsi= plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor)
683    # Note: using ridiculous definition of rectangle distribution, whose width
684    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep
685    # the maximum width to 90.
686    dlimit = DIST_LIMITS[dist]
687    sdtheta = Slider(axdtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta)
688    sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi)
689    sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi)
690
691    ## callback to draw the new view
692    def update(val, axis=None):
693        view = stheta.val, sphi.val, spsi.val
694        jitter = sdtheta.val, sdphi.val, sdpsi.val
695        # set small jitter as 0 if multiple pd dims
696        dims = sum(v > 0 for v in jitter)
697        limit = [0, 0, 0.5, 5][dims]
698        jitter = [0 if v < limit else v for v in jitter]
699        ax.cla()
700        draw_beam(ax, (0, 0))
701        draw_jitter(ax, view, jitter, dist=dist, size=size)
702        #draw_jitter(ax, view, (0,0,0))
703        draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection)
704        draw_scattering(calculator, ax, view, jitter, dist=dist)
705        plt.gcf().canvas.draw()
706
707    ## bind control widgets to view updater
708    stheta.on_changed(lambda v: update(v,'theta'))
709    sphi.on_changed(lambda v: update(v, 'phi'))
710    spsi.on_changed(lambda v: update(v, 'psi'))
711    sdtheta.on_changed(lambda v: update(v, 'dtheta'))
712    sdphi.on_changed(lambda v: update(v, 'dphi'))
713    sdpsi.on_changed(lambda v: update(v, 'dpsi'))
714
715    ## initialize view
716    update(None, 'phi')
717
718    ## go interactive
719    plt.show()
720
721def main():
722    parser = argparse.ArgumentParser(
723        description="Display jitter",
724        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
725        )
726    parser.add_argument('-p', '--projection', choices=PROJECTIONS, default=PROJECTIONS[0], help='coordinate projection')
727    parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths')
728    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='jitter distribution')
729    parser.add_argument('-m', '--mesh', type=int, default=30, help='#points in theta-phi mesh')
730    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape')
731    opts = parser.parse_args()
732    size = tuple(int(v) for v in opts.size.split(','))
733    run(opts.shape, size=size,
734        mesh=opts.mesh, dist=opts.distribution,
735        projection=opts.projection)
736
737if __name__ == "__main__":
738    main()
Note: See TracBrowser for help on using the repository browser.