Index: sasmodels/resolution.py
===================================================================
--- sasmodels/resolution.py	(revision 5258859dcf2d6cbd1ba2ea63ce85ff792c2056bb)
+++ sasmodels/resolution.py	(revision fdc538a09f8fc03aca7f7563f14b6fc5ca0d178f)
@@ -135,5 +135,5 @@
     """
     #print("apply shapes", theory.shape, weight_matrix.shape)
-    Iq = np.dot(theory[None,:], weight_matrix)
+    Iq = np.dot(theory[None, :], weight_matrix)
     #print("result shape",Iq.shape)
     return Iq.flatten()
@@ -153,8 +153,8 @@
     # neither trapezoid nor Simpson's rule improved the accuracy.
     edges = bin_edges(q_calc)
-    edges[edges<0.0] = 0.0 # clip edges below zero
-    G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:] )
+    edges[edges < 0.0] = 0.0 # clip edges below zero
+    G = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :])
     weights = G[1:] - G[:-1]
-    weights /= np.sum(weights, axis=0)[None,:]
+    weights /= np.sum(weights, axis=0)[None, :]
     return weights
 
@@ -287,5 +287,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')
 
@@ -306,10 +306,10 @@
             #print(qi - h, qi + h)
             #print(in_x + abs_x)
-            weights[i,:] = (in_x + abs_x) * np.diff(q_edges) / (2*h)
+            weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h)
         else:
             L = n_height
             for k in range(-L, L+1):
-                weights[i,:] += _q_perp_weights(q_edges, qi+k*h/L, w)
-            weights[i,:] /= 2*L + 1
+                weights[i, :] += _q_perp_weights(q_edges, qi+k*h/L, w)
+            weights[i, :] /= 2*L + 1
 
     return weights.T
@@ -358,5 +358,5 @@
     log-scaled data.
     """
-    if len(x) < 2 or (np.diff(x)<0).any():
+    if len(x) < 2 or (np.diff(x) < 0).any():
         raise ValueError("Expected bins to be an increasing set")
     edges = np.hstack([
@@ -373,13 +373,13 @@
     """
     step = np.diff(q)
-    index = step>max_step
+    index = step > max_step
     if np.any(index):
         inserts = []
-        for q_i,step_i in zip(q[:-1][index],step[index]):
+        for q_i, step_i in zip(q[:-1][index], step[index]):
             n = np.ceil(step_i/max_step)
-            inserts.extend(q_i + np.arange(1,n)*(step_i/n))
+            inserts.extend(q_i + np.arange(1, n)*(step_i/n))
         # Extend a couple of fringes beyond the end of the data
-        inserts.extend(q[-1] + np.arange(1,8)*max_step)
-        q_calc = np.sort(np.hstack((q,inserts)))
+        inserts.extend(q[-1] + np.arange(1, 8)*max_step)
+        q_calc = np.sort(np.hstack((q, inserts)))
     else:
         q_calc = q
@@ -399,10 +399,10 @@
     if q_min < 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
+        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]
     else:
         q_low = []
     if q_max > q[-1]:
-        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1]>q[-2] else 15
+        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15
         q_high = np.linspace(q[-1], q_max, n_high+1)[1:]
     else:
@@ -452,5 +452,5 @@
         if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION
         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]
+        q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1]
     else:
         q_low = []
@@ -489,5 +489,5 @@
     that make it slow to evaluate but give it good accuracy.
     """
-    from scipy.integrate import romberg, dblquad
+    from scipy.integrate import romberg
 
     if any(k not in form.info['defaults'] for k in pars.keys()):
@@ -520,16 +520,15 @@
             result[i] = r/(2*h)
         else:
-            w_grid = np.linspace(0, w, 21)[None,:]
-            h_grid = np.linspace(-h, h, 23)[:,None]
+            w_grid = np.linspace(0, w, 21)[None, :]
+            h_grid = np.linspace(-h, h, 23)[:, None]
             u = sqrt((qi+h_grid)**2 + w_grid**2)
             Iu = np.interp(u, q_calc, Iq)
             #print(np.trapz(Iu, w_grid, axis=1))
-            Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:,0])
+            Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:, 0])
             result[i] = Is / (2*h*w)
-            """
-            r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w,
-                             args=(qi,))
-            result[i] = r/(w*2*h)
-            """
+            # from scipy.integrate import dblquad
+            # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w,
+            #                  args=(qi,))
+            # result[i] = r/(w*2*h)
 
     # r should be [float, ...], but it is [array([float]), array([float]),...]
@@ -553,7 +552,7 @@
 
     _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq)
-    r = [romberg(_fn, max(qi-nsigma*dqi,1e-10*q[0]), qi+nsigma*dqi, args=(qi, dqi),
-                 divmax=100, vec_func=True, tol=0, rtol=1e-8)
-         for qi,dqi in zip(q,q_width)]
+    r = [romberg(_fn, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi,
+                 args=(qi, dqi), divmax=100, vec_func=True, tol=0, rtol=1e-8)
+         for qi, dqi in zip(q, q_width)]
     return np.asarray(r).flatten()
 
@@ -595,6 +594,8 @@
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
-        answer = [ 9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
-                   5.5555, 4.5584, 3.5606, 2.5623, 2.0000 ]
+        answer = [
+            9.0618, 8.6402, 8.1187, 7.1392, 6.1528,
+            5.5555, 4.5584, 3.5606, 2.5623, 2.0000,
+            ]
         np.testing.assert_allclose(output, answer, atol=1e-4)
 
@@ -608,6 +609,8 @@
         theory = 1000*self.Iq(resolution.q_calc**4)
         output = resolution.apply(theory)
-        answer = [ 8.85785, 8.43012, 7.92687, 6.94566, 6.03660,
-                   5.40363, 4.40655, 3.40880, 2.41058, 2.00000 ]
+        answer = [
+            8.85785, 8.43012, 7.92687, 6.94566, 6.03660,
+            5.40363, 4.40655, 3.40880, 2.41058, 2.00000,
+            ]
         np.testing.assert_allclose(output, answer, atol=1e-4)
 
@@ -620,5 +623,7 @@
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
-        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ]
+        answer = [
+            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
+            ]
         np.testing.assert_allclose(output, answer, atol=1e-4)
 
@@ -632,5 +637,7 @@
         theory = self.Iq(resolution.q_calc)
         output = resolution.apply(theory)
-        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ]
+        answer = [
+            11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0,
+            ]
         np.testing.assert_allclose(output, answer, atol=1e-4)
 
@@ -652,7 +659,9 @@
         theory = 12.0-1000.0*resolution.q_calc
         output = resolution.apply(theory)
-        answer = [ 10.44785079, 9.84991299, 8.98101708,
-                  7.99906585, 6.99998311, 6.00001689,
-                  5.00093415, 4.01898292, 3.15008701, 2.55214921]
+        answer = [
+            10.44785079, 9.84991299, 8.98101708,
+            7.99906585, 6.99998311, 6.00001689,
+            5.00093415, 4.01898292, 3.15008701, 2.55214921,
+            ]
         np.testing.assert_allclose(output, answer, atol=1e-8)
 
@@ -783,5 +792,5 @@
             }
         form = load_model('ellipsoid', dtype='double')
-        q = np.logspace(log10(4e-5),log10(2.5e-2), 68)
+        q = np.logspace(log10(4e-5), log10(2.5e-2), 68)
         width, height = 0.117, 0.
         resolution = Slit1D(q, width=width, height=height)
@@ -1071,5 +1080,5 @@
     #h = 0.0277790
     resolution = Slit1D(q, w, h)
-    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w,h))
+    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h))
 
 def demo():
