source: sasmodels/explore/jitter.py @ 06ee63c

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

set projection used in theory to that selected in jitter.py

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