Index: example/sesans_parameters_css-hs.py
===================================================================
--- example/sesans_parameters_css-hs.py	(revision 3415f4fa7cdcc23ef469cfc7387b5cfc2f132505)
+++ example/sesans_parameters_css-hs.py	(revision ca3abe6810a5385238e57da7a70eeeb0c6377ef7)
@@ -22,8 +22,11 @@
 # Initial parameter values (if other than defaults)
 initial_vals = {
-    "sld_core" : 1.0592,
+    "sld_core" : 1.05,
+    "sld_shell" : 2.88*pen-0.05*(1-pen),
     "sld_solvent" : 2.88,
-    "radius" : 890,
-    "thickness" : 130
+    "radius" : 730,
+    "thickness" : 20,
+    "volfraction" : phi,
+    "scale" : (1-phi)
 }
 
Index: sasmodels/kernelcl.py
===================================================================
--- sasmodels/kernelcl.py	(revision c094758ad0df0ba9f8a4aa5b788a7c4f1bdf249f)
+++ sasmodels/kernelcl.py	(revision 8e0d974c44c6bf60f95c3ad8c308c4d0b1399656)
@@ -426,5 +426,5 @@
         self._need_release = [self.loops_b, self.res_b, self.q_input]
 
-    def __call__(self, fixed_pars, pd_pars, cutoff=1e-5):
+    def __call__(self, fixed_pars, pd_pars, cutoff):
         real = (np.float32 if self.q_input.dtype == generate.F32
                 else np.float64 if self.q_input.dtype == generate.F64
Index: sasmodels/models/_spherepy.py
===================================================================
--- sasmodels/models/_spherepy.py	(revision 49da0794022419abc59514fd36d5fed434d4075b)
+++ sasmodels/models/_spherepy.py	(revision 49da0794022419abc59514fd36d5fed434d4075b)
@@ -0,0 +1,121 @@
+r"""
+For information about polarised and magnetic scattering, click here_.
+
+.. _here: polar_mag_help.html
+
+Definition
+----------
+
+The 1D scattering intensity is calculated in the following way (Guinier, 1955)
+
+.. math::
+
+    I(q) = \frac{\text{scale}}{V} \cdot \left[
+        3V(\Delta\rho) \cdot \frac{\sin(qr) - qr\cos(qr))}{(qr)^3}
+        \right]^2 + \text{background}
+
+where *scale* is a volume fraction, $V$ is the volume of the scatterer,
+$r$ is the radius of the sphere, *background* is the background level and
+*sld* and *sld_solvent* are the scattering length densities (SLDs) of the
+scatterer and the solvent respectively.
+
+Note that if your data is in absolute scale, the *scale* should represent
+the volume fraction (which is unitless) if you have a good fit. If not,
+it should represent the volume fraction times a factor (by which your data
+might need to be rescaled).
+
+The 2D scattering intensity is the same as above, regardless of the
+orientation of $\vec q$.
+
+Validation
+----------
+
+Validation of our code was done by comparing the output of the 1D model
+to the output of the software provided by the NIST (Kline, 2006).
+
+
+References
+----------
+
+A Guinier and G. Fournet, *Small-Angle Scattering of X-Rays*,
+John Wiley and Sons, New York, (1955)
+
+*2013/09/09 and 2014/01/06 - Description reviewed by S King and P Parker.*
+"""
+
+import numpy as np
+from numpy import pi, inf, sin, cos, sqrt, log
+
+name = " _sphere (python)"
+title = "PAK testing ideas for Spheres with uniform scattering length density"
+description = """\
+P(q)=(scale/V)*[3V(sld-sld_solvent)*(sin(qr)-qr cos(qr))
+                /(qr)^3]^2 + background
+    r: radius of sphere
+    V: The volume of the scatter
+    sld: the SLD of the sphere
+    sld_solvent: the SLD of the solvent
+"""
+category = "shape:sphere"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["sld", "1e-6/Ang^2", 1, [-inf, inf], "",
+               "Layer scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["radius", "Ang", 50, [0, inf], "volume",
+               "Sphere radius"],
+             ]
+
+
+def form_volume(radius):
+    return 1.333333333333333 * pi * radius ** 3
+
+def Iq(q, sld, sld_solvent, radius):
+    #print "q",q
+    #print "sld,r",sld,sld_solvent,radius
+    qr = q * radius
+    sn, cn = sin(qr), cos(qr)
+    ## The natural expression for the bessel function is the following:
+    ##     bes = 3 * (sn-qr*cn)/qr**3 if qr>0 else 1
+    ## however, to support vector q values we need to handle the conditional
+    ## as a vector, which we do by first evaluating the full expression
+    ## everywhere, then fixing it up where it is broken.   We should probably
+    ## set numpy to ignore the 0/0 error before we do though...
+    bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line
+    bes[qr == 0] = 1
+    fq = bes * (sld - sld_solvent) * form_volume(radius)
+    return 1.0e-4 * fq ** 2
+Iq.vectorized = True  # Iq accepts an array of q values
+
+def Iqxy(qx, qy, sld, sld_solvent, radius):
+    return Iq(sqrt(qx ** 2 + qy ** 2), sld, sld_solvent, radius)
+Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values
+
+def sesans(z, sld, sld_solvent, radius):
+    """
+    Calculate SESANS-correlation function for a solid sphere.
+
+    Wim Bouwman after formulae Timofei Kruglov J.Appl.Cryst. 2003 article
+    """
+    d = z / radius
+    g = np.zeros_like(z)
+    g[d == 0] = 1.
+    low = ((d > 0) & (d < 2))
+    dlow = d[low]
+    dlow2 = dlow ** 2
+    g[low] = sqrt(1 - dlow2 / 4.) * (1 + dlow2 / 8.) + dlow2 / 2.*(1 - dlow2 / 16.) * log(dlow / (2. + sqrt(4. - dlow2)))
+    return g
+sesans.vectorized = True  # sesans accepts an array of z values
+
+def ER(radius):
+    return radius
+
+# VR defaults to 1.0
+
+demo = dict(scale=1, background=0,
+            sld=6, sld_solvent=1,
+            radius=120,
+            radius_pd=.2, radius_pd_n=45)
+oldname = "SphereModel"
+oldpars = dict(sld='sldSph', sld_solvent='sldSolv', radius='radius')
Index: sasmodels/models/core_shell_sphere.py
===================================================================
--- sasmodels/models/core_shell_sphere.py	(revision 3556ad79872d55bd8cff8b5ca27db3a80f326f98)
+++ sasmodels/models/core_shell_sphere.py	(revision 3f167d94d053f25fc53aeb0e9e4e1b6d74f965d1)
@@ -91,4 +91,5 @@
         @param thickness: shell thickness
     """
+    return (1,1)
     whole = 4.0 * pi / 3.0 * pow((radius + thickness), 3)
     core = 4.0 * pi / 3.0 * radius * radius * radius
Index: sasmodels/models/hollow_rectangular_prism.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
+++ sasmodels/models/hollow_rectangular_prism.py	(revision 3d8283b2c240f54cc6f2667de0f1784412a1c426)
@@ -88,5 +88,5 @@
 title = "Hollow rectangular parallelepiped with uniform scattering length density."
 description = """
-    I(q)= scale*V*(sld - solvent_sld)^2*P(q,theta,phi)+background
+    I(q)= scale*V*(sld - sld_solvent)^2*P(q,theta,phi)+background
         P(q,theta,phi) = (2/pi/V^2) * double integral from 0 to pi/2 of ...
            (AP1-AP2)^2(q)*sin(theta)*dtheta*dphi
@@ -103,5 +103,5 @@
 parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "",
                "Parallelepiped scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
+              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "",
                "Solvent scattering length density"],
               ["a_side", "Ang", 35, [0, inf], "volume",
@@ -148,5 +148,5 @@
 # parameters for demo
 demo = dict(scale=1, background=0,
-            sld=6.3e-6, solvent_sld=1.0e-6,
+            sld=6.3e-6, sld_solvent=1.0e-6,
             a_side=35, b2a_ratio=1, c2a_ratio=1, thickness=1,
             a_side_pd=0.1, a_side_pd_n=10,
@@ -158,5 +158,5 @@
 oldname = 'RectangularHollowPrismModel'
 oldpars = dict(a_side='short_side', b2a_ratio='b2a_ratio', c_side='c2a_ratio',
-               thickness='thickness', sld='sldPipe', solvent_sld='sldSolv')
+               thickness='thickness', sld='sldPipe', sld_solvent='sldSolv')
 
 tests = [[{}, 0.2, 0.76687283098],
Index: smodels/models/hollow_rectangular_prism_infinitely_thin_walls.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.c	(revision deb7ee06caa9916d8c0504bdc15ffbc16df623f6)
+++ 	(revision )
@@ -1,97 +1,0 @@
-double form_volume(double a_side, double b2a_ratio, double c2a_ratio);
-double Iq(double q, double sld, double solvent_sld, double a_side, 
-          double b2a_ratio, double c2a_ratio);
-double Iqxy(double qx, double qy, double sld, double solvent_sld, 
-            double a_side, double b2a_ratio, double c2a_ratio);
-
-double form_volume(double a_side, double b2a_ratio, double c2a_ratio)
-{
-    double b_side = a_side * b2a_ratio;
-    double c_side = a_side * c2a_ratio;
-    double vol_shell = 2.0 * (a_side*b_side + a_side*c_side + b_side*c_side);
-    return vol_shell;
-}
-
-double Iq(double q,
-    double sld,
-    double solvent_sld,
-    double a_side,
-    double b2a_ratio,
-    double c2a_ratio)
-{
-    double b_side = a_side * b2a_ratio;
-    double c_side = a_side * c2a_ratio;
-    double a_half = 0.5 * a_side;
-    double b_half = 0.5 * b_side;
-    double c_half = 0.5 * c_side;
-
-   //Integration limits to use in Gaussian quadrature
-    double v1a = 0.0;
-    double v1b = 0.5 * M_PI;  //theta integration limits
-    double v2a = 0.0;
-    double v2b = 0.5 * M_PI;  //phi integration limits
-    
-    //Order of integration
-    int nordi=76;			        
-    int nordj=76;
-
-    double sumi = 0.0;
-    
-    for(int i=0; i<nordi; i++) {
-
-	    double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );	
-        
-        // To check potential problems if denominator goes to zero here !!!
-        double termAL_theta = 8.0*cos(q*c_half*cos(theta)) / (q*q*sin(theta)*sin(theta));
-        double termAT_theta = 8.0*sin(q*c_half*cos(theta)) / (q*q*sin(theta)*cos(theta));
-
-	    double sumj = 0.0;
-        
-	    for(int j=0; j<nordj; j++) {
-
-            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
-            
-            // Amplitude AL from eqn. (7c)
-            double AL = termAL_theta * sin(q*a_half*sin(theta)*sin(phi)) * 
-                sin(q*b_half*sin(theta)*cos(phi)) / (sin(phi)*cos(phi));
-
-            // Amplitude AT from eqn. (9)
-            double AT = termAT_theta * (  (cos(q*a_half*sin(theta)*sin(phi))*sin(q*b_half*sin(theta)*cos(phi))/cos(phi)) 
-                + (cos(q*b_half*sin(theta)*cos(phi))*sin(q*a_half*sin(theta)*sin(phi))/sin(phi)) );
-
-            sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT);
-
-	    }
-
-	    sumj = 0.5 * (v2b-v2a) * sumj;
-	    sumi += Gauss76Wt[i] * sumj * sin(theta);
-
-    }
-
-    double answer = 0.5*(v1b-v1a)*sumi;
-
-    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization)
-    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
-    answer *= (2.0/M_PI);
-
-    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
-    answer *= (sld-solvent_sld)*(sld-solvent_sld);
-
-    // Convert from [1e-12 A-1] to [cm-1]
-    answer *= 1.0e-4;
-
-    return answer;
-    
-}
-
-double Iqxy(double qx, double qy,
-    double sld,
-    double solvent_sld,
-    double a_side,
-    double b2a_ratio,
-    double c2a_ratio)
-{
-    double q = sqrt(qx*qx + qy*qy);
-    double intensity = Iq(q, sld, solvent_sld, a_side, b2a_ratio, c2a_ratio); 
-    return intensity;    
-}
Index: smodels/models/hollow_rectangular_prism_infinitely_thin_walls.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
+++ 	(revision )
@@ -1,145 +1,0 @@
-# rectangular_prism model
-# Note: model title and parameter table are inserted automatically
-r"""
-
-This model provides the form factor, *P(q)*, for a hollow rectangular
-prism with infinitely thin walls. It computes only the 1D scattering, not the 2D.
-
-
-Definition
-----------
-
-The 1D scattering intensity for this model is calculated according to the
-equations given by Nayuk and Huber (Nayuk, 2012).
-
-Assuming a hollow parallelepiped with infinitely thin walls, edge lengths
-:math:`A \le B \le C` and presenting an orientation with respect to the
-scattering vector given by |theta| and |phi|, where |theta| is the angle
-between the *z* axis and the longest axis of the parallelepiped *C*, and
-|phi| is the angle between the scattering vector (lying in the *xy* plane)
-and the *y* axis, the form factor is given by
-
-.. math::
-  P(q) =  \frac{1}{V^2} \frac{2}{\pi} \int_0^{\frac{\pi}{2}}
-  \int_0^{\frac{\pi}{2}} [A_L(q)+A_T(q)]^2 \sin\theta d\theta d\phi
-
-where
-
-.. math::
-  V = 2AB + 2AC + 2BC
-
-.. math::
-  A_L(q) =  8 \times \frac{ \sin \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)
-                              \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr)
-                              \cos \bigl( q \frac{C}{2} \cos\theta \bigr) }
-                            {q^2 \, \sin^2\theta \, \sin\phi \cos\phi}
-
-.. math::
-  A_T(q) =  A_F(q) \times \frac{2 \, \sin \bigl( q \frac{C}{2} \cos\theta \bigr)}{q \, \cos\theta}
-
-and
-
-.. math::
-  A_F(q) =  4 \frac{ \cos \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)
-                       \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) }
-                     {q \, \cos\phi \, \sin\theta} +
-              4 \frac{ \sin \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)
-                       \cos \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) }
-                     {q \, \sin\phi \, \sin\theta}
-
-The 1D scattering intensity is then calculated as
-
-.. math::
-  I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} - \rho_{\mbox{solvent}})^2 \times P(q)
-
-where *V* is the volume of the rectangular prism, :math:`\rho_{\mbox{p}}`
-is the scattering length of the parallelepiped, :math:`\rho_{\mbox{solvent}}`
-is the scattering length of the solvent, and (if the data are in absolute
-units) *scale* represents the volume fraction (which is unitless).
-
-**The 2D scattering intensity is not computed by this model.**
-
-
-Validation
-----------
-
-Validation of the code was conducted  by qualitatively comparing the output
-of the 1D model to the curves shown in (Nayuk, 2012).
-
-
-References
-----------
-
-R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
-
-"""
-
-from numpy import pi, inf, sqrt
-
-name = "hollow_rectangular_prism_infinitely_thin_walls"
-title = "Hollow rectangular parallelepiped with infinitely thin walls."
-description = """
-    I(q)= scale*V*(sld - solvent_sld)^2*P(q)+background
-        with P(q) being the form factor corresponding to a hollow rectangular
-        parallelepiped with infinitely thin walls.
-"""
-category = "shape:parallelepiped"
-
-#             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "",
-               "Parallelepiped scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
-               "Solvent scattering length density"],
-              ["a_side", "Ang", 35, [0, inf], "volume",
-               "Shorter side of the parallelepiped"],
-              ["b2a_ratio", "Ang", 1, [0, inf], "volume",
-               "Ratio sides b/a"],
-              ["c2a_ratio", "Ang", 1, [0, inf], "volume",
-               "Ratio sides c/a"],
-             ]
-
-source = ["lib/gauss76.c", "hollow_rectangular_prism_infinitely_thin_walls.c"]
-
-def ER(a_side, b2a_ratio, c2a_ratio):
-    """
-        Return equivalent radius (ER)
-    """
-    b_side = a_side * b2a_ratio
-    c_side = a_side * c2a_ratio
-
-    # surface average radius (rough approximation)
-    surf_rad = sqrt(a_side * b_side / pi)
-
-    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
-    return 0.5 * (ddd) ** (1. / 3.)
-
-def VR(a_side, b2a_ratio, c2a_ratio):
-    """
-        Return shell volume and total volume
-    """
-    b_side = a_side * b2a_ratio
-    c_side = a_side * c2a_ratio
-    vol_total = a_side * b_side * c_side
-    vol_shell = 2.0 * (a_side*b_side + a_side*c_side + b_side*c_side)
-    return vol_shell, vol_total
-
-
-# parameters for demo
-demo = dict(scale=1, background=0,
-            sld=6.3e-6, solvent_sld=1.0e-6,
-            a_side=35, b2a_ratio=1, c2a_ratio=1,
-            a_side_pd=0.1, a_side_pd_n=10,
-            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
-            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
-
-# For testing against the old sasview models, include the converted parameter
-# names and the target sasview model name.
-oldname = 'RectangularHollowPrismInfThinWallsModel'
-oldpars = dict(a_side='short_side', b2a_ratio='b2a_ratio', c_side='c2a_ratio',
-               sld='sldPipe', solvent_sld='sldSolv')
-
-tests = [[{}, 0.2, 0.837719188592],
-         [{}, [0.2], [0.837719188592]],
-        ]
-
-
Index: sasmodels/models/hollow_rectangular_prism_thin_walls.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision 3d8283b2c240f54cc6f2667de0f1784412a1c426)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision 3d8283b2c240f54cc6f2667de0f1784412a1c426)
@@ -0,0 +1,97 @@
+double form_volume(double a_side, double b2a_ratio, double c2a_ratio);
+double Iq(double q, double sld, double solvent_sld, double a_side, 
+          double b2a_ratio, double c2a_ratio);
+double Iqxy(double qx, double qy, double sld, double solvent_sld, 
+            double a_side, double b2a_ratio, double c2a_ratio);
+
+double form_volume(double a_side, double b2a_ratio, double c2a_ratio)
+{
+    double b_side = a_side * b2a_ratio;
+    double c_side = a_side * c2a_ratio;
+    double vol_shell = 2.0 * (a_side*b_side + a_side*c_side + b_side*c_side);
+    return vol_shell;
+}
+
+double Iq(double q,
+    double sld,
+    double solvent_sld,
+    double a_side,
+    double b2a_ratio,
+    double c2a_ratio)
+{
+    double b_side = a_side * b2a_ratio;
+    double c_side = a_side * c2a_ratio;
+    double a_half = 0.5 * a_side;
+    double b_half = 0.5 * b_side;
+    double c_half = 0.5 * c_side;
+
+   //Integration limits to use in Gaussian quadrature
+    double v1a = 0.0;
+    double v1b = 0.5 * M_PI;  //theta integration limits
+    double v2a = 0.0;
+    double v2b = 0.5 * M_PI;  //phi integration limits
+    
+    //Order of integration
+    int nordi=76;			        
+    int nordj=76;
+
+    double sumi = 0.0;
+    
+    for(int i=0; i<nordi; i++) {
+
+	    double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );	
+        
+        // To check potential problems if denominator goes to zero here !!!
+        double termAL_theta = 8.0*cos(q*c_half*cos(theta)) / (q*q*sin(theta)*sin(theta));
+        double termAT_theta = 8.0*sin(q*c_half*cos(theta)) / (q*q*sin(theta)*cos(theta));
+
+	    double sumj = 0.0;
+        
+	    for(int j=0; j<nordj; j++) {
+
+            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
+            
+            // Amplitude AL from eqn. (7c)
+            double AL = termAL_theta * sin(q*a_half*sin(theta)*sin(phi)) * 
+                sin(q*b_half*sin(theta)*cos(phi)) / (sin(phi)*cos(phi));
+
+            // Amplitude AT from eqn. (9)
+            double AT = termAT_theta * (  (cos(q*a_half*sin(theta)*sin(phi))*sin(q*b_half*sin(theta)*cos(phi))/cos(phi)) 
+                + (cos(q*b_half*sin(theta)*cos(phi))*sin(q*a_half*sin(theta)*sin(phi))/sin(phi)) );
+
+            sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT);
+
+	    }
+
+	    sumj = 0.5 * (v2b-v2a) * sumj;
+	    sumi += Gauss76Wt[i] * sumj * sin(theta);
+
+    }
+
+    double answer = 0.5*(v1b-v1a)*sumi;
+
+    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization)
+    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
+    answer *= (2.0/M_PI);
+
+    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
+    answer *= (sld-solvent_sld)*(sld-solvent_sld);
+
+    // Convert from [1e-12 A-1] to [cm-1]
+    answer *= 1.0e-4;
+
+    return answer;
+    
+}
+
+double Iqxy(double qx, double qy,
+    double sld,
+    double solvent_sld,
+    double a_side,
+    double b2a_ratio,
+    double c2a_ratio)
+{
+    double q = sqrt(qx*qx + qy*qy);
+    double intensity = Iq(q, sld, solvent_sld, a_side, b2a_ratio, c2a_ratio); 
+    return intensity;    
+}
Index: sasmodels/models/hollow_rectangular_prism_thin_walls.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism_thin_walls.py	(revision 3d8283b2c240f54cc6f2667de0f1784412a1c426)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.py	(revision 3d8283b2c240f54cc6f2667de0f1784412a1c426)
@@ -0,0 +1,145 @@
+# rectangular_prism model
+# Note: model title and parameter table are inserted automatically
+r"""
+
+This model provides the form factor, *P(q)*, for a hollow rectangular
+prism with infinitely thin walls. It computes only the 1D scattering, not the 2D.
+
+
+Definition
+----------
+
+The 1D scattering intensity for this model is calculated according to the
+equations given by Nayuk and Huber (Nayuk, 2012).
+
+Assuming a hollow parallelepiped with infinitely thin walls, edge lengths
+:math:`A \le B \le C` and presenting an orientation with respect to the
+scattering vector given by |theta| and |phi|, where |theta| is the angle
+between the *z* axis and the longest axis of the parallelepiped *C*, and
+|phi| is the angle between the scattering vector (lying in the *xy* plane)
+and the *y* axis, the form factor is given by
+
+.. math::
+  P(q) =  \frac{1}{V^2} \frac{2}{\pi} \int_0^{\frac{\pi}{2}}
+  \int_0^{\frac{\pi}{2}} [A_L(q)+A_T(q)]^2 \sin\theta d\theta d\phi
+
+where
+
+.. math::
+  V = 2AB + 2AC + 2BC
+
+.. math::
+  A_L(q) =  8 \times \frac{ \sin \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)
+                              \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr)
+                              \cos \bigl( q \frac{C}{2} \cos\theta \bigr) }
+                            {q^2 \, \sin^2\theta \, \sin\phi \cos\phi}
+
+.. math::
+  A_T(q) =  A_F(q) \times \frac{2 \, \sin \bigl( q \frac{C}{2} \cos\theta \bigr)}{q \, \cos\theta}
+
+and
+
+.. math::
+  A_F(q) =  4 \frac{ \cos \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)
+                       \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) }
+                     {q \, \cos\phi \, \sin\theta} +
+              4 \frac{ \sin \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)
+                       \cos \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) }
+                     {q \, \sin\phi \, \sin\theta}
+
+The 1D scattering intensity is then calculated as
+
+.. math::
+  I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} - \rho_{\mbox{solvent}})^2 \times P(q)
+
+where *V* is the volume of the rectangular prism, :math:`\rho_{\mbox{p}}`
+is the scattering length of the parallelepiped, :math:`\rho_{\mbox{solvent}}`
+is the scattering length of the solvent, and (if the data are in absolute
+units) *scale* represents the volume fraction (which is unitless).
+
+**The 2D scattering intensity is not computed by this model.**
+
+
+Validation
+----------
+
+Validation of the code was conducted  by qualitatively comparing the output
+of the 1D model to the curves shown in (Nayuk, 2012).
+
+
+References
+----------
+
+R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
+
+"""
+
+from numpy import pi, inf, sqrt
+
+name = "hollow_rectangular_prism_thin_walls"
+title = "Hollow rectangular parallelepiped with thin walls."
+description = """
+    I(q)= scale*V*(sld - sld_solvent)^2*P(q)+background
+        with P(q) being the form factor corresponding to a hollow rectangular
+        parallelepiped with infinitely thin walls.
+"""
+category = "shape:parallelepiped"
+
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "",
+               "Parallelepiped scattering length density"],
+              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "",
+               "Solvent scattering length density"],
+              ["a_side", "Ang", 35, [0, inf], "volume",
+               "Shorter side of the parallelepiped"],
+              ["b2a_ratio", "Ang", 1, [0, inf], "volume",
+               "Ratio sides b/a"],
+              ["c2a_ratio", "Ang", 1, [0, inf], "volume",
+               "Ratio sides c/a"],
+             ]
+
+source = ["lib/gauss76.c", "hollow_rectangular_prism_thin_walls.c"]
+
+def ER(a_side, b2a_ratio, c2a_ratio):
+    """
+        Return equivalent radius (ER)
+    """
+    b_side = a_side * b2a_ratio
+    c_side = a_side * c2a_ratio
+
+    # surface average radius (rough approximation)
+    surf_rad = sqrt(a_side * b_side / pi)
+
+    ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad))
+    return 0.5 * (ddd) ** (1. / 3.)
+
+def VR(a_side, b2a_ratio, c2a_ratio):
+    """
+        Return shell volume and total volume
+    """
+    b_side = a_side * b2a_ratio
+    c_side = a_side * c2a_ratio
+    vol_total = a_side * b_side * c_side
+    vol_shell = 2.0 * (a_side*b_side + a_side*c_side + b_side*c_side)
+    return vol_shell, vol_total
+
+
+# parameters for demo
+demo = dict(scale=1, background=0,
+            sld=6.3e-6, sld_solvent=1.0e-6,
+            a_side=35, b2a_ratio=1, c2a_ratio=1,
+            a_side_pd=0.1, a_side_pd_n=10,
+            b2a_ratio_pd=0.1, b2a_ratio_pd_n=1,
+            c2a_ratio_pd=0.1, c2a_ratio_pd_n=1)
+
+# For testing against the old sasview models, include the converted parameter
+# names and the target sasview model name.
+oldname = 'RectangularHollowPrismInfThinWallsModel'
+oldpars = dict(a_side='short_side', b2a_ratio='b2a_ratio', c_side='c2a_ratio',
+               sld='sldPipe', sld_solvent='sldSolv')
+
+tests = [[{}, 0.2, 0.837719188592],
+         [{}, [0.2], [0.837719188592]],
+        ]
+
+
Index: sasmodels/models/parallelepiped.py
===================================================================
--- sasmodels/models/parallelepiped.py	(revision 6ef42936b03a4811d7838d1ba343eeb957f3f445)
+++ sasmodels/models/parallelepiped.py	(revision 76e50417db49089f63d6c1f3177f458fb6b749ce)
@@ -151,5 +151,7 @@
 ----------
 
-None.
+P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211
+
+R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854
 """
 
@@ -160,5 +162,5 @@
 title = "Rectangular parallelepiped with uniform scattering length density."
 description = """
-    I(q)= scale*V*(sld - solvent_sld)^2*P(q,alpha)+background
+    I(q)= scale*V*(sld - sld_solvent)^2*P(q,alpha)+background
         P(q,alpha) = integral from 0 to 1 of ...
            phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma
@@ -177,5 +179,5 @@
 parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "",
                "Parallelepiped scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
+              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "",
                "Solvent scattering length density"],
               ["a_side", "Ang", 35, [0, inf], "volume",
@@ -210,5 +212,5 @@
 # parameters for demo
 demo = dict(scale=1, background=0,
-            sld=6.3e-6, solvent_sld=1.0e-6,
+            sld=6.3e-6, sld_solvent=1.0e-6,
             a_side=35, b_side=75, c_side=400,
             theta=45, phi=30, psi=15,
@@ -225,5 +227,5 @@
 oldpars = dict(theta='parallel_theta', phi='parallel_phi', psi='parallel_psi',
                a_side='short_a', b_side='short_b', c_side='long_c',
-               sld='sldPipe', solvent_sld='sldSolv')
+               sld='sldPipe', sld_solvent='sldSolv')
 
 
Index: sasmodels/models/rectangular_prism.py
===================================================================
--- sasmodels/models/rectangular_prism.py	(revision 43b7eea5794b4eb3362114ebde67ff2fb1c38fdc)
+++ sasmodels/models/rectangular_prism.py	(revision 76e50417db49089f63d6c1f3177f458fb6b749ce)
@@ -9,5 +9,5 @@
 The only difference is that the way the relevant
 parameters are defined here (*a*, *b/a*, *c/a* instead of *a*, *b*, *c*)
-allows to use polydispersity with this model while keeping the shape of
+which allows use of polydispersity with this model while keeping the shape of
 the prism (e.g. setting *b/a* = 1 and *c/a* = 1 and applying polydispersity
 to *a* will generate a distribution of cubes of different sizes).
@@ -86,5 +86,5 @@
 title = "Rectangular parallelepiped with uniform scattering length density."
 description = """
-    I(q)= scale*V*(sld - solvent_sld)^2*P(q,theta,phi)+background
+    I(q)= scale*V*(sld - sld_solvent)^2*P(q,theta,phi)+background
         P(q,theta,phi) = (2/pi) * double integral from 0 to pi/2 of ...
            AP^2(q)*sin(theta)*dtheta*dphi
@@ -97,5 +97,5 @@
 parameters = [["sld", "1e-6/Ang^2", 6.3, [-inf, inf], "",
                "Parallelepiped scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "",
+              ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "",
                "Solvent scattering length density"],
               ["a_side", "Ang", 35, [0, inf], "volume",
@@ -125,5 +125,5 @@
 # parameters for demo
 demo = dict(scale=1, background=0,
-            sld=6.3e-6, solvent_sld=1.0e-6,
+            sld=6.3e-6, sld_solvent=1.0e-6,
             a_side=35, b2a_ratio=1, c2a_ratio=1,
             a_side_pd=0.1, a_side_pd_n=10,
@@ -135,5 +135,5 @@
 oldname = 'RectangularPrismModel'
 oldpars = dict(a_side='short_side', b2a_ratio='b2a_ratio', c_side='c2a_ratio',
-               sld='sldPipe', solvent_sld='sldSolv')
+               sld='sldPipe', sld_solvent='sldSolv')
 
 tests = [[{}, 0.2, 0.375248406825],
Index: sasmodels/models/sphere.py
===================================================================
--- sasmodels/models/sphere.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ sasmodels/models/sphere.py	(revision 49da0794022419abc59514fd36d5fed434d4075b)
@@ -17,5 +17,5 @@
 where *scale* is a volume fraction, $V$ is the volume of the scatterer,
 $r$ is the radius of the sphere, *background* is the background level and
-*sld* and *solvent_sld* are the scattering length densities (SLDs) of the
+*sld* and *sld_solvent* are the scattering length densities (SLDs) of the
 scatterer and the solvent respectively.
 
@@ -49,10 +49,10 @@
 title = "Spheres with uniform scattering length density"
 description = """\
-P(q)=(scale/V)*[3V(sld-solvent_sld)*(sin(qr)-qr cos(qr))
+P(q)=(scale/V)*[3V(sld-sld_solvent)*(sin(qr)-qr cos(qr))
                 /(qr)^3]^2 + background
     r: radius of sphere
     V: The volume of the scatter
     sld: the SLD of the sphere
-    solvent_sld: the SLD of the solvent
+    sld_solvent: the SLD of the solvent
 """
 category = "shape:sphere"
@@ -61,5 +61,5 @@
 parameters = [["sld", "1e-6/Ang^2", 1, [-inf, inf], "",
                "Layer scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 6, [-inf, inf], "",
+              ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "",
                "Solvent scattering length density"],
               ["radius", "Ang", 50, [0, inf], "volume",
@@ -76,5 +76,5 @@
 
 Iq = """
-    return sphere_form(q, radius, sld, solvent_sld);
+    return sphere_form(q, radius, sld, sld_solvent);
     """
 
@@ -82,5 +82,5 @@
     // never called since no orientation or magnetic parameters.
     //return -1.0;
-    return Iq(sqrt(qx*qx + qy*qy), sld, solvent_sld, radius);
+    return Iq(sqrt(qx*qx + qy*qy), sld, sld_solvent, radius);
     """
 
@@ -94,7 +94,7 @@
 
 demo = dict(scale=1, background=0,
-            sld=6, solvent_sld=1,
+            sld=6, sld_solvent=1,
             radius=120,
             radius_pd=.2, radius_pd_n=45)
 oldname = "SphereModel"
-oldpars = dict(sld='sldSph', solvent_sld='sldSolv', radius='radius')
+oldpars = dict(sld='sldSph', sld_solvent='sldSolv', radius='radius')
Index: smodels/models/spherepy.py
===================================================================
--- sasmodels/models/spherepy.py	(revision aa2edb2b2fb550caed0506dcf63bc7eac7e1779a)
+++ 	(revision )
@@ -1,121 +1,0 @@
-r"""
-For information about polarised and magnetic scattering, click here_.
-
-.. _here: polar_mag_help.html
-
-Definition
-----------
-
-The 1D scattering intensity is calculated in the following way (Guinier, 1955)
-
-.. math::
-
-    I(q) = \frac{\text{scale}}{V} \cdot \left[
-        3V(\Delta\rho) \cdot \frac{\sin(qr) - qr\cos(qr))}{(qr)^3}
-        \right]^2 + \text{background}
-
-where *scale* is a volume fraction, $V$ is the volume of the scatterer,
-$r$ is the radius of the sphere, *background* is the background level and
-*sld* and *solvent_sld* are the scattering length densities (SLDs) of the
-scatterer and the solvent respectively.
-
-Note that if your data is in absolute scale, the *scale* should represent
-the volume fraction (which is unitless) if you have a good fit. If not,
-it should represent the volume fraction times a factor (by which your data
-might need to be rescaled).
-
-The 2D scattering intensity is the same as above, regardless of the
-orientation of $\vec q$.
-
-Validation
-----------
-
-Validation of our code was done by comparing the output of the 1D model
-to the output of the software provided by the NIST (Kline, 2006).
-
-
-References
-----------
-
-A Guinier and G. Fournet, *Small-Angle Scattering of X-Rays*,
-John Wiley and Sons, New York, (1955)
-
-*2013/09/09 and 2014/01/06 - Description reviewed by S King and P Parker.*
-"""
-
-import numpy as np
-from numpy import pi, inf, sin, cos, sqrt, log
-
-name = "sphere (python)"
-title = "Spheres with uniform scattering length density"
-description = """\
-P(q)=(scale/V)*[3V(sld-solvent_sld)*(sin(qr)-qr cos(qr))
-                /(qr)^3]^2 + background
-    r: radius of sphere
-    V: The volume of the scatter
-    sld: the SLD of the sphere
-    solvent_sld: the SLD of the solvent
-"""
-category = "shape:sphere"
-
-#             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["sld", "1e-6/Ang^2", 1, [-inf, inf], "",
-               "Layer scattering length density"],
-              ["solvent_sld", "1e-6/Ang^2", 6, [-inf, inf], "",
-               "Solvent scattering length density"],
-              ["radius", "Ang", 50, [0, inf], "volume",
-               "Sphere radius"],
-             ]
-
-
-def form_volume(radius):
-    return 1.333333333333333 * pi * radius ** 3
-
-def Iq(q, sld, solvent_sld, radius):
-    #print "q",q
-    #print "sld,r",sld,solvent_sld,radius
-    qr = q * radius
-    sn, cn = sin(qr), cos(qr)
-    ## The natural expression for the bessel function is the following:
-    ##     bes = 3 * (sn-qr*cn)/qr**3 if qr>0 else 1
-    ## however, to support vector q values we need to handle the conditional
-    ## as a vector, which we do by first evaluating the full expression
-    ## everywhere, then fixing it up where it is broken.   We should probably
-    ## set numpy to ignore the 0/0 error before we do though...
-    bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line
-    bes[qr == 0] = 1
-    fq = bes * (sld - solvent_sld) * form_volume(radius)
-    return 1.0e-4 * fq ** 2
-Iq.vectorized = True  # Iq accepts an array of q values
-
-def Iqxy(qx, qy, sld, solvent_sld, radius):
-    return Iq(sqrt(qx ** 2 + qy ** 2), sld, solvent_sld, radius)
-Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values
-
-def sesans(z, sld, solvent_sld, radius):
-    """
-    Calculate SESANS-correlation function for a solid sphere.
-
-    Wim Bouwman after formulae Timofei Kruglov J.Appl.Cryst. 2003 article
-    """
-    d = z / radius
-    g = np.zeros_like(z)
-    g[d == 0] = 1.
-    low = ((d > 0) & (d < 2))
-    dlow = d[low]
-    dlow2 = dlow ** 2
-    g[low] = sqrt(1 - dlow2 / 4.) * (1 + dlow2 / 8.) + dlow2 / 2.*(1 - dlow2 / 16.) * log(dlow / (2. + sqrt(4. - dlow2)))
-    return g
-sesans.vectorized = True  # sesans accepts an array of z values
-
-def ER(radius):
-    return radius
-
-# VR defaults to 1.0
-
-demo = dict(scale=1, background=0,
-            sld=6, solvent_sld=1,
-            radius=120,
-            radius_pd=.2, radius_pd_n=45)
-oldname = "SphereModel"
-oldpars = dict(sld='sldSph', solvent_sld='sldSolv', radius='radius')
Index: sasmodels/models/triaxial_ellipsoid.py
===================================================================
--- sasmodels/models/triaxial_ellipsoid.py	(revision 6ef42936b03a4811d7838d1ba343eeb957f3f445)
+++ sasmodels/models/triaxial_ellipsoid.py	(revision 49da0794022419abc59514fd36d5fed434d4075b)
@@ -91,5 +91,5 @@
                "Solvent scattering length density"],
               ["req_minor", "Ang", 20, [0, inf], "volume",
-               "Minor equitorial radius"],
+               "Minor equatorial radius"],
               ["req_major", "Ang", 400, [0, inf], "volume",
                "Major equatorial radius"],
Index: sasmodels/product.py
===================================================================
--- sasmodels/product.py	(revision 72a081dea6caa8fa493057133a6232f6b118c89a)
+++ sasmodels/product.py	(revision 8e0d974c44c6bf60f95c3ad8c308c4d0b1399656)
@@ -163,6 +163,6 @@
         s_pd[0] = [effect_radius], [1.0]
 
-        p_res = self.p_kernel(p_fixed, p_pd)
-        s_res = self.s_kernel(s_fixed, s_pd)
+        p_res = self.p_kernel(p_fixed, p_pd, cutoff)
+        s_res = self.s_kernel(s_fixed, s_pd, cutoff)
         #print s_fixed, s_pd, p_fixed, p_pd
 
Index: sasmodels/sesans.py
===================================================================
--- sasmodels/sesans.py	(revision c94577f038b813ab23072f3ec552add5da36539e)
+++ sasmodels/sesans.py	(revision 2684d455a49939d488bff52e8a235d181fcbb4be)
@@ -29,5 +29,8 @@
 
     q_min = dq = 0.1 * 2*pi / Rmax
-    return np.arange(q_min, Converter("1/A")(q_max[0], units=q_max[1]), dq)
+    return np.arange(q_min,
+                     Converter(q_max[1])(q_max[0],
+                                         units="1/A"),
+                     dq)
     
 def make_all_q(data):
@@ -75,6 +78,8 @@
 
 def call_hankel(data, q_calc, Iq_calc):
-    return hankel(data.x, data.lam * 1e-9,
-                  data.sample.thickness / 10,
+    return hankel((data.x, data.x_unit),
+                  (data.lam, data.lam_unit),
+                  (data.sample.thickness,
+                   data.sample.thickness_unit),
                   q_calc, Iq_calc)
   
@@ -182,4 +187,11 @@
     *I* [cm$^{-1}$] is the value of the SANS model at *q*
     """
+
+    from sas.sascalc.data_util.nxsunit import Converter
+    wavelength = Converter(wavelength[1])(wavelength[0],"A")
+    thickness = Converter(thickness[1])(thickness[0],"A")
+    Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters
+    SElength = Converter(SElength[1])(SElength[0],"A")
+
     G = np.zeros_like(SElength, 'd')
 #==============================================================================
@@ -192,9 +204,9 @@
 
     # [m^-1] step size in q, needed for integration
-    dq = (q[1]-q[0])*1e10
+    dq = (q[1]-q[0])
 
     # integration step, convert q into [m**-1] and 2 pi circle integration
-    G *= dq*1e10*2*pi
-    G0 *= dq*1e10*2*pi
+    G *= dq*2*pi
+    G0 = np.sum(Iq*q)*dq*2*np.pi
 
     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0))
