Index: doc/ref/index.rst
===================================================================
--- doc/ref/index.rst	(revision c34a31f9923326734c823a9bf9c47652ff0475e7)
+++ doc/ref/index.rst	(revision 9f12fbed66bab9196fec52997bd80899d65bc9a2)
@@ -10,4 +10,5 @@
    refs.rst
    gpu/gpu_computations.rst
+   gpu/opencl_installation.rst
    magnetism/magnetism.rst
    sesans/sans_to_sesans.rst
Index: sasmodels/models/core_shell_parallelepiped.py
===================================================================
--- sasmodels/models/core_shell_parallelepiped.py	(revision 5810f00b5c1f61f1386ab0be011468bf48fbf8c0)
+++ sasmodels/models/core_shell_parallelepiped.py	(revision cb0dc22aeaff717fb8130fd48c3134fa935d7ceb)
@@ -4,6 +4,12 @@
 
 Calculates the form factor for a rectangular solid with a core-shell structure.
-**The thickness and the scattering length density of the shell or "rim"
-can be different on all three (pairs) of faces.**
+The thickness and the scattering length density of the shell or 
+"rim" can be different on each (pair) of faces. However at this time
+the model does **NOT** actually calculate a c face rim despite the presence of
+the parameter.
+
+.. note::
+   This model was originally ported from NIST IGOR macros. However,t is not
+   yet fully understood by the SasView developers and is currently review.
 
 The form factor is normalized by the particle volume $V$ such that
@@ -35,5 +41,6 @@
     V = ABC + 2t_ABC + 2t_BAC + 2t_CAB
 
-**meaning that there are "gaps" at the corners of the solid.**
+**meaning that there are "gaps" at the corners of the solid.**  Again note that
+$t_C = 0$ currently. 
 
 The intensity calculated follows the :ref:`parallelepiped` model, with the
@@ -41,12 +48,4 @@
 amplitudes of the core and shell, in the same manner as a core-shell model.
 
-.. math::
-
-    F_{a}(Q,\alpha,\beta)=
-    \Bigg(\frac{sin(Q(L_A+2t_A)/2sin\alpha sin\beta)}{Q(L_A+2t_A)/2sin\alpha
-    sin\beta)}
-    - \frac{sin(QL_A/2sin\alpha sin\beta)}{QL_A/2sin\alpha sin\beta)} \Bigg)
-    + \frac{sin(QL_B/2sin\alpha sin\beta)}{QL_B/2sin\alpha sin\beta)}
-    + \frac{sin(QL_C/2sin\alpha sin\beta)}{QL_C/2sin\alpha sin\beta)}
 
 .. note::
@@ -93,10 +92,4 @@
     detector plane.
 
-Validation
-----------
-
-The model uses the form factor calculations implemented in a c-library provided
-by the NIST Center for Neutron Research (Kline, 2006).
-
 References
 ----------
@@ -113,6 +106,7 @@
 
 * **Author:** NIST IGOR/DANSE **Date:** pre 2010
-* **Last Modified by:** Paul Butler **Date:** September 30, 2016
-* **Last Reviewed by:** Miguel Gonzales **Date:** March 21, 2016
+* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016
+* **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017
+* **Currently Under review by:** Paul Butler
 """
 
Index: sasmodels/models/fractal.c
===================================================================
--- sasmodels/models/fractal.c	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/fractal.c	(revision 478882278efdd578c67d41a151c1209805308f3f)
@@ -1,3 +1,9 @@
 #define INVALID(p) (p.fractal_dim < 0.0)
+
+ static double
+ form_volume(double radius)
+ {
+     return M_4PI_3 * cube(radius);
+ }
 
 static double
@@ -13,6 +19,6 @@
 
     //calculate P(q) for the spherical subunits
-    const double V = M_4PI_3*cube(radius);
-    const double pq = V * square((sld_block-sld_solvent)*sas_3j1x_x(q*radius));
+    const double pq = square(form_volume(radius) * (sld_block-sld_solvent)
+                      *sas_3j1x_x(q*radius));
 
     // scale to units cm-1 sr-1 (assuming data on absolute scale)
Index: sasmodels/models/fractal.py
===================================================================
--- sasmodels/models/fractal.py	(revision 925ad6edca2969db8fb763eb02beee2c5a5f2872)
+++ sasmodels/models/fractal.py	(revision df89d77375bb6b013466ad838acfb111b7a0ab90)
@@ -20,10 +20,7 @@
 .. math::
 
-    P(q)&= F(qR_0)^2
-
-    F(q)&= \frac{3 (\sin x - x \cos x)}{x^3}
-
-    V_\text{particle} &= \frac{4}{3}\ \pi R_0
-
+    P(q)&= F(qR_0)^2 \\
+    F(q)&= \frac{3 (\sin x - x \cos x)}{x^3} \\
+    V_\text{particle} &= \frac{4}{3}\ \pi R_0 \\
     S(q) &= 1 + \frac{D_f\  \Gamma\!(D_f-1)}{[1+1/(q \xi)^2\  ]^{(D_f -1)/2}}
     \frac{\sin[(D_f-1) \tan^{-1}(q \xi) ]}{(q R_0)^{D_f}}
@@ -32,5 +29,5 @@
 is the fractal dimension, representing the self similarity of the structure. 
 Note that S(q) here goes negative if $D_f$ is too large, and the Gamma function 
-diverges at $D_f$=0 and $D_f$=1.  
+diverges at $D_f=0$ and $D_f=1$.
 
 **Polydispersity on the radius is provided for.**
@@ -47,11 +44,13 @@
 ----------
 
-J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785
+.. [#] J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785
 
-**Author:** NIST IGOR/DANSE **on:** pre 2010
+Authorship and Verification
+----------------------------
 
-**Last Modified by:** Paul Butler **on:** March 20, 2016
-
-**Last Reviewed by:** Paul Butler **on:** March 20, 2016
+* **Author:** NIST IGOR/DANSE **Date:** pre 2010
+* **Converted to sasmodels by:** Paul Butler **Date:** March 19, 2016
+* **Last Modified by:** Paul Butler **Date:** March 12, 2017
+* **Last Reviewed by:** Paul Butler **Date:** March 12, 2017
 
 """
@@ -84,5 +83,5 @@
 parameters = [["volfraction", "", 0.05, [0.0, 1], "",
                "volume fraction of blocks"],
-              ["radius",    "Ang",  5.0, [0.0, inf], "",
+              ["radius",    "Ang",  5.0, [0.0, inf], "volume",
                "radius of particles"],
               ["fractal_dim",      "",  2.0, [0.0, 6.0], "",
Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision b397165ea71fa5026c955689596b9f7baf7bd132)
+++ sasmodels/resolution.py	(revision b32caaba897db3d9a33220584b9dca61f82df589)
@@ -17,9 +17,5 @@
 
 MINIMUM_RESOLUTION = 1e-8
-
-
-# When extrapolating to -q, what is the minimum positive q relative to q_min
-# that we wish to calculate?
-MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION = 0.01
+MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data
 
 class Resolution(object):
@@ -82,6 +78,14 @@
         self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma)
                        if q_calc is None else np.sort(q_calc))
+
+        # Protect against models which are not defined for very low q.  Limit
+        # the smallest q value evaluated (in absolute) to 0.02*min
+        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
+        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff]
+
+        # Build weight matrix from calculated q values
         self.weight_matrix = pinhole_resolution(self.q_calc, self.q,
                                 np.maximum(q_width, MINIMUM_RESOLUTION))
+        self.q_calc = abs(self.q_calc)
 
     def apply(self, theory):
@@ -123,6 +127,14 @@
         self.q_calc = slit_extend_q(q, qx_width, qy_width) \
             if q_calc is None else np.sort(q_calc)
+
+        # Protect against models which are not defined for very low q.  Limit
+        # the smallest q value evaluated (in absolute) to 0.02*min
+        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q)
+        self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff]
+
+        # Build weight matrix from calculated q values
         self.weight_matrix = \
             slit_resolution(self.q_calc, self.q, qx_width, qy_width)
+        self.q_calc = abs(self.q_calc)
 
     def apply(self, theory):
@@ -153,5 +165,5 @@
     # neither trapezoid nor Simpson's rule improved the accuracy.
     edges = bin_edges(q_calc)
-    edges[edges < 0.0] = 0.0 # clip edges below zero
+    #edges[edges < 0.0] = 0.0 # clip edges below zero
     cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :])
     weights = cdf[1:] - cdf[:-1]
@@ -286,5 +298,5 @@
     # The current algorithm is a midpoint rectangle rule.
     q_edges = bin_edges(q_calc) # Note: requires q > 0
-    q_edges[q_edges < 0.0] = 0.0 # clip edges below zero
+    #q_edges[q_edges < 0.0] = 0.0 # clip edges below zero
     weights = np.zeros((len(q), len(q_calc)), 'd')
 
@@ -392,9 +404,8 @@
     interval.
 
-    if *q_min* is zero or less then *q[0]/10* is used instead.
+    Note that extrapolated values may be negative.
     """
     q = np.sort(q)
     if q_min + 2*MINIMUM_RESOLUTION < q[0]:
-        if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION
         n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15
         q_low = np.linspace(q_min, q[0], n_low+1)[:-1]
@@ -448,5 +459,5 @@
         log_delta_q = log(10.) / points_per_decade
     if q_min < q[0]:
-        if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION
+        if q_min < 0: q_min = q[0]*MINIMUM_ABSOLUTE_Q
         n_low = log_delta_q * (log(q[0])-log(q_min))
         q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1]
