Changes in sasmodels/models/parallelepiped.py [3330bb4:34a9e4e] in sasmodels

File:

• sasmodels/models/parallelepiped.py (modified) (11 diffs)

sasmodels/models/parallelepiped.py

@@ -9 +9 @@
 ----------
 
-| This model calculates the scattering from a rectangular parallelepiped
-| (\:numref:`parallelepiped-image`\).
-| If you need to apply polydispersity, see also :ref:`rectangular-prism`.
+ This model calculates the scattering from a rectangular parallelepiped
+ (\:numref:`parallelepiped-image`\).
+ If you need to apply polydispersity, see also :ref:`rectangular-prism`.
 
 .. _parallelepiped-image:
 
+
 .. figure:: img/parallelepiped_geometry.jpg
 
@@ -21 +22 @@
 .. note::
 
-   The edge of the solid must satisfy the condition that $A < B < C$.
-   This requirement is not enforced in the model, so it is up to the
-   user to check this during the analysis.
+The three dimensions of the parallelepiped (strictly here a cuboid) may be given in 
+$any$ size order. To avoid multiple fit solutions, especially
+with Monte-Carlo fit methods, it may be advisable to restrict their ranges. There may 
+be a number of closely similar "best fits", so some trial and error, or fixing of some 
+dimensions at expected values, may help.
 
 The 1D scattering intensity $I(q)$ is calculated as:
@@ -67 +70 @@
     \mu &= qB
 
-
 The scattering intensity per unit volume is returned in units of |cm^-1|.
 
 NB: The 2nd virial coefficient of the parallelepiped is calculated based on
-the averaged effective radius $(=\sqrt{A B / \pi})$ and
+the averaged effective radius, after appropriately sorting the three
+dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and
 length $(= C)$ values, and used as the effective radius for
 $S(q)$ when $P(q) \cdot S(q)$ is applied.
@@ -102 +105 @@
 .. _parallelepiped-orientation:
 
-.. figure:: img/parallelepiped_angle_definition.jpg
-
-    Definition of the angles for oriented parallelepipeds.
-
-.. figure:: img/parallelepiped_angle_projection.jpg
-
-    Examples of the angles for oriented parallelepipeds against the
+.. figure:: img/parallelepiped_angle_definition.png
+
+    Definition of the angles for oriented parallelepiped, shown with $A<B<C$.
+
+.. figure:: img/parallelepiped_angle_projection.png
+
+    Examples of the angles for an oriented parallelepiped against the
     detector plane.
 
+On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will
+appear. These are actually rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, perpendicular to the $a$ x $c$ and $b$ x $c$ faces. 
+(When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) The third orientation distribution, in $\psi$, is 
+about the $c$ axis of the particle, perpendicular to the $a$ x $b$ face. Some experimentation may be required to 
+understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation 
+distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.) 
+
+    
 For a given orientation of the parallelepiped, the 2D form factor is
 calculated as
@@ -116 +127 @@
 .. math::
 
-    P(q_x, q_y) = \left[\frac{\sin(qA\cos\alpha/2)}{(qA\cos\alpha/2)}\right]^2
-                  \left[\frac{\sin(qB\cos\beta/2)}{(qB\cos\beta/2)}\right]^2
-                  \left[\frac{\sin(qC\cos\gamma/2)}{(qC\cos\gamma/2)}\right]^2
+    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2
+                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2
+                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2
 
 with
@@ -154 +165 @@
 angles.
 
-This model is based on form factor calculations implemented in a c-library
-provided by the NIST Center for Neutron Research (Kline, 2006).
 
 References
@@ -163 +172 @@
 
 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
+
+Authorship and Verification
+----------------------------
+
+* **Author:** This model is based on form factor calculations implemented
+    in a c-library provided by the NIST Center for Neutron Research (Kline, 2006).
+* **Last Modified by:**  Paul Kienzle **Date:** April 05, 2017
+* **Last Reviewed by:**  Richard Heenan **Date:** April 06, 2017
+
 """
 
 import numpy as np
-from numpy import pi, inf, sqrt
+from numpy import pi, inf, sqrt, sin, cos
 
 name = "parallelepiped"
@@ -180 +198 @@
             mu = q*B
         V: Volume of the rectangular parallelepiped
-        alpha: angle between the long axis of the 
+        alpha: angle between the long axis of the
             parallelepiped and the q-vector for 1D
 """
@@ -196 +214 @@
               ["length_c", "Ang", 400, [0, inf], "volume",
                "Larger side of the parallelepiped"],
-              ["theta", "degrees", 60, [-inf, inf], "orientation",
-               "In plane angle"],
-              ["phi", "degrees", 60, [-inf, inf], "orientation",
-               "Out of plane angle"],
-              ["psi", "degrees", 60, [-inf, inf], "orientation",
-               "Rotation angle around its own c axis against q plane"],
+              ["theta", "degrees", 60, [-360, 360], "orientation",
+               "c axis to beam angle"],
+              ["phi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about beam"],
+              ["psi", "degrees", 60, [-360, 360], "orientation",
+               "rotation about c axis"],
              ]
 
@@ -208 +226 @@
 def ER(length_a, length_b, length_c):
     """
-        Return effective radius (ER) for P(q)*S(q)
+    Return effective radius (ER) for P(q)*S(q)
     """
-
+    # now that axes can be in any size order, need to sort a,b,c where a~b and c is either much smaller
+    # or much larger
+    abc = np.vstack((length_a, length_b, length_c))
+    abc = np.sort(abc, axis=0)
+    selector = (abc[1] - abc[0]) > (abc[2] - abc[1])
+    length = np.where(selector, abc[0], abc[2])
     # surface average radius (rough approximation)
-    surf_rad = sqrt(length_a * length_b / pi)
-
-    ddd = 0.75 * surf_rad * (2 * surf_rad * length_c + (length_c + surf_rad) * (length_c + pi * surf_rad))
+    radius = np.sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2]) / pi)
+
+    ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius))
     return 0.5 * (ddd) ** (1. / 3.)
 
@@ -230 +253 @@
             phi_pd=10, phi_pd_n=1,
             psi_pd=10, psi_pd_n=10)
-
-qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5)
+# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
+qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
 tests = [[{}, 0.2, 0.17758004974],
          [{}, [0.2], [0.17758004974]],
-         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.00560296014],
-         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.00560296014]],
+         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0089517140475],
+         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0089517140475]],
         ]
 del qx, qy  # not necessary to delete, but cleaner