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;
-    
 }
