Update orientation viewer

explore/angles.py

@@ -223,8 +223,8 @@ def print_steps(jitter, jitter_inv, view, view_inv, qc_only):
     )
 
     comment("\n**** qab from qc ****")
-    # The indirect calculation of qab is better than directly c
-    # alculating qab^2 = qa^2 + qb^2 since qc can be computed
+    # The indirect calculation of qab is better than directly
+    # calculating qab^2 = qa^2 + qb^2 since qc can be computed
     # as qc = M31*qx + M32*qy, thus requiring only two elements
     # of the rotation matrix.
     #vprint(qab, sqrt(qa**2 + qb**2), "Direct calculation of qab")

sasmodels/TwoYukawa/CalTYSk.py

@@ -1,8 +1,8 @@
 import numpy as np
-from numpy import pi, mean
+from numpy import mean, pi
 
-from .Ecoefficient import TYCoeff
 from .CalcRealRoot import CalRealRoot
+from .Ecoefficient import TYCoeff
 from .TInvFourier import TInvFourier
 
 # Supplied Q vector must go out to at least this value to calculate g(r).

sasmodels/TwoYukawa/CalcRealRoot.py

@@ -3,6 +3,7 @@
 from .Ecoefficient import TYCoeff
 from .Epoly import make_poly
 
+
 def CalRealRoot(coeff: TYCoeff):
 
     Ecoefficient = coeff.Ecoefficient

sasmodels/TwoYukawa/Ecoefficient.py

@@ -1,7 +1,7 @@
 from typing import Callable, Tuple
 
 import numpy as np
-from numpy import exp, pi, cos, sin, cosh
+from numpy import cos, cosh, exp, pi, sin
 from numpy.typing import NDArray
 
 # CalCoeff.m
@@ -194,17 +194,17 @@ def ABC12(d1, d2):
               Ccd2_22*Cdd1_12*d1*d2**2 -
               Ccd1_21*Cdd2_22*d1*d2**2 -
               Ccd2_22*d2*k1 + Ccd1_21*d1*k2)))/
-        ((d1*((
+        (d1*(
               (-Ccd1_21)*Ccd2_12*d2 +
-                Ccd1_11*Ccd2_22*d2)))))
+                Ccd1_11*Ccd2_22*d2)))
         C2 = ((-((Ccd2_12*d2*
                 (((-Cd1_1)*d1 - Cdd1_11*d1**2 -
                     Cdd1_12*d1*d2 + k1)) -
               Ccd1_11*d1*
                 (((-Cd2_2)*d2 - Cdd2_12*d1*d2 -
                     Cdd2_22*d2**2 + k2))))) /
-        (((-Ccd1_21)*Ccd2_12*d1*d2 +
-            Ccd1_11*Ccd2_22*d1*d2)))
+        ((-Ccd1_21)*Ccd2_12*d1*d2 +
+            Ccd1_11*Ccd2_22*d1*d2))
         return A, B, C1, C2
     self.ABC12 = ABC12
 
@@ -370,10 +370,10 @@ def tau_d21(s):
 
     E1d02 = 12*c1F01*phi*sigma_d01(z1) - 12*c1F01*exp(-z1)*phi*tau_d01(z1)
 
-    E1d11 = (((12*c1F10*phi*sigma_d01(z1) + \
+    E1d11 = (12*c1F10*phi*sigma_d01(z1) + \
           12*c1F01*phi*sigma_d10(z1) - \
           12*c1F10*exp(-z1)*phi*tau_d01(z1) - \
-          12*c1F01*exp(-z1)*phi*tau_d10(z1))))
+          12*c1F01*exp(-z1)*phi*tau_d10(z1))
     E1d12 = (-c1F01 + 12*c1F11*phi*sigma_d01(z1) +
              12*c1F01*phi*sigma_d11(z1) -
              12*c1F11*exp(-z1)*phi*tau_d01(z1) -

sasmodels/TwoYukawa/TFourier.py

@@ -1,6 +1,7 @@
-from numpy import exp, pi, arange, linspace, abs, ceil, log2, interp
+from numpy import abs, arange, ceil, exp, interp, linspace, log2, pi
 from numpy.fft import fft
 
+
 def TFourier(x, y, deltaX):
     """
     Compute the Fourier transform of a function y(x) with sampling interval deltaX.

sasmodels/TwoYukawa/TInvFourier.py

@@ -1,6 +1,7 @@
-from numpy import exp, ceil, log2, pi, arange, interp
+from numpy import arange, ceil, exp, interp, log2, pi
 from numpy.fft import fft
 
+
 def TInvFourier(x, y, deltaX):
     """
     Inverse Fourier transform implementation

sasmodels/jitter.py

@@ -46,6 +46,7 @@
 """
 
 import argparse
+import warnings
 
 import numpy as np
 from numpy import arccos, arctan2, cos, degrees, exp, log, pi, radians, sin, sqrt
@@ -569,8 +570,7 @@ def draw_labels(axes, view, jitter, text):
 
     # TODO: zdir for labels is broken, and labels aren't appearing.
     for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)):
-        zdir = np.asarray(zdir).flatten()
-        axes.text(p[0], p[1], p[2], label, zdir=zdir)
+        axes.text(p[0], p[1], p[2], label, zdir=tuple(zdir))
 
 # Definition of rotation matrices comes from wikipedia:
 #    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
@@ -615,15 +615,17 @@ def apply_jitter(jitter, points):
     Apply the jitter transform to a set of points.
 
     Points are stored in a 3 x n numpy matrix, not a numpy array or tuple.
+
+    Jitter is (dtheta, dphi, dpsi, convention) where convention is "zyz" or "xyz".
+    If convention is missing it defaults to "xyz".
     """
     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
+    convention = jitter[3] if len(jitter) == 4 else "xyz"
+    dtheta, dphi, dpsi = jitter[:3]
+    if convention == "zyz":
         points = Rz(dphi)@Ry(dtheta)@Rz(dpsi)@points
     else:
-        dtheta, dphi, dpsi = jitter
         points = Rx(dphi)@Ry(dtheta)@Rz(dpsi)@points
     return points
 
@@ -808,12 +810,15 @@ def draw_scattering(calculator, axes, view, jitter, dist='gaussian'):
     ## increase pd_n for testing jitter integration rather than simple viz
     #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)]
 
-    pars = dict(
-        theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n,
-        phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n,
-        psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n,
-    )
-    pars.update(calculator.pars)
+    # Some models do not accept theta, phi or psi
+    pars = calculator.pars.copy()
+    if 'theta' in pars:
+        pars.update(theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n)
+    if 'phi' in pars:
+        pars.update(phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n)
+    if 'psi' in pars:
+        pars.update(psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n)
+
 
     # compute the pattern
     qx, qy = calculator.qxqy
@@ -849,7 +854,7 @@ def draw_scattering(calculator, axes, view, jitter, dist='gaussian'):
     else:
         axes.pcolormesh(qx, qy, Iqxy)
 
-def build_model(model_name, n=150, qmax=0.5, **pars):
+def build_model(model_name, n=150, num_angles=3, qmax=0.5, **pars):
     """
     Build a calculator for the given shape.
 
@@ -881,10 +886,15 @@ def build_model(model_name, n=150, qmax=0.5, **pars):
     # stuff the values for non-orientation parameters into the calculator
     calculator.pars = pars.copy()
     calculator.pars.setdefault('background', 1e-3)
+    if num_angles > 0:
+        calculator.pars.update(theta=0, phi=0)
+    if num_angles == 3:
+        calculator.pars.update(psi=0)
+
 
     # fix the data limits so that we can see if the pattern fades
     # under rotation or angular dispersion
-    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars)
+    Iqxy = calculator(**calculator.pars)
     Iqxy = log(Iqxy)
     vmin, vmax = clipped_range(Iqxy, 0.95, mode='top')
     calculator.limits = vmin, vmax+1
@@ -908,52 +918,54 @@ def select_calculator(model_name, n=150, size=(10, 40, 100)):
     d_factor = 0.06  # for paracrystal models
     if model_name == 'sphere':
         calculator = build_model(
-            'sphere', n=n, radius=c)
+            'sphere', n=n, num_angles=0, radius=c)
         a = b = c
     elif model_name == 'sc_paracrystal':
         a = b = c
         dnn = c
         radius = 0.5*c
         calculator = build_model(
-            'sc_paracrystal', n=n,
+            'sc_paracrystal', n=n, num_angles=0,
             dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius,
             background=0)
     elif model_name == 'fcc_paracrystal':
         a = b = c
         # nearest neigbour distance dnn should be 2 radius, but I think the
         # model uses lattice spacing rather than dnn in its calculations
+        warnings.warn("nearest neighbour distance may be wrong here")
         dnn = 0.5*c
         radius = sqrt(2)/4 * c
         calculator = build_model(
-            'fcc_paracrystal', n=n,
+            'fcc_paracrystal', n=n, num_angles=0,
             dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius,
             background=0)
     elif model_name == 'bcc_paracrystal':
         a = b = c
         # nearest neigbour distance dnn should be 2 radius, but I think the
         # model uses lattice spacing rather than dnn in its calculations
+        warnings.warn("nearest neighbour distance may be wrong here")
         dnn = 0.5*c
         radius = sqrt(3)/2 * c
         calculator = build_model(
-            'bcc_paracrystal', n=n,
+            'bcc_paracrystal', n=n, num_angles=0,
             dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius,
             background=0)
     elif model_name == 'cylinder':
         calculator = build_model(
-            'cylinder', n=n, qmax=0.3, radius=b, length=c)
+            'cylinder', n=n, num_angles=2, qmax=0.3, radius=b, length=c)
         a = b
     elif model_name == 'ellipsoid':
         calculator = build_model(
-            'ellipsoid', n=n, qmax=1.0,
+            'ellipsoid', n=n, num_angles=2, qmax=1.0,
             radius_polar=c, radius_equatorial=b)
         a = b
     elif model_name == 'triaxial_ellipsoid':
         calculator = build_model(
-            'triaxial_ellipsoid', n=n, qmax=0.5,
+            'triaxial_ellipsoid', n=n, num_angles=3, qmax=0.5,
             radius_equat_minor=a, radius_equat_major=b, radius_polar=c)
     elif model_name == 'parallelepiped':
         calculator = build_model(
-            'parallelepiped', n=n, length_a=a, length_b=b, length_c=c)
+            'parallelepiped', n=n, num_angles=3, length_a=a, length_b=b, length_c=c)
     else:
         raise ValueError("unknown model %s"%model_name)
 
@@ -1092,6 +1104,9 @@ def _update(val, axis=None):
         limit = [0, 0.5, 5, 5][dims]
         jitter = [0 if v < limit else v for v in jitter]
         axes.cla()
+        axes.set_xlabel('x axis')
+        axes.set_ylabel('y axis')
+        axes.set_zlabel('z axis')
 
         ## Visualize as person on globe
         #draw_sphere(axes, radius=0.5)

sasmodels/models/two_yukawa.py

@@ -94,7 +94,7 @@
 from numpy import inf
 
 # TODO: pep8 says packages and modules should not use camel case
-from sasmodels.TwoYukawa.CalTYSk import CalTYSk, K_MIN, Z_MIN, Z_MIN_DIFF
+from sasmodels.TwoYukawa.CalTYSk import K_MIN, Z_MIN, Z_MIN_DIFF, CalTYSk
 
 # If you want a customized version of two_yukawa as a plugin (for example,
 # because you want to use the high precision polynomial root solver from mpmath)
@@ -104,9 +104,9 @@
 #     https://github.com/SasView/sasmodels/tree/master/sasmodels/models/two_yukawa.py
 #     https://github.com/SasView/sasmodels/tree/master/sasmodels/TwoYukawa
 if 0:
+    import importlib.util
     import sys
     from pathlib import Path
-    import importlib.util
 
     # Remove existing TwoYukawa from sys.modules to force a reload
     remove = [modname for modname in sys.modules if modname.startswith('TwoYukawa.') or modname == 'TwoYukawa']
@@ -121,7 +121,7 @@
     sys.modules['TwoYukawa'] = module
 
     # Override sasmodels library symbols with the local symbols.
-    from TwoYukawa.CalTYSk import CalTYSk, K_MIN, Z_MIN, Z_MIN_DIFF
+    from TwoYukawa.CalTYSk import K_MIN, Z_MIN, Z_MIN_DIFF, CalTYSk
 
 name = "two_yukawa"
 title = "User model for two Yukawa structure factor (S(q))"