Index: sasmodels/models/hollow_rectangular_prism.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/hollow_rectangular_prism.c	(revision ab2aea81f7fb392e5f24f1c76dd989914d8f7022)
@@ -5,11 +5,11 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness)
 {
-    double b_side = length_a * b2a_ratio;
-    double c_side = length_a * c2a_ratio;
+    double length_b = length_a * b2a_ratio;
+    double length_c = length_a * c2a_ratio;
     double a_core = length_a - 2.0*thickness;
-    double b_core = b_side - 2.0*thickness;
-    double c_core = c_side - 2.0*thickness;
+    double b_core = length_b - 2.0*thickness;
+    double c_core = length_c - 2.0*thickness;
     double vol_core = a_core * b_core * c_core;
-    double vol_total = length_a * b_side * c_side;
+    double vol_total = length_a * length_b * length_c;
     double vol_shell = vol_total - vol_core;
     return vol_shell;
@@ -26,11 +26,13 @@
     double termA1, termA2, termB1, termB2, termC1, termC2;
     
-    double b_side = length_a * b2a_ratio;
-    double c_side = length_a * c2a_ratio;
+    double length_b = length_a * b2a_ratio;
+    double length_c = length_a * c2a_ratio;
     double a_half = 0.5 * length_a;
-    double b_half = 0.5 * b_side;
-    double c_half = 0.5 * c_side;
+    double b_half = 0.5 * length_b;
+    double c_half = 0.5 * length_c;
+    double vol_total = length_a * length_b * length_c;
+    double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness);
 
-   //Integration limits to use in Gaussian quadrature
+    //Integration limits to use in Gaussian quadrature
     double v1a = 0.0;
     double v1b = M_PI_2;  //theta integration limits
@@ -38,22 +40,16 @@
     double v2b = M_PI_2;  //phi integration limits
     
-    //Order of integration
-    int nordi=76;			        
-    int nordj=76;
+    double outer_sum = 0.0;
+    
+    for(int i=0; i<76; i++) {
 
-    double sumi = 0.0;
-    
-    for(int i=0; i<nordi; i++) {
+        double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );    
 
-	    double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );	
+        double termC1 = sinc(q * c_half * cos(theta));
+        double termC2 = sinc(q * (c_half-thickness)*cos(theta));
 
-	    double arg = q * c_half * cos(theta);
-	    if (fabs(arg) > 1.e-16) {termC1 = sin(arg)/arg;} else {termC1 = 1.0;}
-	    arg = q * (c_half-thickness)*cos(theta);
-	    if (fabs(arg) > 1.e-16) {termC2 = sin(arg)/arg;} else {termC2 = 1.0;}
-
-	    double sumj = 0.0;
+        double inner_sum = 0.0;
         
-	    for(int j=0; j<nordj; j++) {
+        for(int j=0; j<76; j++) {
 
             double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
@@ -61,28 +57,23 @@
             // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0
 
-	        arg = q * a_half * sin(theta) * sin(phi);
-	        if (fabs(arg) > 1.e-16) {termA1 = sin(arg)/arg;} else {termA1 = 1.0;}
-	        arg = q * (a_half-thickness) * sin(theta) * sin(phi);
-	        if (fabs(arg) > 1.e-16) {termA2 = sin(arg)/arg;} else {termA2 = 1.0;}
+            termA1 = sinc(q * a_half * sin(theta) * sin(phi));
+            termA2 = sinc(q * (a_half-thickness) * sin(theta) * sin(phi));
 
-	        arg = q * b_half * sin(theta) * cos(phi);
-	        if (fabs(arg) > 1.e-16) {termB1 = sin(arg)/arg;} else {termB1 = 1.0;}
-	        arg = q * (b_half-thickness) * sin(theta) * cos(phi);
-	        if (fabs(arg) > 1.e-16) {termB2 = sin(arg)/arg;} else {termB2 = 1.0;}
+            termB1 = sinc(q * b_half * sin(theta) * cos(phi));
+            termB2 = sinc(q * (b_half-thickness) * sin(theta) * cos(phi));
 
-            double AP1 = (length_a*b_side*c_side) * termA1 * termB1 * termC1;
-            double AP2 = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness) * termA2 * termB2 * termC2;
-            double AP = AP1 - AP2;
+            double AP1 = vol_total * termA1 * termB1 * termC1;
+            double AP2 = vol_core * termA2 * termB2 * termC2;
 
-	        sumj += Gauss76Wt[j] * (AP*AP);
+            inner_sum += Gauss76Wt[j] * square(AP1-AP2);
 
-	    }
+        }
 
-	    sumj = 0.5 * (v2b-v2a) * sumj;
-	    sumi += Gauss76Wt[i] * sumj * sin(theta);
+        inner_sum = 0.5 * (v2b-v2a) * inner_sum;
+        outer_sum += Gauss76Wt[i] * inner_sum * sin(theta);
 
     }
 
-    double answer = 0.5*(v1b-v1a)*sumi;
+    double answer = 0.5*(v1b-v1a)*outer_sum;
 
     // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization)
@@ -91,5 +82,5 @@
 
     // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
-    answer *= (sld-solvent_sld)*(sld-solvent_sld);
+    answer *= square(sld-solvent_sld);
 
     // Convert from [1e-12 A-1] to [cm-1]
Index: sasmodels/models/hollow_rectangular_prism.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/hollow_rectangular_prism.py	(revision ab2aea81f7fb392e5f24f1c76dd989914d8f7022)
@@ -13,5 +13,5 @@
 the difference of the amplitudes of two massive parallelepipeds
 differing in their outermost dimensions in each direction by the
-same length increment :math:`2\Delta` (Nayuk, 2012).
+same length increment $2\Delta$ (Nayuk, 2012).
 
 As in the case of the massive parallelepiped model (:ref:`rectangular-prism`),
@@ -58,9 +58,9 @@
 
 .. math::
-  I(q) = \text{scale} \times V \times (\rho_{\text{p}} -
-  \rho_{\text{solvent}})^2 \times P(q) + \text{background}
+  I(q) = \text{scale} \times V \times (\rho_\text{p} -
+  \rho_\text{solvent})^2 \times P(q) + \text{background}
 
-where $\rho_{\text{p}}$ is the scattering length of the parallelepiped,
-$\rho_{\text{solvent}}$ is the scattering length of the solvent,
+where $\rho_\text{p}$ is the scattering length of the parallelepiped,
+$\rho_\text{solvent}$ is the scattering length of the solvent,
 and (if the data are in absolute units) *scale* represents the volume fraction
 (which is unitless).
@@ -148,5 +148,5 @@
 # parameters for demo
 demo = dict(scale=1, background=0,
-            sld=6.3e-6, sld_solvent=1.0e-6,
+            sld=6.3, sld_solvent=1.0,
             length_a=35, b2a_ratio=1, c2a_ratio=1, thickness=1,
             length_a_pd=0.1, length_a_pd_n=10,
Index: sasmodels/models/hollow_rectangular_prism_thin_walls.c
===================================================================
--- sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.c	(revision ab2aea81f7fb392e5f24f1c76dd989914d8f7022)
@@ -5,7 +5,7 @@
 double form_volume(double length_a, double b2a_ratio, double c2a_ratio)
 {
-    double b_side = length_a * b2a_ratio;
-    double c_side = length_a * c2a_ratio;
-    double vol_shell = 2.0 * (length_a*b_side + length_a*c_side + b_side*c_side);
+    double length_b = length_a * b2a_ratio;
+    double length_c = length_a * c2a_ratio;
+    double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c);
     return vol_shell;
 }
@@ -18,61 +18,63 @@
     double c2a_ratio)
 {
-    double b_side = length_a * b2a_ratio;
-    double c_side = length_a * c2a_ratio;
-    double a_half = 0.5 * length_a;
-    double b_half = 0.5 * b_side;
-    double c_half = 0.5 * c_side;
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double a_half = 0.5 * length_a;
+    const double b_half = 0.5 * length_b;
+    const double c_half = 0.5 * length_c;
 
    //Integration limits to use in Gaussian quadrature
-    double v1a = 0.0;
-    double v1b = M_PI_2;  //theta integration limits
-    double v2a = 0.0;
-    double v2b = M_PI_2;  //phi integration limits
+    const double v1a = 0.0;
+    const double v1b = M_PI_2;  //theta integration limits
+    const double v2a = 0.0;
+    const double v2b = M_PI_2;  //phi integration limits
     
-    //Order of integration
-    int nordi=76;			        
-    int nordj=76;
+    double outer_sum = 0.0;
+    for(int i=0; i<76; i++) {
+        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
 
-    double sumi = 0.0;
-    
-    for(int i=0; i<nordi; i++) {
+        double sin_theta, cos_theta;
+        double sin_c, cos_c;
+        SINCOS(theta, sin_theta, cos_theta);
+        SINCOS(q*c_half*cos_theta, sin_c, cos_c);
 
-	    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));
+        const double termAL_theta = 8.0 * cos_c / (q*q*sin_theta*sin_theta);
+        const double termAT_theta = 8.0 * sin_c / (q*q*sin_theta*cos_theta);
 
-	    double sumj = 0.0;
-        
-	    for(int j=0; j<nordj; j++) {
+        double inner_sum = 0.0;
+        for(int j=0; j<76; j++) {
+            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
 
-            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
-            
+            double sin_phi, cos_phi;
+            double sin_a, cos_a;
+            double sin_b, cos_b;
+            SINCOS(phi, sin_phi, cos_phi);
+            SINCOS(q*a_half*sin_theta*sin_phi, sin_a, cos_a);
+            SINCOS(q*b_half*sin_theta*cos_phi, sin_b, cos_b);
+
             // 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));
+            const double AL = termAL_theta
+                * sin_a*sin_b / (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)) );
+            const double AT = termAT_theta
+                * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi );
 
-            sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT);
+            inner_sum += Gauss76Wt[j] * square(AL+AT);
+        }
 
-	    }
-
-	    sumj = 0.5 * (v2b-v2a) * sumj;
-	    sumi += Gauss76Wt[i] * sumj * sin(theta);
-
+        inner_sum *= 0.5 * (v2b-v2a);
+        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta;
     }
 
-    double answer = 0.5*(v1b-v1a)*sumi;
+    outer_sum *= 0.5*(v1b-v1a);
 
     // 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 /= M_PI_2;
+    double answer = outer_sum/M_PI_2;
 
     // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
-    answer *= (sld-solvent_sld)*(sld-solvent_sld);
+    answer *= square(sld-solvent_sld);
 
     // Convert from [1e-12 A-1] to [cm-1]
@@ -80,4 +82,3 @@
 
     return answer;
-    
 }
Index: sasmodels/models/hollow_rectangular_prism_thin_walls.py
===================================================================
--- sasmodels/models/hollow_rectangular_prism_thin_walls.py	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/hollow_rectangular_prism_thin_walls.py	(revision ab2aea81f7fb392e5f24f1c76dd989914d8f7022)
@@ -3,5 +3,5 @@
 r"""
 
-This model provides the form factor, *P(q)*, for a hollow rectangular
+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.
 
@@ -14,36 +14,37 @@
 
 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
+$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
+
+    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}
+    V &= 2AB + 2AC + 2BC \\
+    A_L(q) &=  8 \times \frac{
+            \sin \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right)
+            \sin \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right)
+            \cos \left( \tfrac{1}{2} q C \cos\theta \right)
+        }{q^2 \, \sin^2\theta \, \sin\phi \cos\phi} \\
+    A_T(q) &=  A_F(q) \times
+      \frac{2\,\sin \left( \tfrac{1}{2} q C \cos\theta \right)}{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) }
+
+  A_F(q) =  4 \frac{ \cos \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right)
+                       \sin \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) }
                      {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) }
+              4 \frac{ \sin \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right)
+                       \cos \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) }
                      {q \, \sin\phi \, \sin\theta}
 
@@ -51,8 +52,9 @@
 
 .. 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}}`
+  I(q) = \text{scale} \times V \times (\rho_\text{p} - \rho_\text{solvent})^2 \times P(q)
+
+where $V$ is the volume of the rectangular prism, $\rho_\text{p}$
+is the scattering length of the parallelepiped, $\rho_\text{solvent}$
 is the scattering length of the solvent, and (if the data are in absolute
 units) *scale* represents the volume fraction (which is unitless).
@@ -127,5 +129,5 @@
 # parameters for demo
 demo = dict(scale=1, background=0,
-            sld=6.3e-6, sld_solvent=1.0e-6,
+            sld=6.3, sld_solvent=1.0,
             length_a=35, b2a_ratio=1, c2a_ratio=1,
             length_a_pd=0.1, length_a_pd_n=10,
Index: sasmodels/models/rectangular_prism.c
===================================================================
--- sasmodels/models/rectangular_prism.c	(revision 3a48772f136eb2795903c63707474e7527b4a325)
+++ sasmodels/models/rectangular_prism.c	(revision ab2aea81f7fb392e5f24f1c76dd989914d8f7022)
@@ -15,58 +15,41 @@
     double c2a_ratio)
 {
-    double termA, termB, termC;
-    
-    double b_side = length_a * b2a_ratio;
-    double c_side = length_a * c2a_ratio;
-    double volume = length_a * b_side * c_side;
-    double a_half = 0.5 * length_a;
-    double b_half = 0.5 * b_side;
-    double c_half = 0.5 * c_side;
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double a_half = 0.5 * length_a;
+    const double b_half = 0.5 * length_b;
+    const double c_half = 0.5 * length_c;
 
    //Integration limits to use in Gaussian quadrature
-    double v1a = 0.0;
-    double v1b = M_PI_2;  //theta integration limits
-    double v2a = 0.0;
-    double v2b = M_PI_2;  //phi integration limits
+    const double v1a = 0.0;
+    const double v1b = M_PI_2;  //theta integration limits
+    const double v2a = 0.0;
+    const double v2b = M_PI_2;  //phi integration limits
     
-    //Order of integration
-    int nordi=76;			        
-    int nordj=76;
+    double outer_sum = 0.0;
+    for(int i=0; i<76; i++) {
+        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );
+        double sin_theta, cos_theta;
+        SINCOS(theta, sin_theta, cos_theta);
 
-    double sumi = 0.0;
-    
-    for(int i=0; i<nordi; i++) {
+        const double termC = sinc(q * c_half * cos_theta);
 
-	    double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );	
+        double inner_sum = 0.0;
+        for(int j=0; j<76; j++) {
+            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b );
+            double sin_phi, cos_phi;
+            SINCOS(phi, sin_phi, cos_phi);
 
-	    double arg = q * c_half * cos(theta);
-	    if (fabs(arg) > 1.e-16) {termC = sin(arg)/arg;} else {termC = 1.0;}  
-
-	    double sumj = 0.0;
-        
-	    for(int j=0; j<nordj; j++) {
-
-            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
-
-	        // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0
-
-	        arg = q * a_half * sin(theta) * sin(phi); 
-	        if (fabs(arg) > 1.e-16) {termA = sin(arg)/arg;} else {termA = 1.0;}
-	       
-	        arg = q * b_half * sin(theta) * cos(phi); 
-	        if (fabs(arg) > 1.e-16) {termB = sin(arg)/arg;} else {termB = 1.0;}	  
-               
-	        double AP = termA * termB * termC;  
-
-	        sumj += Gauss76Wt[j] * (AP*AP);
-
-	    }
-
-	    sumj = 0.5 * (v2b-v2a) * sumj;
-	    sumi += Gauss76Wt[i] * sumj * sin(theta);
-
+            // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0
+            const double termA = sinc(q * a_half * sin_theta * sin_phi);
+            const double termB = sinc(q * b_half * sin_theta * cos_phi);
+            const double AP = termA * termB * termC;
+            inner_sum += Gauss76Wt[j] * AP * AP;
+        }
+        inner_sum = 0.5 * (v2b-v2a) * inner_sum;
+        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta;
     }
 
-    double answer = 0.5*(v1b-v1a)*sumi;
+    double answer = 0.5*(v1b-v1a)*outer_sum;
 
     // Normalize by Pi (Eqn. 16). 
@@ -78,5 +61,6 @@
 
     // Multiply by contrast^2 and volume^2
-    answer *= (sld-solvent_sld)*(sld-solvent_sld)*volume*volume;
+    const double volume = length_a * length_b * length_c;
+    answer *= square((sld-solvent_sld)*volume);
 
     // Convert from [1e-12 A-1] to [cm-1] 
@@ -84,4 +68,3 @@
 
     return answer;
-    
 }
Index: sasmodels/models/rectangular_prism.py
===================================================================
--- sasmodels/models/rectangular_prism.py	(revision 2941abfa5112eb51b9d041f501b1a9d3f1bc5bd4)
+++ sasmodels/models/rectangular_prism.py	(revision ab2aea81f7fb392e5f24f1c76dd989914d8f7022)
@@ -128,5 +128,5 @@
 # parameters for demo
 demo = dict(scale=1, background=0,
-            sld=6.3e-6, sld_solvent=1.0e-6,
+            sld=6.3, sld_solvent=1.0,
             length_a=35, b2a_ratio=1, c2a_ratio=1,
             length_a_pd=0.1, length_a_pd_n=10,
