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

File:

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

sasmodels/models/parallelepiped.py

@@ -15 +15 @@
 .. _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 edge of the solid used to have to satisfy the condition that $A < B < C$.
+   After some improvements to the effective radius calculation, used with
+   an S(Q), it is beleived that this is no longer the case.
 
 The 1D scattering intensity $I(q)$ is calculated as:
@@ -71 +72 @@
 
 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 +104 @@
 .. _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.
 
@@ -116 +118 @@
 .. 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 +156 @@
 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 +163 @@
 
 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 +189 @@
             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
 """
@@ -208 +217 @@
 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 +244 @@
             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