source: sasmodels/explore/jitter.py @ 8cfb486

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

update orientation docs

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