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

Files:

• explore/beta/data_files/richard_test.txt (added)
• explore/beta/data_files/richard_test2.txt (added)
• explore/beta/data_files/richard_test3.txt (added)
• explore/beta/data_files/richard_test4.txt (added)
• explore/beta/data_files/richard_test5.txt (added)
• explore/beta/data_files/richard_test6.txt (added)
• explore/beta/data_files/sasfit_ellipsoid_schulz_IQBD.txt (deleted)
• explore/beta/data_files/sasfit_ellipsoid_schulz_IQD.txt (deleted)
• explore/beta/data_files/sasfit_ellipsoid_schulz_IQSD.txt (deleted)
• explore/beta/data_files/sasfit_sphere_schulz_IQBD.txt (deleted)
• explore/beta/data_files/sasfit_sphere_schulz_IQD.txt (deleted)
• explore/beta/data_files/sasfit_sphere_schulz_IQSD.txt (deleted)
• 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('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'))
+    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'))
     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('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'))
+    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'))
     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
@@ -6 +5 @@
 
 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
@@ -50 +10 @@
 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), alpha=0.5, steps=25):
+def draw_beam(axes, view=(0, 0)):
     """
     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=alpha)
-
+    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8)
+
+    steps = 25
     u = np.linspace(0, 2 * np.pi, steps)
-    v = np.linspace(-1, 1, 2)
+    v = np.linspace(-1, 1, steps)
 
     r = 0.02
@@ -73 +42 @@
     points = Rz(phi)*Ry(theta)*points
     x, y, z = [v.reshape(shape) for v in points]
-    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)
-
+
+    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5)
 
 def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1):
@@ -97 +55 @@
     x, y, z = transform_xyz(view, jitter, x, y, z)
 
-    axes.plot_surface(x, y, z, color='w', alpha=alpha)
+    axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha)
 
     draw_labels(axes, view, jitter, [
@@ -166 +124 @@
     return atoms
 
-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):
+def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1):
     """Draw a parallelepiped."""
     a, b, c = size
@@ -200 +142 @@
 
     x, y, z = transform_xyz(view, jitter, x, y, z)
-    axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha,
-                      linewidth=0)
-
-    # Colour the c+ face of the box.
+    axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha)
+
+    # Draw pink face on 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 0:
-        color = (1, 0.6, 0.6)  # pink
+    if 1:
         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=color, alpha=alpha)
+        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha)
 
     draw_labels(axes, view, jitter, [
@@ -224 +164 @@
     ])
 
-def draw_sphere(axes, radius=0.5, steps=25, center=(0,0,0), color='w', alpha=1.):
+def draw_sphere(axes, radius=10., steps=100):
     """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)) + 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):
+    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):
     """
     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
-    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)
+    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)
     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 = [
@@ -329 +238 @@
     'azimuthal_equal_area',
 ]
-def get_projection(projection):
-
-    """
+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.
+
     jitter projections
     <https://en.wikipedia.org/wiki/List_of_map_projections>
@@ -387 +299 @@
     # 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 _project(theta_i, phi_j, psi):
-            latitude, longitude = theta_i, phi_j
-            return latitude, longitude, psi
-            #return Rx(phi_j)*Ry(theta_i)
+        def _rotate(theta_i, phi_j):
+            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 _project(theta_i, phi_j, psi):
+        def _rotate(theta_i, phi_j):
             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 latitude, longitude, psi
-            #return Rx(longitude)*Ry(latitude)
-        def _project(theta_i, phi_j, w_i, w_j):
+            return Rx(longitude)*Ry(latitude)
+        def _weight(theta_i, phi_j, w_i, w_j):
             latitude = theta_i
             scale = cos(radians(latitude))
@@ -408 +330 @@
             return active*w_i*w_j
     elif projection == 'guyou':  #define PROJECTION 3  (eventually?)
-        def _project(theta_i, phi_j, psi):
+        def _rotate(theta_i, phi_j):
             from .guyou import guyou_invert
             #latitude, longitude = guyou_invert([theta_i], [phi_j])
             longitude, latitude = guyou_invert([phi_j], [theta_i])
-            return latitude, longitude, psi
-            #return Rx(longitude[0])*Ry(latitude[0])
+            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 that calculates angles for Rz Ry rather than Rx Ry
-        def _project(theta_i, phi_j, psi):
+    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry
+        def _rotate(theta_i, phi_j):
             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 latitude, longitude, psi-longitude, 'zyz'
-            #R = Rz(longitude)*Ry(latitude)*Rz(psi)
-            #return R_to_xyz(R)
-            #return Rz(longitude)*Ry(latitude)
+            return Rz(longitude)*Ry(latitude)
         def _weight(theta_i, phi_j, w_i, w_j):
             # Weighting for each point comes from the integral:
@@ -459 +376 @@
             return weight*w_i*w_j if latitude < 180 else 0
     elif projection == 'azimuthal_equal_area':
-        # Note that calculates angles for Rz Ry rather than Rx Ry
-        def _project(theta_i, phi_j, psi):
+        def _rotate(theta_i, phi_j):
             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 latitude, longitude, psi, "zyz"
-            #R = Rz(longitude)*Ry(latitude)*Rz(psi)
-            #return R_to_xyz(R)
-            #return Rz(longitude)*Ry(latitude)
+            return Rz(longitude)*Ry(latitude)
         def _weight(theta_i, phi_j, w_i, w_j):
             latitude = sqrt(theta_i**2 + phi_j**2)
@@ -476 +389 @@
         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])
@@ -602 +470 @@
     Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
     """
-    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
+    dtheta, dphi, dpsi = jitter
+    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points
     return points
 
@@ -620 +481 @@
     """
     theta, phi, psi = view
-    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
+    points = Rz(phi)*Ry(theta)*Rz(psi)*points
     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
@@ -777 +556 @@
         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)
-        #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)
+        axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors)
     elif 1:
         axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1,
@@ -920 +692 @@
 }
 
-
 def run(model_name='parallelepiped', size=(10, 40, 100),
-        view=(0, 0, 0), jitter=(0, 0, 0),
         dist='gaussian', mesh=30,
         projection='equirectangular'):
@@ -932 +702 @@
 
     *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.
@@ -955 +721 @@
     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
@@ -961 +726 @@
     calculator.limits = None
 
-    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
-    
+    ## 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
+
     ## create the plot window
     #plt.hold(True)
@@ -985 +747 @@
         pass
 
-
     # CRUFT: use axisbg instead of facecolor for matplotlib<2
     facecolor_prop = 'facecolor' if mpl.__version__ > '2' else 'axisbg'
@@ -1001 +762 @@
     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, 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)
+    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)
+
 
     ## callback to draw the new view
@@ -1026 +777 @@
         # set small jitter as 0 if multiple pd dims
         dims = sum(v > 0 for v in jitter)
-        limit = [0, 0.5, 5, 5][dims]
+        limit = [0, 0.5, 5][dims]
         jitter = [0 if v < limit else v for v in jitter]
         axes.cla()
-
-        ## 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_beam(axes, (0, 0))
+        draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape)
+        #draw_jitter(axes, view, (0,0,0))
         draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection)
         draw_scattering(calculator, axes, view, jitter, dist=dist)
-
         plt.gcf().canvas.draw()
 
@@ -1060 +800 @@
     ## 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():
@@ -1330 +811 @@
     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],
@@ -1342 +819 @@
                         help='oriented shape')
     opts = parser.parse_args()
-    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,
+    size = tuple(int(v) for v in opts.size.split(','))
+    run(opts.shape, size=size,
         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 = pown(si, num_pearls);
+    const double pow_si = pow(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);