Changes in / [8064d5e:4057e06] in sasmodels

Files:

• explore/beta/data_files/richard_test.txt (deleted)
• explore/beta/data_files/richard_test2.txt (deleted)
• explore/beta/data_files/richard_test3.txt (deleted)
• explore/beta/data_files/richard_test4.txt (deleted)
• explore/beta/data_files/richard_test5.txt (deleted)
• explore/beta/data_files/richard_test6.txt (deleted)
• explore/beta/data_files/sasfit_ellipsoid_schulz_IQBD.txt (added)
• explore/beta/data_files/sasfit_ellipsoid_schulz_IQD.txt (added)
• explore/beta/data_files/sasfit_ellipsoid_schulz_IQSD.txt (added)
• explore/beta/data_files/sasfit_sphere_schulz_IQBD.txt (added)
• explore/beta/data_files/sasfit_sphere_schulz_IQD.txt (added)
• explore/beta/data_files/sasfit_sphere_schulz_IQSD.txt (added)
• explore/beta/sasfit_compare.py (modified) (2 diffs)
• sasmodels/jitter.py (modified) (26 diffs)
• sasmodels/models/pearl_necklace.c (modified) (2 diffs)

explore/beta/sasfit_compare.py

@@ -505 +505 @@
     }
 
-    Q, IQ = load_sasfit(data_file('richard_test.txt'))
-    Q, IQSD = load_sasfit(data_file('richard_test2.txt'))
-    Q, IQBD = load_sasfit(data_file('richard_test3.txt'))
+    Q, IQ = load_sasfit(data_file('sasfit_sphere_schulz_IQD.txt'))
+    Q, IQSD = load_sasfit(data_file('sasfit_sphere_schulz_IQSD.txt'))
+    Q, IQBD = load_sasfit(data_file('sasfit_sphere_schulz_IQBD.txt'))
     target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
     actual = sphere_r(Q, norm="sasfit", **pars)
@@ -526 +526 @@
     }
 
-    Q, IQ = load_sasfit(data_file('richard_test4.txt'))
-    Q, IQSD = load_sasfit(data_file('richard_test5.txt'))
-    Q, IQBD = load_sasfit(data_file('richard_test6.txt'))
+    Q, IQ = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQD.txt'))
+    Q, IQSD = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQSD.txt'))
+    Q, IQBD = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQBD.txt'))
     target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD)
     actual = ellipsoid_pe(Q, norm="sasfit", **pars)

sasmodels/jitter.py

@@ -1 +1 @@
 #!/usr/bin/env python
+# -*- coding: utf-8 -*-
 """
 Jitter Explorer
@@ -5 +6 @@
 
 Application to explore orientation angle and angular dispersity.
+
+From the command line::
+
+    # Show docs
+    python -m sasmodels.jitter --help
+
+    # Guyou projection jitter, uniform over 20 degree theta and 10 in phi
+    python -m sasmodels.jitter --projection=guyou --dist=uniform --jitter=20,10,0
+
+From a jupyter cell::
+
+    import ipyvolume as ipv
+    from sasmodels import jitter
+    import importlib; importlib.reload(jitter)
+    jitter.set_plotter("ipv")
+
+    size = (10, 40, 100)
+    view = (20, 0, 0)
+
+    #size = (15, 15, 100)
+    #view = (60, 60, 0)
+
+    dview = (0, 0, 0)
+    #dview = (5, 5, 0)
+    #dview = (15, 180, 0)
+    #dview = (180, 15, 0)
+
+    projection = 'equirectangular'
+    #projection = 'azimuthal_equidistance'
+    #projection = 'guyou'
+    #projection = 'sinusoidal'
+    #projection = 'azimuthal_equal_area'
+
+    dist = 'uniform'
+    #dist = 'gaussian'
+
+    jitter.run(size=size, view=view, jitter=dview, dist=dist, projection=projection)
+    #filename = projection+('_theta' if dview[0] == 180 else '_phi' if dview[1] == 180 else '')
+    #ipv.savefig(filename+'.png')
 """
 from __future__ import division, print_function
@@ -10 +50 @@
 import argparse
 
-try: # CRUFT: travis-ci does not support mpl_toolkits.mplot3d
-    import mpl_toolkits.mplot3d  # Adds projection='3d' option to subplot
-except ImportError:
-    pass
-
-import matplotlib as mpl
-import matplotlib.pyplot as plt
-from matplotlib.widgets import Slider
 import numpy as np
 from numpy import pi, cos, sin, sqrt, exp, degrees, radians
 
-def draw_beam(axes, view=(0, 0)):
+def draw_beam(axes, view=(0, 0), alpha=0.5, steps=25):
     """
     Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1)
     """
     #axes.plot([0,0],[0,0],[1,-1])
-    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8)
-
-    steps = 25
+    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=alpha)
+
     u = np.linspace(0, 2 * np.pi, steps)
-    v = np.linspace(-1, 1, steps)
+    v = np.linspace(-1, 1, 2)
 
     r = 0.02
@@ -42 +73 @@
     points = Rz(phi)*Ry(theta)*points
     x, y, z = [v.reshape(shape) for v in points]
-
-    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5)
+    axes.plot_surface(x, y, z, color='yellow', alpha=alpha)
+
+    # TODO: draw endcaps on beam
+    ## Drawing tiny balls on the end will work
+    #draw_sphere(axes, radius=0.02, center=(0, 0, 1.3), color='yellow', alpha=alpha)
+    #draw_sphere(axes, radius=0.02, center=(0, 0, -1.3), color='yellow', alpha=alpha)
+    ## The following does not work
+    #triangles = [(0, i+1, i+2) for i in range(steps-2)]
+    #x_cap, y_cap = x[:, 0], y[:, 0]
+    #for z_cap in z[:, 0], z[:, -1]:
+    #    axes.plot_trisurf(x_cap, y_cap, z_cap, triangles,
+    #                      color='yellow', alpha=alpha)
+
 
 def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1):
@@ -55 +97 @@
     x, y, z = transform_xyz(view, jitter, x, y, z)
 
-    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
+    axes.plot_surface(x, y, z, color='w', alpha=alpha)
 
     draw_labels(axes, view, jitter, [
@@ -124 +166 @@
     return atoms
 
-def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1):
+def draw_box(axes, size, view):
+    a, b, c = size
+    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1])
+    y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1])
+    z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1])
+    x, y, z = transform_xyz(view, None, x, y, z)
+    def draw(i, j):
+        axes.plot([x[i],x[j]], [y[i], y[j]], [z[i], z[j]], color='black')
+    draw(0, 1)
+    draw(0, 2)
+    draw(0, 3)
+    draw(7, 4)
+    draw(7, 5)
+    draw(7, 6)
+
+def draw_parallelepiped(axes, size, view, jitter, steps=None,
+                        color=(0.6, 1.0, 0.6), alpha=1):
     """Draw a parallelepiped."""
     a, b, c = size
@@ -142 +200 @@
 
     x, y, z = transform_xyz(view, jitter, x, y, z)
-    axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
-
-    # Draw pink face on box.
+    axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha,
+                      linewidth=0)
+
+    # Colour the c+ face of the box.
     # Since I can't control face color, instead draw a thin box situated just
     # in front of the "c+" face.  Use the c face so that rotations about psi
     # rotate that face.
-    if 1:
+    if 0:
+        color = (1, 0.6, 0.6)  # pink
         x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1])
         y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1])
         z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1])
         x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001)
-        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha)
+        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha)
 
     draw_labels(axes, view, jitter, [
@@ -164 +224 @@
     ])
 
-def draw_sphere(axes, radius=10., steps=100):
+def draw_sphere(axes, radius=0.5, steps=25, center=(0,0,0), color='w', alpha=1.):
     """Draw a sphere"""
     u = np.linspace(0, 2 * np.pi, steps)
     v = np.linspace(0, np.pi, steps)
 
-    x = radius * np.outer(np.cos(u), np.sin(v))
-    y = radius * np.outer(np.sin(u), np.sin(v))
-    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))
-    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w')
-
-def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0),
-                draw_shape=draw_parallelepiped):
+    x = radius * np.outer(np.cos(u), np.sin(v)) + center[0]
+    y = radius * np.outer(np.sin(u), np.sin(v)) + center[1]
+    z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) + center[2]
+    axes.plot_surface(x, y, z, color=color, alpha=alpha)
+    #axes.plot_wireframe(x, y, z)
+
+def draw_axes(axes, origin=(-1, -1, -1), length=(2, 2, 2)):
+    x, y, z = origin
+    dx, dy, dz = length
+    axes.plot([x, x+dx], [y, y], [z, z], color='black')
+    axes.plot([x, x], [y, y+dy], [z, z], color='black')
+    axes.plot([x, x], [y, y], [z, z+dz], color='black')
+
+def draw_person_on_sphere(axes, view, height=0.5, radius=0.5):
+    limb_offset = height * 0.05
+    head_radius = height * 0.10
+    head_height = height - head_radius
+    neck_length = head_radius * 0.50
+    shoulder_height = height - 2*head_radius - neck_length
+    torso_length = shoulder_height * 0.55
+    torso_radius = torso_length * 0.30
+    leg_length = shoulder_height - torso_length
+    arm_length = torso_length * 0.90
+
+    def _draw_part(x, z):
+        y = np.zeros_like(x)
+        xp, yp, zp = transform_xyz(view, None, x, y, z + radius)
+        axes.plot(xp, yp, zp, color='k')
+
+    # circle for head
+    u = np.linspace(0, 2 * np.pi, 40)
+    x = head_radius * np.cos(u)
+    z = head_radius * np.sin(u) + head_height
+    _draw_part(x, z)
+
+    # rectangle for body
+    x = np.array([-torso_radius, torso_radius, torso_radius, -torso_radius, -torso_radius])
+    z = np.array([0., 0, torso_length, torso_length, 0]) + leg_length
+    _draw_part(x, z)
+
+    # arms
+    x = np.array([-torso_radius - limb_offset, -torso_radius - limb_offset, -torso_radius])
+    z = np.array([shoulder_height - arm_length, shoulder_height, shoulder_height])
+    _draw_part(x, z)
+    _draw_part(-x, z)
+
+    # legs
+    x = np.array([-torso_radius + limb_offset, -torso_radius + limb_offset])
+    z = np.array([0, leg_length])
+    _draw_part(x, z)
+    _draw_part(-x, z)
+
+    limits = [-radius-height, radius+height]
+    axes.set_xlim(limits)
+    axes.set_ylim(limits)
+    axes.set_zlim(limits)
+    axes.set_axis_off()
+
+def draw_jitter(axes, view, jitter, dist='gaussian',
+                size=(0.1, 0.4, 1.0),
+                draw_shape=draw_parallelepiped,
+                projection='equirectangular',
+                alpha=0.8,
+                views=None):
     """
     Represent jitter as a set of shapes at different orientations.
     """
+    project, project_weight = get_projection(projection)
+
     # set max diagonal to 0.95
     scale = 0.95/sqrt(sum(v**2 for v in size))
     size = tuple(scale*v for v in size)
 
-    #np.random.seed(10)
-    #cloud = np.random.randn(10,3)
-    cloud = [
-        [-1, -1, -1],
-        [-1, -1, +0],
-        [-1, -1, +1],
-        [-1, +0, -1],
-        [-1, +0, +0],
-        [-1, +0, +1],
-        [-1, +1, -1],
-        [-1, +1, +0],
-        [-1, +1, +1],
-        [+0, -1, -1],
-        [+0, -1, +0],
-        [+0, -1, +1],
-        [+0, +0, -1],
-        [+0, +0, +0],
-        [+0, +0, +1],
-        [+0, +1, -1],
-        [+0, +1, +0],
-        [+0, +1, +1],
-        [+1, -1, -1],
-        [+1, -1, +0],
-        [+1, -1, +1],
-        [+1, +0, -1],
-        [+1, +0, +0],
-        [+1, +0, +1],
-        [+1, +1, -1],
-        [+1, +1, +0],
-        [+1, +1, +1],
-    ]
     dtheta, dphi, dpsi = jitter
-    if dtheta == 0:
-        cloud = [v for v in cloud if v[0] == 0]
-    if dphi == 0:
-        cloud = [v for v in cloud if v[1] == 0]
-    if dpsi == 0:
-        cloud = [v for v in cloud if v[2] == 0]
-    draw_shape(axes, size, view, [0, 0, 0], steps=100, alpha=0.8)
-    scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist]
-    for point in cloud:
-        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]]
-        draw_shape(axes, size, view, delta, alpha=0.8)
+    base = {'gaussian':3, 'rectangle':sqrt(3), 'uniform':1}[dist]
+    def steps(delta):
+        if views is None:
+            n = max(3, min(25, 2*int(base*delta/5)))
+        else:
+            n = views
+        return base*delta*np.linspace(-1, 1, n) if delta > 0 else [0.]
+    for theta in steps(dtheta):
+        for phi in steps(dphi):
+            for psi in steps(dpsi):
+                w = project_weight(theta, phi, 1.0, 1.0)
+                if w > 0:
+                    dview = project(theta, phi, psi)
+                    draw_shape(axes, size, view, dview, alpha=alpha)
     for v in 'xyz':
         a, b, c = size
         lim = np.sqrt(a**2 + b**2 + c**2)
         getattr(axes, 'set_'+v+'lim')([-lim, lim])
-        getattr(axes, v+'axis').label.set_text(v)
+        #getattr(axes, v+'axis').label.set_text(v)
 
 PROJECTIONS = [
@@ -238 +329 @@
     'azimuthal_equal_area',
 ]
-def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian',
-              projection='equirectangular'):
-    """
-    Draw the dispersion mesh showing the theta-phi orientations at which
-    the model will be evaluated.
-
+def get_projection(projection):
+
+    """
     jitter projections
     <https://en.wikipedia.org/wiki/List_of_map_projections>
@@ -299 +387 @@
     # TODO: try Kent distribution instead of a gaussian warped by projection
 
-    dist_x = np.linspace(-1, 1, n)
-    weights = np.ones_like(dist_x)
-    if dist == 'gaussian':
-        dist_x *= 3
-        weights = exp(-0.5*dist_x**2)
-    elif dist == 'rectangle':
-        # Note: uses sasmodels ridiculous definition of rectangle width
-        dist_x *= sqrt(3)
-    elif dist == 'uniform':
-        pass
-    else:
-        raise ValueError("expected dist to be gaussian, rectangle or uniform")
-
     if projection == 'equirectangular':  #define PROJECTION 1
-        def _rotate(theta_i, phi_j):
-            return Rx(phi_j)*Ry(theta_i)
+        def _project(theta_i, phi_j, psi):
+            latitude, longitude = theta_i, phi_j
+            return latitude, longitude, psi
+            #return Rx(phi_j)*Ry(theta_i)
         def _weight(theta_i, phi_j, w_i, w_j):
             return w_i*w_j*abs(cos(radians(theta_i)))
     elif projection == 'sinusoidal':  #define PROJECTION 2
-        def _rotate(theta_i, phi_j):
+        def _project(theta_i, phi_j, psi):
             latitude = theta_i
             scale = cos(radians(latitude))
             longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0
             #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
-            return Rx(longitude)*Ry(latitude)
-        def _weight(theta_i, phi_j, w_i, w_j):
+            return latitude, longitude, psi
+            #return Rx(longitude)*Ry(latitude)
+        def _project(theta_i, phi_j, w_i, w_j):
             latitude = theta_i
             scale = cos(radians(latitude))
@@ -330 +408 @@
             return active*w_i*w_j
     elif projection == 'guyou':  #define PROJECTION 3  (eventually?)
-        def _rotate(theta_i, phi_j):
+        def _project(theta_i, phi_j, psi):
             from .guyou import guyou_invert
             #latitude, longitude = guyou_invert([theta_i], [phi_j])
             longitude, latitude = guyou_invert([phi_j], [theta_i])
-            return Rx(longitude[0])*Ry(latitude[0])
+            return latitude, longitude, psi
+            #return Rx(longitude[0])*Ry(latitude[0])
         def _weight(theta_i, phi_j, w_i, w_j):
             return w_i*w_j
-    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry
-        def _rotate(theta_i, phi_j):
+    elif projection == 'azimuthal_equidistance':
+        # Note that calculates angles for Rz Ry rather than Rx Ry
+        def _project(theta_i, phi_j, psi):
             latitude = sqrt(theta_i**2 + phi_j**2)
             longitude = degrees(np.arctan2(phi_j, theta_i))
             #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
-            return Rz(longitude)*Ry(latitude)
+            return latitude, longitude, psi-longitude, 'zyz'
+            #R = Rz(longitude)*Ry(latitude)*Rz(psi)
+            #return R_to_xyz(R)
+            #return Rz(longitude)*Ry(latitude)
         def _weight(theta_i, phi_j, w_i, w_j):
             # Weighting for each point comes from the integral:
@@ -376 +459 @@
             return weight*w_i*w_j if latitude < 180 else 0
     elif projection == 'azimuthal_equal_area':
-        def _rotate(theta_i, phi_j):
+        # Note that calculates angles for Rz Ry rather than Rx Ry
+        def _project(theta_i, phi_j, psi):
             radius = min(1, sqrt(theta_i**2 + phi_j**2)/180)
             latitude = 180-degrees(2*np.arccos(radius))
             longitude = degrees(np.arctan2(phi_j, theta_i))
             #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude))
-            return Rz(longitude)*Ry(latitude)
+            return latitude, longitude, psi, "zyz"
+            #R = Rz(longitude)*Ry(latitude)*Rz(psi)
+            #return R_to_xyz(R)
+            #return Rz(longitude)*Ry(latitude)
         def _weight(theta_i, phi_j, w_i, w_j):
             latitude = sqrt(theta_i**2 + phi_j**2)
@@ -389 +476 @@
         raise ValueError("unknown projection %r"%projection)
 
+    return _project, _weight
+
+def R_to_xyz(R):
+    """
+    Return phi, theta, psi Tait-Bryan angles corresponding to the given rotation matrix.
+
+    Extracting Euler Angles from a Rotation Matrix
+    Mike Day, Insomniac Games
+    https://d3cw3dd2w32x2b.cloudfront.net/wp-content/uploads/2012/07/euler-angles1.pdf
+    Based on: Shoemake’s "Euler Angle Conversion", Graphics Gems IV, pp.  222-229
+    """
+    phi = np.arctan2(R[1, 2], R[2, 2])
+    theta = np.arctan2(-R[0, 2], np.sqrt(R[0, 0]**2 + R[0, 1]**2))
+    psi = np.arctan2(R[0, 1], R[0, 0])
+    return np.degrees(phi), np.degrees(theta), np.degrees(psi)
+
+def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian',
+              projection='equirectangular'):
+    """
+    Draw the dispersion mesh showing the theta-phi orientations at which
+    the model will be evaluated.
+    """
+
+    _project, _weight = get_projection(projection)
+    def _rotate(theta, phi, z):
+        dview = _project(theta, phi, 0.)
+        if len(dview) == 4: # hack for zyz coords
+            return Rz(dview[1])*Ry(dview[0])*z
+        else:
+            return Rx(dview[1])*Ry(dview[0])*z
+
+
+    dist_x = np.linspace(-1, 1, n)
+    weights = np.ones_like(dist_x)
+    if dist == 'gaussian':
+        dist_x *= 3
+        weights = exp(-0.5*dist_x**2)
+    elif dist == 'rectangle':
+        # Note: uses sasmodels ridiculous definition of rectangle width
+        dist_x *= sqrt(3)
+    elif dist == 'uniform':
+        pass
+    else:
+        raise ValueError("expected dist to be gaussian, rectangle or uniform")
+
     # mesh in theta, phi formed by rotating z
     dtheta, dphi, dpsi = jitter
     z = np.matrix([[0], [0], [radius]])
-    points = np.hstack([_rotate(theta_i, phi_j)*z
+    points = np.hstack([_rotate(theta_i, phi_j, z)
                         for theta_i in dtheta*dist_x
                         for phi_j in dphi*dist_x])
@@ -470 +602 @@
     Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
     """
-    dtheta, dphi, dpsi = jitter
-    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
+    if jitter is None:
+        return points
+    # Hack to deal with the fact that azimuthal_equidistance uses euler angles
+    if len(jitter) == 4:
+        dtheta, dphi, dpsi, _ = jitter
+        points = Rz(dphi)*Ry(dtheta)*Rz(dpsi)*points
+    else:
+        dtheta, dphi, dpsi = jitter
+        points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
     return points
 
@@ -481 +620 @@
     """
     theta, phi, psi = view
-    points = Rz(phi)*Ry(theta)*Rz(psi)*points
+    points = Rz(phi)*Ry(theta)*Rz(psi)*points # viewing angle
+    #points = Rz(psi)*Ry(pi/2-theta)*Rz(phi)*points # 1-D integration angles
+    #points = Rx(phi)*Ry(theta)*Rz(psi)*points  # angular dispersion angle
     return points
+
+def orient_relative_to_beam_quaternion(view, points):
+    """
+    Apply the view transform to a set of points.
+
+    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
+
+    This variant uses quaternions rather than rotation matrices for the
+    computation.  It works but it is not used because it doesn't solve
+    any problems.  The challenge of mapping theta/phi/psi to SO(3) does
+    not disappear by calculating the transform differently.
+    """
+    theta, phi, psi = view
+    x, y, z = [1, 0, 0], [0, 1, 0], [0, 0, 1]
+    q = Quaternion(1, [0, 0, 0])
+    ## Compose a rotation about the three axes by rotating
+    ## the unit vectors before applying the rotation.
+    #q = Quaternion.from_angle_axis(theta, q.rot(x)) * q
+    #q = Quaternion.from_angle_axis(phi, q.rot(y)) * q
+    #q = Quaternion.from_angle_axis(psi, q.rot(z)) * q
+    ## The above turns out to be equivalent to reversing
+    ## the order of application, so ignore it and use below.
+    q = q * Quaternion.from_angle_axis(theta, x)
+    q = q * Quaternion.from_angle_axis(phi, y)
+    q = q * Quaternion.from_angle_axis(psi, z)
+    ## Reverse the order by post-multiply rather than pre-multiply
+    #q = Quaternion.from_angle_axis(theta, x) * q
+    #q = Quaternion.from_angle_axis(phi, y) * q
+    #q = Quaternion.from_angle_axis(psi, z) * q
+    #print("axes psi", q.rot(np.matrix([x, y, z]).T))
+    return q.rot(points)
+#orient_relative_to_beam = orient_relative_to_beam_quaternion
+
+# Simple stand-alone quaternion class
+import numpy as np
+from copy import copy
+class Quaternion(object):
+    def __init__(self, w, r):
+         self.w = w
+         self.r = np.asarray(r, dtype='d')
+    @staticmethod
+    def from_angle_axis(theta, r):
+         theta = np.radians(theta)/2
+         r = np.asarray(r)
+         w = np.cos(theta)
+         r = np.sin(theta)*r/np.dot(r,r)
+         return Quaternion(w, r)
+    def __mul__(self, other):
+        if isinstance(other, Quaternion):
+            w = self.w*other.w - np.dot(self.r, other.r)
+            r = self.w*other.r + other.w*self.r + np.cross(self.r, other.r)
+            return Quaternion(w, r)
+    def rot(self, v):
+        v = np.asarray(v).T
+        use_transpose = (v.shape[-1] != 3)
+        if use_transpose: v = v.T
+        v = v + np.cross(2*self.r, np.cross(self.r, v) + self.w*v)
+        #v = v + 2*self.w*np.cross(self.r, v) + np.cross(2*self.r, np.cross(self.r, v))
+        if use_transpose: v = v.T
+        return v.T
+    def conj(self):
+        return Quaternion(self.w, -self.r)
+    def inv(self):
+        return self.conj()/self.norm()**2
+    def norm(self):
+        return np.sqrt(self.w**2 + np.sum(self.r**2))
+    def __str__(self):
+        return "%g%+gi%+gj%+gk"%(self.w, self.r[0], self.r[1], self.r[2])
+def test_qrot():
+    # Define rotation of 60 degrees around an axis in y-z that is 60 degrees from y.
+    # The rotation axis is determined by rotating the point [0, 1, 0] about x.
+    ax = Quaternion.from_angle_axis(60, [1, 0, 0]).rot([0, 1, 0])
+    q = Quaternion.from_angle_axis(60, ax)
+    # Set the point to be rotated, and its expected rotated position.
+    p = [1, -1, 2]
+    target = [(10+4*np.sqrt(3))/8, (1+2*np.sqrt(3))/8, (14-3*np.sqrt(3))/8]
+    #print(q, q.rot(p) - target)
+    assert max(abs(q.rot(p) - target)) < 1e-14
+#test_qrot()
+#import sys; sys.exit()
 
 # translate between number of dimension of dispersity and the number of
@@ -556 +777 @@
         vmin = vmax*10**-7
         #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top')
+    #vmin, vmax = Iqxy.min(), Iqxy.max()
     #print("range",(vmin,vmax))
     #qx, qy = np.meshgrid(qx, qy)
     if 0:
+        from matplotlib import cm
         level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i')
         level[level < 0] = 0
         colors = plt.get_cmap()(level)
-        axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors)
+        #colors = cm.coolwarm(level)
+        #colors = cm.gist_yarg(level)
+        #colors = cm.Wistia(level)
+        colors[level<=0, 3] = 0.  # set floor to transparent
+        x, y = np.meshgrid(qx/qx.max(), qy/qy.max())
+        axes.plot_surface(x, y, -1.1*np.ones_like(x), facecolors=colors)
     elif 1:
         axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1,
@@ -692 +920 @@
 }
 
+
 def run(model_name='parallelepiped', size=(10, 40, 100),
+        view=(0, 0, 0), jitter=(0, 0, 0),
         dist='gaussian', mesh=30,
         projection='equirectangular'):
@@ -702 +932 @@
 
     *size* gives the dimensions (a, b, c) of the shape.
+
+    *view* gives the initial view (theta, phi, psi) of the shape.
+
+    *view* gives the initial jitter (dtheta, dphi, dpsi) of the shape.
 
     *dist* is the type of dispersition: gaussian, rectangle, or uniform.
@@ -721 +955 @@
     calculator, size = select_calculator(model_name, n=150, size=size)
     draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped)
+    #draw_shape = draw_fcc
 
     ## uncomment to set an independent the colour range for every view
@@ -726 +961 @@
     calculator.limits = None
 
-    ## initial view
-    #theta, dtheta = 70., 10.
-    #phi, dphi = -45., 3.
-    #psi, dpsi = -45., 3.
-    theta, phi, psi = 0, 0, 0
-    dtheta, dphi, dpsi = 0, 0, 0
-
+    PLOT_ENGINE(calculator, draw_shape, size, view, jitter, dist, mesh, projection)
+
+def mpl_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection):
+    # Note: travis-ci does not support mpl_toolkits.mplot3d, but this shouldn't be
+    # an issue since we are lazy-loading the package on a path that isn't tested.
+    import mpl_toolkits.mplot3d  # Adds projection='3d' option to subplot
+    import matplotlib as mpl
+    import matplotlib.pyplot as plt
+    from matplotlib.widgets import Slider
+    
     ## create the plot window
     #plt.hold(True)
@@ -747 +985 @@
         pass
 
+
     # CRUFT: use axisbg instead of facecolor for matplotlib<2
     facecolor_prop = 'facecolor' if mpl.__version__ > '2' else 'axisbg'
@@ -762 +1001 @@
     axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props)
     axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props)
+
     # Note: using ridiculous definition of rectangle distribution, whose width
     # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep
     # the maximum width to 90.
     dlimit = DIST_LIMITS[dist]
-    sdtheta = Slider(axes_dtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta)
-    sdphi = Slider(axes_dphi, 'dPhi', 0, 2*dlimit, valinit=dphi)
-    sdpsi = Slider(axes_dpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi)
-
+    sdtheta = Slider(axes_dtheta, u'Δξ', 0, 2*dlimit, valinit=0)
+    sdphi = Slider(axes_dphi, u'Δφ', 0, 2*dlimit, valinit=0)
+    sdpsi = Slider(axes_dpsi, u'Δψ', 0, 2*dlimit, valinit=0)
+
+    ## initial view and jitter
+    theta, phi, psi = view
+    stheta.set_val(theta)
+    sphi.set_val(phi)
+    spsi.set_val(psi)
+    dtheta, dphi, dpsi = jitter
+    sdtheta.set_val(dtheta)
+    sdphi.set_val(dphi)
+    sdpsi.set_val(dpsi)
 
     ## callback to draw the new view
@@ -777 +1026 @@
         # set small jitter as 0 if multiple pd dims
         dims = sum(v > 0 for v in jitter)
-        limit = [0, 0.5, 5][dims]
+        limit = [0, 0.5, 5, 5][dims]
         jitter = [0 if v < limit else v for v in jitter]
         axes.cla()
-        draw_beam(axes, (0, 0))
-        draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape)
-        #draw_jitter(axes, view, (0,0,0))
+
+        ## Visualize as person on globe
+        #draw_sphere(axes)
+        #draw_person_on_sphere(axes, view)
+
+        ## Move beam instead of shape
+        #draw_beam(axes, -view[:2])
+        #draw_jitter(axes, (0,0,0), (0,0,0), views=3)
+
+        ## Move shape and draw scattering
+        draw_beam(axes, (0, 0), alpha=1.)
+        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1.,
+                    draw_shape=draw_shape, projection=projection, views=3)
         draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection)
         draw_scattering(calculator, axes, view, jitter, dist=dist)
+
         plt.gcf().canvas.draw()
 
@@ -800 +1060 @@
     ## go interactive
     plt.show()
+
+
+def map_colors(z, kw):
+    from matplotlib import cm
+
+    cmap = kw.pop('cmap', cm.coolwarm)
+    alpha = kw.pop('alpha', None)
+    vmin = kw.pop('vmin', z.min())
+    vmax = kw.pop('vmax', z.max())
+    c = kw.pop('c', None)
+    color = kw.pop('color', c)
+    if color is None:
+        znorm = ((z - vmin) / (vmax - vmin)).clip(0, 1)
+        color = cmap(znorm)
+    elif isinstance(color, np.ndarray) and color.shape == z.shape:
+        color = cmap(color)
+    if alpha is None:
+        if isinstance(color, np.ndarray):
+            color = color[..., :3]
+    else:
+        color[..., 3] = alpha
+    kw['color'] = color
+
+def make_vec(*args):
+    #return [np.asarray(v, 'd').flatten() for v in args]
+    return [np.asarray(v, 'd') for v in args]
+
+def make_image(z, kw):
+    import PIL.Image
+    from matplotlib import cm
+
+    cmap = kw.pop('cmap', cm.coolwarm)
+
+    znorm = (z-z.min())/z.ptp()
+    c = cmap(znorm)
+    c = c[..., :3]
+    rgb = np.asarray(c*255, 'u1')
+    image = PIL.Image.fromarray(rgb, mode='RGB')
+    return image
+
+
+_IPV_MARKERS = {
+    'o': 'sphere',
+}
+_IPV_COLORS = {
+    'w': 'white',
+    'k': 'black',
+    'c': 'cyan',
+    'm': 'magenta',
+    'y': 'yellow',
+    'r': 'red',
+    'g': 'green',
+    'b': 'blue',
+}
+def ipv_fix_color(kw):
+    alpha = kw.pop('alpha', None)
+    color = kw.get('color', None)
+    if isinstance(color, str):
+        color = _IPV_COLORS.get(color, color)
+        kw['color'] = color
+    if alpha is not None:
+        color = kw['color']
+        #TODO: convert color to [r, g, b, a] if not already
+        if isinstance(color, (tuple, list)):
+            if len(color) == 3:
+                color = (color[0], color[1], color[2], alpha)
+            else:
+                color = (color[0], color[1], color[2], alpha*color[3])
+            color = np.array(color)
+        if isinstance(color, np.ndarray) and color.shape[-1] == 4:
+            color[..., 3] = alpha
+            kw['color'] = color
+
+def ipv_set_transparency(kw, obj):
+    color = kw.get('color', None)
+    if (isinstance(color, np.ndarray)
+            and color.shape[-1] == 4
+            and (color[..., 3] != 1.0).any()):
+        obj.material.transparent = True
+        obj.material.side = "FrontSide"
+
+def ipv_axes():
+    import ipyvolume as ipv
+
+    class Plotter:
+        # transparency can be achieved by setting the following:
+        #    mesh.color = [r, g, b, alpha]
+        #    mesh.material.transparent = True
+        #    mesh.material.side = "FrontSide"
+        # smooth(ish) rotation can be achieved by setting:
+        #    slide.continuous_update = True
+        #    figure.animation = 0.
+        #    mesh.material.x = x
+        #    mesh.material.y = y
+        #    mesh.material.z = z
+        # maybe need to synchronize update of x/y/z to avoid shimmy when moving
+        def plot(self, x, y, z, **kw):
+            ipv_fix_color(kw)
+            x, y, z = make_vec(x, y, z)
+            ipv.plot(x, y, z, **kw)
+        def plot_surface(self, x, y, z, **kw):
+            facecolors = kw.pop('facecolors', None)
+            if facecolors is not None:
+                kw['color'] = facecolors
+            ipv_fix_color(kw)
+            x, y, z = make_vec(x, y, z)
+            h = ipv.plot_surface(x, y, z, **kw)
+            ipv_set_transparency(kw, h)
+            #h.material.side = "DoubleSide"
+            return h
+        def plot_trisurf(self, x, y, triangles=None, Z=None, **kw):
+            kw.pop('linewidth', None)
+            ipv_fix_color(kw)
+            x, y, z = make_vec(x, y, Z)
+            if triangles is not None:
+                triangles = np.asarray(triangles)
+            h = ipv.plot_trisurf(x, y, z, triangles=triangles, **kw)
+            ipv_set_transparency(kw, h)
+            return h
+        def scatter(self, x, y, z, **kw):
+            x, y, z = make_vec(x, y, z)
+            map_colors(z, kw)
+            marker = kw.get('marker', None)
+            kw['marker'] = _IPV_MARKERS.get(marker, marker)
+            h = ipv.scatter(x, y, z, **kw)
+            ipv_set_transparency(kw, h)
+            return h
+        def contourf(self, x, y, v, zdir='z', offset=0, levels=None, **kw):
+            # Don't use contour for now (although we might want to later)
+            self.pcolor(x, y, v, zdir='z', offset=offset, **kw)
+        def pcolor(self, x, y, v, zdir='z', offset=0, **kw):
+            x, y, v = make_vec(x, y, v)
+            image = make_image(v, kw)
+            xmin, xmax = x.min(), x.max()
+            ymin, ymax = y.min(), y.max()
+            x = np.array([[xmin, xmax], [xmin, xmax]])
+            y = np.array([[ymin, ymin], [ymax, ymax]])
+            z = x*0 + offset
+            u = np.array([[0., 1], [0, 1]])
+            v = np.array([[0., 0], [1, 1]])
+            h = ipv.plot_mesh(x, y, z, u=u, v=v, texture=image, wireframe=False)
+            ipv_set_transparency(kw, h)
+            h.material.side = "DoubleSide"
+            return h
+        def text(self, *args, **kw):
+            pass
+        def set_xlim(self, limits):
+            ipv.xlim(*limits)
+        def set_ylim(self, limits):
+            ipv.ylim(*limits)
+        def set_zlim(self, limits):
+            ipv.zlim(*limits)
+        def set_axes_on(self):
+            ipv.style.axis_on()
+        def set_axis_off(self):
+            ipv.style.axes_off()
+    return Plotter()
+
+def ipv_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection):
+    import ipywidgets as widgets
+    import ipyvolume as ipv
+
+    axes = ipv_axes()
+
+    def draw(view, jitter):
+        camera = ipv.gcf().camera
+        #print(ipv.gcf().__dict__.keys())
+        #print(dir(ipv.gcf()))
+        ipv.figure(animation=0.)  # no animation when updating object mesh
+
+        # set small jitter as 0 if multiple pd dims
+        dims = sum(v > 0 for v in jitter)
+        limit = [0, 0.5, 5, 5][dims]
+        jitter = [0 if v < limit else v for v in jitter]
+
+        ## Visualize as person on globe
+        #draw_beam(axes, (0, 0))
+        #draw_sphere(axes)
+        #draw_person_on_sphere(axes, view)
+
+        ## Move beam instead of shape
+        #draw_beam(axes, view=(-view[0], -view[1]))
+        #draw_jitter(axes, view=(0,0,0), jitter=(0,0,0))
+
+        ## Move shape and draw scattering
+        draw_beam(axes, (0, 0), steps=25)
+        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1.0,
+                    draw_shape=draw_shape, projection=projection)
+        draw_mesh(axes, view, jitter, dist=dist, n=mesh, radius=0.95,
+                  projection=projection)
+        draw_scattering(calculator, axes, view, jitter, dist=dist)
+
+        draw_axes(axes, origin=(-1, -1, -1.1))
+        ipv.style.box_off()
+        ipv.style.axes_off()
+        ipv.xyzlabel(" ", " ", " ")
+
+        ipv.gcf().camera = camera
+        ipv.show()
+
+
+    trange, prange = (-180., 180., 1.), (-180., 180., 1.)
+    dtrange, dprange = (0., 180., 1.), (0., 180., 1.)
+
+    ## Super simple interfaca, but uses non-ascii variable namese
+    # Ξ φ ψ Δξ Δφ Δψ
+    #def update(**kw):
+    #    view = kw['Ξ'], kw['φ'], kw['ψ']
+    #    jitter = kw['Δξ'], kw['Δφ'], kw['Δψ']
+    #    draw(view, jitter)
+    #widgets.interact(update, Ξ=trange, φ=prange, ψ=prange, Δξ=dtrange, Δφ=dprange, Δψ=dprange)
+
+    def update(theta, phi, psi, dtheta, dphi, dpsi):
+        draw(view=(theta, phi, psi), jitter=(dtheta, dphi, dpsi))
+
+    def slider(name, slice, init=0.):
+        return widgets.FloatSlider(
+            value=init,
+            min=slice[0],
+            max=slice[1],
+            step=slice[2],
+            description=name,
+            disabled=False,
+            #continuous_update=True,
+            continuous_update=False,
+            orientation='horizontal',
+            readout=True,
+            readout_format='.1f',
+            )
+    theta = slider(u'Ξ', trange, view[0])
+    phi = slider(u'φ', prange, view[1])
+    psi = slider(u'ψ', prange, view[2])
+    dtheta = slider(u'Δξ', dtrange, jitter[0])
+    dphi = slider(u'Δφ', dprange, jitter[1])
+    dpsi = slider(u'Δψ', dprange, jitter[2])
+    fields = {
+        'theta': theta, 'phi': phi, 'psi': psi,
+        'dtheta': dtheta, 'dphi': dphi, 'dpsi': dpsi,
+    }
+    ui = widgets.HBox([
+        widgets.VBox([theta, phi, psi]),
+        widgets.VBox([dtheta, dphi, dpsi])
+    ])
+
+    out = widgets.interactive_output(update, fields)
+    display(ui, out)
+
+
+_ENGINES = {
+    "matplotlib": mpl_plot,
+    "mpl": mpl_plot,
+    #"plotly": plotly_plot,
+    "ipvolume": ipv_plot,
+    "ipv": ipv_plot,
+}
+PLOT_ENGINE = _ENGINES["matplotlib"]
+def set_plotter(name):
+    global PLOT_ENGINE
+    PLOT_ENGINE = _ENGINES[name]
 
 def main():
@@ -811 +1330 @@
     parser.add_argument('-s', '--size', type=str, default='10,40,100',
                         help='a,b,c lengths')
+    parser.add_argument('-v', '--view', type=str, default='0,0,0',
+                        help='initial view angles')
+    parser.add_argument('-j', '--jitter', type=str, default='0,0,0',
+                        help='initial angular dispersion')
     parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS,
                         default=DISTRIBUTIONS[0],
@@ -819 +1342 @@
                         help='oriented shape')
     opts = parser.parse_args()
-    size = tuple(int(v) for v in opts.size.split(','))
-    run(opts.shape, size=size,
+    size = tuple(float(v) for v in opts.size.split(','))
+    view = tuple(float(v) for v in opts.view.split(','))
+    jitter = tuple(float(v) for v in opts.jitter.split(','))
+    run(opts.shape, size=size, view=view, jitter=jitter,
         mesh=opts.mesh, dist=opts.distribution,
         projection=opts.projection)

sasmodels/models/pearl_necklace.c

@@ -40 +40 @@
     const double si = sas_sinx_x(q*A_s);
     const double omsi = 1.0 - si;
-    const double pow_si = pow(si, num_pearls);
+    const double pow_si = pown(si, num_pearls);
 
     // form factor for num_pearls
@@ -81 +81 @@
 radius_from_volume(double radius, double edge_sep, double thick_string, double fp_num_pearls)
 {
-    const int num_pearls = (int) fp_num_pearls +0.5;
     const double vol_tot = form_volume(radius, edge_sep, thick_string, fp_num_pearls);
     return cbrt(vol_tot/M_4PI_3);