diff --git a/wafo/integrate.py b/wafo/integrate.py
index c69382a..a92e9bd 100644
--- a/wafo/integrate.py
+++ b/wafo/integrate.py
@@ -7,14 +7,14 @@ import scipy.special.orthogonal as ort
 from scipy import special as sp
 
 from scipy.integrate import simps, trapz
-from .plotbackend import plotbackend as plt
-from .demos import humps
-from .misc import dea3
-from .dctpack import dct
+from wafo.plotbackend import plotbackend as plt
+from wafo.demos import humps
+from wafo.misc import dea3
+from wafo.dctpack import dct
 # from pychebfun import Chebfun
 
 _EPS = np.finfo(float).eps
-_POINTS_AND_WEIGHTS = {}
+_NODES_AND_WEIGHTS = {}
 
 __all__ = ['dea3', 'clencurt', 'romberg',
            'h_roots', 'j_roots', 'la_roots', 'p_roots', 'qrule',
@@ -31,7 +31,7 @@ def _assert_warn(cond, msg):
         warnings.warn(msg)
 
 
-def clencurt(fun, a, b, n0=5, trace=False):
+def clencurt(fun, a, b, n=5, trace=False):
     """
     Numerical evaluation of an integral, Clenshaw-Curtis method.
 
@@ -78,31 +78,27 @@ def clencurt(fun, a, b, n0=5, trace=False):
     [2] Clenshaw, C.W. and Curtis, A.R. (1960),
     Numerische Matematik, Vol. 2, pp. 197--205
     """
-
-    # make sure n is even
-    n = 2 * int(n0)
+    # make sure n_2 is even
+    n_2 = 2 * int(n)
     a, b = np.atleast_1d(a, b)
     a_shape = a.shape
-    af = a.ravel()
-    bf = b.ravel()
+    a = a.ravel()
+    b = b.ravel()
 
-    na = np.prod(a_shape)
+    a_size = np.prod(a_shape)
 
-    s = np.r_[0:n + 1]
-    s2 = np.r_[0:n + 1:2]
-    s2.shape = (-1, 1)
-    x1 = np.cos(np.pi * s / n)
-    x1.shape = (-1, 1)
-    x = x1 * (bf - af) / 2. + (bf + af) / 2
+    s = np.c_[0:n_2 + 1:1]
+    s_2 = np.c_[0:n_2 + 1:2]
+    x = np.cos(np.pi * s / n_2) * (b - a) / 2. + (b + a) / 2
 
     if hasattr(fun, '__call__'):
         f = fun(x)
     else:
-        x0 = np.flipud(fun[:, 0])
-        n = len(x0) - 1
-        if abs(x - x0) > 1e-8:
+        x_0 = np.flipud(fun[:, 0])
+        n_2 = len(x_0) - 1
+        if abs(x - x_0) > 1e-8:
             raise ValueError(
-                'Input vector x must equal cos(pi*s/n)*(b-a)/2+(b+a)/2')
+                'Input vector x must equal cos(pi*s/n_2)*(b-a)/2+(b+a)/2')
 
         f = np.flipud(fun[:, 1::])
 
@@ -112,27 +108,27 @@ def clencurt(fun, a, b, n0=5, trace=False):
     # using a Gauss-Lobatto variant, i.e., first and last
     # term f(a) and f(b) is multiplied with 0.5
     f[0, :] = f[0, :] / 2
-    f[n, :] = f[n, :] / 2
+    f[n_2, :] = f[n_2, :] / 2
 
-    # x = cos(pi*0:n/n)
+    # x = cos(pi*0:n_2/n_2)
     # f = f(x)
     #
     #               N+1
     #  c(k) = (2/N) sum  f''(n)*cos(pi*(2*k-2)*(n-1)/N), 1 <= k <= N/2+1.
     #               n=1
-    n0 = n // 2
+    n = n_2 // 2
     fft = np.fft.fft
-    tmp = np.real(fft(f[:n, :], axis=0))
-    c = 2 / n * (tmp[0:n0 + 1, :] + np.cos(np.pi * s2) * f[n, :])
+    tmp = np.real(fft(f[:n_2, :], axis=0))
+    c = 2 / n_2 * (tmp[0:n + 1, :] + np.cos(np.pi * s_2) * f[n_2, :])
     c[0, :] = c[0, :] / 2
-    c[n0, :] = c[n0, :] / 2
+    c[n, :] = c[n, :] / 2
 
-    c = c[0:n0 + 1, :] / ((s2 - 1) * (s2 + 1))
-    q_val = (af - bf) * np.sum(c, axis=0)
+    c = c[0:n + 1, :] / ((s_2 - 1) * (s_2 + 1))
+    q_val = (a - b) * np.sum(c, axis=0)
 
-    abserr = (bf - af) * np.abs(c[n0, :])
+    abserr = (b - a) * np.abs(c[n, :])
 
-    if na > 1:
+    if a_size > 1:
         abserr = np.reshape(abserr, a_shape)
         q_val = np.reshape(q_val, a_shape)
     return q_val, abserr
@@ -181,51 +177,54 @@ def romberg(fun, a, b, releps=1e-3, abseps=1e-3):
 
     rom[0, 0] = h * (fun(a) + fun(b)) / 2
     ipower = 1
-    fp = ones(table_limit) * 4
+    f_p = ones(table_limit) * 4
 
-    # ih1 = 0
-    ih2 = 0.
-    ih4 = rom[0, 0]
-    abserr = ih4
+    # q_val1 = 0
+    q_val2 = 0.
+    q_val4 = rom[0, 0]
+    abserr = q_val4
     # epstab = zeros(1,decdigs+7)
     # newflg = 1
-    # [res,abserr,epstab,newflg] = dea(newflg,ih4,abserr,epstab)
+    # [res,abserr,epstab,newflg] = dea(newflg,q_val4,abserr,epstab)
     two = 1
     one = 0
+    converged = False
     for i in range(1, table_limit):
         h *= 0.5
-        un5 = np.sum(fun(a + np.arange(1, 2 * ipower, 2) * h)) * h
+        u_n5 = np.sum(fun(a + np.arange(1, 2 * ipower, 2) * h)) * h
 
         #     trapezoidal approximations
-        # T2n = 0.5 * (Tn + Un) = 0.5*Tn + un5
-        rom[two, 0] = 0.5 * rom[one, 0] + un5
+        # T2n = 0.5 * (Tn + Un) = 0.5*Tn + u_n5
+        rom[two, 0] = 0.5 * rom[one, 0] + u_n5
 
-        fp[i] = 4 * fp[i - 1]
+        f_p[i] = 4 * f_p[i - 1]
         #   Richardson extrapolation
         for k in range(i):
             rom[two, k + 1] = (rom[two, k] +
-                               (rom[two, k] - rom[one, k]) / (fp[k] - 1))
+                               (rom[two, k] - rom[one, k]) / (f_p[k] - 1))
 
-        ih1 = ih2
-        ih2 = ih4
-        ih4 = rom[two, i]
+        q_val1 = q_val2
+        q_val2 = q_val4
+        q_val4 = rom[two, i]
 
         if 2 <= i:
-            res, abserr = dea3(ih1, ih2, ih4)
-            # ih4 = res
-            if abserr <= max(abseps, releps * abs(res)):
+            res, abserr = dea3(q_val1, q_val2, q_val4)
+            # q_val4 = res
+            converged = abserr <= max(abseps, releps * abs(res))
+            if converged:
                 break
 
         # rom(1,1:i) = rom(2,1:i)
         two = one
         one = (one + 1) % 2
         ipower *= 2
+    _assert(converged, "Integral did not converge to the required accuracy!")
     return res, abserr
 
 
 def _h_roots_newton(n, releps=3e-14, max_iter=10):
-    # PIM4=0.7511255444649425
-    PIM4 = np.pi ** (-1. / 4)
+    # pim4=0.7511255444649425
+    pim4 = np.pi ** (-1. / 4)
 
     # The roots are symmetric about the origin, so we have to
     # find only half of them.
@@ -237,49 +236,42 @@ def _h_roots_newton(n, releps=3e-14, max_iter=10):
     theta = _get_theta(rhs)
     z = sqrt(anu) * np.cos(theta)
 
-    L = zeros((3, len(z)))
-    k0 = 0
-    kp1 = 1
+    p = zeros((3, len(z)))
+    k_0 = 0
+    k_p1 = 1
 
-    for i in range(max_iter):  # @UnusedVariable
+    for _i in range(max_iter):
         # Newtons method carried out simultaneously on the roots.
-        L[k0, :] = 0
-        L[kp1, :] = PIM4
+        p[k_0, :] = 0
+        p[k_p1, :] = pim4
 
         for j in range(1, n + 1):
             # Loop up the recurrence relation to get the Hermite
             # polynomials evaluated at z.
-            km1 = k0
-            k0 = kp1
-            kp1 = np.mod(kp1 + 1, 3)
+            k_m1 = k_0
+            k_0 = k_p1
+            k_p1 = np.mod(k_p1 + 1, 3)
 
-            L[kp1, :] = (z * sqrt(2 / j) * L[k0, :] -
-                         np.sqrt((j - 1) / j) * L[km1, :])
+            p[k_p1, :] = (z * sqrt(2 / j) * p[k_0, :] -
+                          sqrt((j - 1) / j) * p[k_m1, :])
 
-        # L now contains the desired Hermite polynomials.
-        # We next compute pp, the derivatives,
+        # p now contains the desired Hermite polynomials.
+        # We next compute p_deriv, the derivatives,
         # by the relation (4.5.21) using p2, the polynomials
         # of one lower order.
 
-        pp = sqrt(2 * n) * L[k0, :]
-        dz = L[kp1, :] / pp
+        p_deriv = sqrt(2 * n) * p[k_0, :]
+        d_z = p[k_p1, :] / p_deriv
 
-        z = z - dz  # Newtons formula.
+        z = z - d_z  # Newtons formula.
 
-        converged = not np.any(abs(dz) > releps)
+        converged = not np.any(abs(d_z) > releps)
         if converged:
             break
 
     _assert_warn(converged, 'Newton iteration did not converge!')
-    w = 2. / pp ** 2
-    return _expand_roots(z, w, n, m)
-#     x = np.empty(n)
-#     w = np.empty(n)
-#     x[0:m] = z      # Store the root
-#     x[n - 1:n - m - 1:-1] = -z     # and its symmetric counterpart.
-#     w[0:m] = 2. / pp ** 2    # Compute the weight
-#     w[n - 1:n - m - 1:-1] = w[0:m]  # and its symmetric counterpart.
-#     return x, w
+    weights = 2. / p_deriv ** 2
+    return _expand_roots(z, weights, n, m)
 
 
 def h_roots(n, method='newton'):
@@ -337,48 +329,50 @@ def _j_roots_newton(n, alpha, beta, releps=3e-14, max_iter=10):
     z = np.cos(np.pi * (np.arange(1, n + 1) - 0.25 + 0.5 * alpha) /
                (n + 0.5 * (alfbet + 1)))
 
-    L = zeros((3, len(z)))
-    k0 = 0
-    kp1 = 1
-    for i in range(max_iter):
+    p = zeros((3, len(z)))
+    k_0 = 0
+    k_p1 = 1
+    for _i in range(max_iter):
         # Newton's method carried out simultaneously on the roots.
         tmp = 2 + alfbet
-        L[k0, :] = 1
-        L[kp1, :] = (alpha - beta + tmp * z) / 2
+        p[k_0, :] = 1
+        p[k_p1, :] = (alpha - beta + tmp * z) / 2
 
         for j in range(2, n + 1):
             # Loop up the recurrence relation to get the Jacobi
             # polynomials evaluated at z.
-            km1 = k0
-            k0 = kp1
-            kp1 = np.mod(kp1 + 1, 3)
+            k_m1 = k_0
+            k_0 = k_p1
+            k_p1 = np.mod(k_p1 + 1, 3)
 
             a = 2. * j * (j + alfbet) * tmp
             tmp = tmp + 2
             c = 2 * (j - 1 + alpha) * (j - 1 + beta) * tmp
             b = (tmp - 1) * (alpha ** 2 - beta ** 2 + tmp * (tmp - 2) * z)
 
-            L[kp1, :] = (b * L[k0, :] - c * L[km1, :]) / a
+            p[k_p1, :] = (b * p[k_0, :] - c * p[k_m1, :]) / a
 
-        # L now contains the desired Jacobi polynomials.
-        # We next compute pp, the derivatives with a standard
+        # p now contains the desired Jacobi polynomials.
+        # We next compute p_deriv, the derivatives with a standard
         # relation involving the polynomials of one lower order.
 
-        pp = ((n * (alpha - beta - tmp * z) * L[kp1, :] +
-               2 * (n + alpha) * (n + beta) * L[k0, :]) / (tmp * (1 - z ** 2)))
-        dz = L[kp1, :] / pp
-        z = z - dz  # Newton's formula.
+        p_deriv = ((n * (alpha - beta - tmp * z) * p[k_p1, :] +
+                    2 * (n + alpha) * (n + beta) * p[k_0, :]) /
+                   (tmp * (1 - z ** 2)))
+        d_z = p[k_p1, :] / p_deriv
+        z = z - d_z  # Newton's formula.
 
-        if not any(abs(dz) > releps * abs(z)):
+        converged = not any(abs(d_z) > releps * abs(z))
+        if converged:
             break
-    else:
-        warnings.warn('too many iterations in jrule')
+
+    _assert_warn(converged, 'too many iterations in jrule')
 
     x = z  # Store the root and the weight.
     f = (sp.gammaln(alpha + n) + sp.gammaln(beta + n) -
          sp.gammaln(n + 1) - sp.gammaln(alpha + beta + n + 1))
-    w = (np.exp(f) * tmp * 2 ** alfbet / (pp * L[k0, :]))
-    return x, w
+    weights = (np.exp(f) * tmp * 2 ** alfbet / (p_deriv * p[k_0, :]))
+    return x, weights
 
 
 def j_roots(n, alpha, beta, method='newton'):
@@ -483,10 +477,10 @@ def la_roots(n, alpha=0, method='newton'):
 
 
 def _get_theta(rhs):
-    r3 = rhs ** (1. / 3)
-    r2 = r3 ** 2
-    C = [9.084064e-01, 5.214976e-02, 2.579930e-03, 3.986126e-03]
-    theta = r3 * (C[0] + r2 * (C[1] + r2 * (C[2] + r2 * C[3])))
+    r_3 = rhs ** (1. / 3)
+    r_2 = r_3 ** 2
+    c = [9.084064e-01, 5.214976e-02, 2.579930e-03, 3.986126e-03]
+    theta = r_3 * (c[0] + r_2 * (c[1] + r_2 * (c[2] + r_2 * c[3])))
     return theta
 
 
@@ -498,57 +492,57 @@ def _la_roots_newton(n, alpha, releps=3e-14, max_iter=10):
     theta = _get_theta(rhs)
     z = anu * np.cos(theta) ** 2
 
-    dz = zeros(len(z))
-    L = zeros((3, len(z)))
-    Lp = zeros((1, len(z)))
-    pp = zeros((1, len(z)))
-    k0 = 0
-    kp1 = 1
+    d_z = zeros(len(z))
+    p = zeros((3, len(z)))
+    p_previous = zeros((1, len(z)))
+    p_deriv = zeros((1, len(z)))
+    k_0 = 0
+    k_p1 = 1
     k = slice(len(z))
-    for _its in range(max_iter):
+    for _i in range(max_iter):
         # Newton's method carried out simultaneously on the roots.
-        L[k0, k] = 0.
-        L[kp1, k] = 1.
+        p[k_0, k] = 0.
+        p[k_p1, k] = 1.
 
-        for jj in range(1, n + 1):
+        for j in range(1, n + 1):
             # Loop up the recurrence relation to get the Laguerre
             # polynomials evaluated at z.
-            km1 = k0
-            k0 = kp1
-            kp1 = np.mod(kp1 + 1, 3)
+            km1 = k_0
+            k_0 = k_p1
+            k_p1 = np.mod(k_p1 + 1, 3)
 
-            L[kp1, k] = ((2 * jj - 1 + alpha - z[k]) * L[k0, k] -
-                         (jj - 1 + alpha) * L[km1, k]) / jj
+            p[k_p1, k] = ((2 * j - 1 + alpha - z[k]) * p[k_0, k] -
+                          (j - 1 + alpha) * p[km1, k]) / j
         # end
-        # L now contains the desired Laguerre polynomials.
-        # We next compute pp, the derivatives with a standard
-        #  relation involving the polynomials of one lower order.
+        # p now contains the desired Laguerre polynomials.
+        # We next compute p_deriv, the derivatives with a standard
+        # relation involving the polynomials of one lower order.
 
-        Lp[k] = L[k0, k]
-        pp[k] = (n * L[kp1, k] - (n + alpha) * Lp[k]) / z[k]
+        p_previous[k] = p[k_0, k]
+        p_deriv[k] = (n * p[k_p1, k] - (n + alpha) * p_previous[k]) / z[k]
 
-        dz[k] = L[kp1, k] / pp[k]
-        z[k] = z[k] - dz[k]  # % Newton?s formula.
-        # k = find((abs(dz) > releps.*z))
+        d_z[k] = p[k_p1, k] / p_deriv[k]
+        z[k] = z[k] - d_z[k]  # Newton?s formula.
+        # k = find((abs(d_z) > releps.*z))
 
-        converged = not np.any(abs(dz) > releps)
+        converged = not np.any(abs(d_z) > releps)
         if converged:
             break
 
-    _assert_warn('too many iterations!')
+    _assert_warn(converged, 'too many iterations!')
 
-    x = z
-    w = -np.exp(sp.gammaln(alpha + n) - sp.gammaln(n)) / (pp * n * Lp)
-    return x, w
+    nodes = z
+    weights = -np.exp(sp.gammaln(alpha + n) -
+                      sp.gammaln(n)) / (p_deriv * n * p_previous)
+    return nodes, weights
 
 
 def _p_roots_newton_start(n):
     m = int(np.fix((n + 1) / 2))
-    mm = 4 * m - 1
-    t = (np.pi / (4 * n + 2)) * np.arange(3, mm + 1, 4)
-    nn = 1 - (1 - 1 / n) / (8 * n * n)
-    xo = nn * np.cos(t)
-    return m, xo
+    t = (np.pi / (4 * n + 2)) * np.arange(3, 4 * m, 4)
+    a = 1 - (1 - 1 / n) / (8 * n * n)
+    x = a * np.cos(t)
+    return m, x
 
 
 def _p_roots_newton(n, ):
@@ -556,52 +550,52 @@ def _p_roots_newton(n, ):
     Algorithm given by Davis and Rabinowitz in 'Methods
     of Numerical Integration', page 365, Academic Press, 1975.
     """
-    m, xo = _p_roots_newton_start(n)
+    m, x = _p_roots_newton_start(n)
 
-    e1 = n * (n + 1)
+    e_1 = n * (n + 1)
     for _j in range(2):
-        pkm1 = 1
-        pk = xo
+        p_km1 = 1
+        p_k = x
         for k in range(2, n + 1):
-            t1 = xo * pk
-            pkp1 = t1 - pkm1 - (t1 - pkm1) / k + t1
-            pkm1 = pk
-            pk = pkp1
-
-        den = 1. - xo * xo
-        d1 = n * (pkm1 - xo * pk)
-        dpn = d1 / den
-        d2pn = (2. * xo * dpn - e1 * pk) / den
-        d3pn = (4. * xo * d2pn + (2 - e1) * dpn) / den
-        d4pn = (6. * xo * d3pn + (6 - e1) * d2pn) / den
-        u = pk / dpn
-        v = d2pn / dpn
-        h = -u * (1 + (.5 * u) * (v + u * (v * v - u * d3pn / (3 * dpn))))
-        p = pk + h * (dpn + (.5 * h) * (d2pn +
-                                        (h / 3) * (d3pn + .25 * h * d4pn)))
-        dp = dpn + h * (d2pn + (.5 * h) * (d3pn + h * d4pn / 3))
-        h = h - p / dp
-        xo = xo + h
-
-    x = -xo - h
-    fx = d1 - h * e1 * (pk + (h / 2) * (dpn + (h / 3) *
-                                        (d2pn + (h / 4) *
-                                         (d3pn + (.2 * h) * d4pn))))
-    w = 2 * (1 - x ** 2) / (fx ** 2)
-    return _expand_roots(x, w, n, m)
+            t_1 = x * p_k
+            p_kp1 = t_1 - p_km1 - (t_1 - p_km1) / k + t_1
+            p_km1 = p_k
+            p_k = p_kp1
+
+        den = 1. - x * x
+        d_1 = n * (p_km1 - x * p_k)
+        d_pn = d_1 / den
+        d_2pn = (2. * x * d_pn - e_1 * p_k) / den
+        d_3pn = (4. * x * d_2pn + (2 - e_1) * d_pn) / den
+        d_4pn = (6. * x * d_3pn + (6 - e_1) * d_2pn) / den
+        u = p_k / d_pn
+        v = d_2pn / d_pn
+        h = -u * (1 + (.5 * u) * (v + u * (v * v - u * d_3pn / (3 * d_pn))))
+        p = p_k + h * (d_pn + (.5 * h) * (d_2pn + (h / 3) *
+                                          (d_3pn + .25 * h * d_4pn)))
+        d_p = d_pn + h * (d_2pn + (.5 * h) * (d_3pn + h * d_4pn / 3))
+        h = h - p / d_p
+        x = x + h
+
+    nodes = -x - h
+    f_x = d_1 - h * e_1 * (p_k + (h / 2) * (d_pn + (h / 3) *
+                                            (d_2pn + (h / 4) *
+                                             (d_3pn + (.2 * h) * d_4pn))))
+    weights = 2 * (1 - nodes ** 2) / (f_x ** 2)
+    return _expand_roots(nodes, weights, n, m)
 
 
 def _p_roots_newton1(n, releps=1e-15, max_iter=100):
-    m, xo = _p_roots_newton_start(n)
+    m, x = _p_roots_newton_start(n)
     # Compute the zeros of the N+1 Legendre Polynomial
     # using the recursion relation and the Newton-Raphson method
 
     # Legendre-Gauss Polynomials
-    L = zeros((3, m))
+    p = zeros((3, m))
 
     # Derivative of LGP
-    Lp = zeros((m,))
-    dx = zeros((m,))
+    p_deriv = zeros((m,))
+    d_x = zeros((m,))
 
     # Compute the zeros of the N+1 Legendre Polynomial
     # using the recursion relation and the Newton-Raphson method
@@ -609,32 +603,33 @@ def _p_roots_newton1(n, releps=1e-15, max_iter=100):
     # Iterate until new points are uniformly within epsilon of old
     # points
     k = slice(m)
-    k0 = 0
-    kp1 = 1
+    k_0 = 0
+    k_p1 = 1
     for _ix in range(max_iter):
-        L[k0, k] = 1
-        L[kp1, k] = xo[k]
-
-        for jj in range(2, n + 1):
-            km1 = k0
-            k0 = kp1
-            kp1 = np.mod(k0 + 1, 3)
-            L[kp1, k] = ((2 * jj - 1) * xo[k] * L[k0, k] -
-                         (jj - 1) * L[km1, k]) / jj
-
-        Lp[k] = n * (L[k0, k] - xo[k] * L[kp1, k]) / (1 - xo[k] ** 2)
-
-        dx[k] = L[kp1, k] / Lp[k]
-        xo[k] = xo[k] - dx[k]
-        k, = np.nonzero((abs(dx) > releps * np.abs(xo)))
-        if len(k) == 0:
+        p[k_0, k] = 1
+        p[k_p1, k] = x[k]
+
+        for j in range(2, n + 1):
+            k_m1 = k_0
+            k_0 = k_p1
+            k_p1 = np.mod(k_0 + 1, 3)
+            p[k_p1, k] = ((2 * j - 1) * x[k] * p[k_0, k] -
+                          (j - 1) * p[k_m1, k]) / j
+
+        p_deriv[k] = n * (p[k_0, k] - x[k] * p[k_p1, k]) / (1 - x[k] ** 2)
+
+        d_x[k] = p[k_p1, k] / p_deriv[k]
+        x[k] = x[k] - d_x[k]
+        k, = np.nonzero((abs(d_x) > releps * np.abs(x)))
+        converged = len(k) == 0
+        if converged:
             break
-    else:
-        warnings.warn('Too many iterations!')
 
-    x = -xo
-    w = 2. / ((1 - x ** 2) * (Lp ** 2))
-    return _expand_roots(x, w, n, m)
+    _assert(converged, 'Too many iterations!')
+
+    nodes = -x
+    weights = 2. / ((1 - nodes ** 2) * (p_deriv ** 2))
+    return _expand_roots(nodes, weights, n, m)
 
 
 def _expand_roots(x, w, n, m):
@@ -650,7 +645,7 @@ def _expand_roots(x, w, n, m):
 def p_roots(n, method='newton', a=-1, b=1):
     """
     Returns the roots (x) of the nth order Legendre polynomial, P_n(x),
-    and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting
+    and weights to use in Gaussian Quadrature over [-1,1] with weighting
     function 1.
 
     Parameters
@@ -664,9 +659,9 @@ def p_roots(n, method='newton', a=-1, b=1):
 
     Returns
     -------
-    x : ndarray
+    nodes : ndarray
         roots
-    w : ndarray
+    weights : ndarray
         weights
 
 
@@ -674,8 +669,8 @@ def p_roots(n, method='newton', a=-1, b=1):
     -------
     Integral of exp(x) from a = 0 to b = 3 is: exp(3)-exp(0)=
     >>> import numpy as np
-    >>> x, w = p_roots(11, a=0, b=3)
-    >>> np.allclose(np.sum(np.exp(x)*w), 19.085536923187668)
+    >>> nodes, weights = p_roots(11, a=0, b=3)
+    >>> np.allclose(np.sum(np.exp(nodes) * weights), 19.085536923187668)
     True
 
     See also
@@ -697,74 +692,74 @@ def p_roots(n, method='newton', a=-1, b=1):
     """
 
     if not method.startswith('n'):
-        x, w = ort.p_roots(n)
+        nodes, weights = ort.p_roots(n)
     else:
         if method.endswith('1'):
-            x, w = _p_roots_newton1(n)
+            nodes, weights = _p_roots_newton1(n)
         else:
-            x, w = _p_roots_newton(n)
+            nodes, weights = _p_roots_newton(n)
 
     if (a != -1) | (b != 1):
         # Linear map from[-1,1] to [a,b]
-        dh = (b - a) / 2
-        x = dh * (x + 1) + a
-        w = w * dh
+        d_h = (b - a) / 2
+        nodes = d_h * (nodes + 1) + a
+        weights = weights * d_h
 
-    return x, w
+    return nodes, weights
 
 
 def q5_roots(n):
     """
     5      : p(x) = 1/sqrt((x-a)*(b-x)), a =-1,   b = 1 Chebyshev 1'st kind
     """
-    jj = np.arange(1, n + 1)
-    wf = ones(n) * np.pi / n
-    bp = np.cos((2 * jj - 1) * np.pi / (2 * n))
-    return bp, wf
+    j = np.arange(1, n + 1)
+    weights = ones(n) * np.pi / n
+    nodes = np.cos((2 * j - 1) * np.pi / (2 * n))
+    return nodes, weights
 
 
 def q6_roots(n):
     """
     6      : p(x) = sqrt((x-a)*(b-x)),   a =-1,   b = 1 Chebyshev 2'nd kind
     """
-    jj = np.arange(1, n + 1)
-    xj = jj * np.pi / (n + 1)
-    wf = np.pi / (n + 1) * np.sin(xj) ** 2
-    bp = np.cos(xj)
-    return bp, wf
+    j = np.arange(1, n + 1)
+    x_j = j * np.pi / (n + 1)
+    weights = np.pi / (n + 1) * np.sin(x_j) ** 2
+    nodes = np.cos(x_j)
+    return nodes, weights
 
 
 def q7_roots(n):
     """
     7 : p(x) = sqrt((x-a)/(b-x)),   a = 0,   b = 1
     """
-    jj = np.arange(1, n + 1)
-    xj = (jj - 0.5) * pi / (2 * n + 1)
-    bp = np.cos(xj) ** 2
-    wf = 2 * np.pi * bp / (2 * n + 1)
-    return bp, wf
+    j = np.arange(1, n + 1)
+    x_j = (j - 0.5) * pi / (2 * n + 1)
+    nodes = np.cos(x_j) ** 2
+    weights = 2 * np.pi * nodes / (2 * n + 1)
+    return nodes, weights
 
 
 def q8_roots(n):
     """
     8 : p(x) = 1/sqrt(b-x),         a = 0,   b = 1
     """
-    bp1, wf1 = p_roots(2 * n)
-    k, = np.where(0 <= bp1)
-    wf = 2 * wf1[k]
-    bp = 1 - bp1[k] ** 2
-    return bp, wf
+    nodes_1, weights_1 = p_roots(2 * n)
+    k, = np.where(0 <= nodes_1)
+    weights = 2 * weights_1[k]
+    nodes = 1 - nodes_1[k] ** 2
+    return nodes, weights
 
 
 def q9_roots(n):
     """
      9 : p(x) = sqrt(b-x),           a = 0,   b = 1
     """
-    bp1, wf1 = p_roots(2 * n + 1)
-    k, = np.where(0 < bp1)
-    wf = 2 * bp1[k] ** 2 * wf1[k]
-    bp = 1 - bp1[k] ** 2
-    return bp, wf
+    nodes_1, weights_1 = p_roots(2 * n + 1)
+    k, = np.where(0 < nodes_1)
+    weights = 2 * nodes_1[k] ** 2 * weights_1[k]
+    nodes = 1 - nodes_1[k] ** 2
+    return nodes, weights
 
 
 def qrule(n, wfun=1, alpha=0, beta=0):
@@ -805,18 +800,18 @@ def qrule(n, wfun=1, alpha=0, beta=0):
     >>> import numpy as np
 
     # integral of x^2 from a = -1 to b = 1
-    >>> [bp,wf] = qrule(10)
+    >>> bp, wf = qrule(10)
     >>> np.allclose(sum(bp**2*wf), 0.66666666666666641)
     True
 
-    # integral of exp(-x.^2)*x.^2 from a = -inf to b = inf
-    >>> [bp,wf] = qrule(10,2)
-    >>> np.allclose(sum(bp**2*wf), 0.88622692545275772)
+    # integral of exp(-x**2)*x**2 from a = -inf to b = inf
+    >>> bp, wf = qrule(10,2)
+    >>> np.allclose(sum(bp ** 2 * wf), 0.88622692545275772)
     True
 
-    # integral of (x+1)*(1-x)^2 from  a = -1 to b = 1
-    >>> [bp,wf] = qrule(10,4,1,2)
-    >>> np.allclose((bp*wf).sum(), 0.26666666666666755)
+    # integral of (x+1)*(1-x)**2 from  a = -1 to b = 1
+    >>> bp, wf = qrule(10,4,1,2)
+    >>> np.allclose((bp * wf).sum(), 0.26666666666666755)
     True
 
     See also
@@ -928,22 +923,22 @@ class _Gaussq(object):
         def fun1(x):
             return x
         if wfun == 4:
-            dx = jacob ** (alpha + beta + 1)
+            d_x = jacob ** (alpha + beta + 1)
         else:
-            dx = [None, fun1, fun1, fun1, None, lambda x: ones(np.shape(x)),
-                  lambda x: x ** 2, fun1, lambda x: sqrt(x),
-                  lambda x: sqrt(x) ** 3][wfun](jacob)
-        return dx.ravel()
+            d_x = [None, fun1, fun1, fun1, None, lambda x: ones(np.shape(x)),
+                   lambda x: x ** 2, fun1, sqrt,
+                   lambda x: sqrt(x) ** 3][wfun](jacob)
+        return d_x.ravel()
 
     @staticmethod
-    def _points_and_weights(gn, wfun, alpha, beta):
-        global _POINTS_AND_WEIGHTS
-        name = 'wfun{:d}_{:d}_{:g}_{:g}'.format(wfun, gn, alpha, beta)
-        x_and_w = _POINTS_AND_WEIGHTS.setdefault(name, [])
-        if len(x_and_w) == 0:
-            x_and_w.extend(qrule(gn, wfun, alpha, beta))
-        xn, w = x_and_w
-        return xn, w
+    def _nodes_and_weights(num_nodes, wfun, alpha, beta):
+        global _NODES_AND_WEIGHTS
+        name = 'wfun{:d}_{:d}_{:g}_{:g}'.format(wfun, num_nodes, alpha, beta)
+        nodes_and_weights = _NODES_AND_WEIGHTS.setdefault(name, [])
+        if len(nodes_and_weights) == 0:
+            nodes_and_weights.extend(qrule(num_nodes, wfun, alpha, beta))
+        nodes, weights = nodes_and_weights
+        return nodes, weights
 
     def _initialize_trace(self, max_iter):
         if self.trace:
@@ -963,24 +958,23 @@ class _Gaussq(object):
             plt.plot(np.hstack(self.x_trace), np.hstack(self.y_trace), '+')
 
     @staticmethod
-    def _get_jacob(wfun, A, B):
+    def _get_jacob(wfun, a, b):
         if wfun in [2, 3]:
-            nk = np.size(A)
-            jacob = ones((nk, 1))
+            jacob = ones((np.size(a), 1))
         else:
-            jacob = (B - A) * 0.5
+            jacob = (b - a) * 0.5
             if wfun in [7, 8, 9]:
-                jacob = jacob * 2
+                jacob *= 2
         return jacob
 
     @staticmethod
     def _warn(k, a_shape):
-        nk = len(k)
-        if nk > 1:
-            if nk == np.prod(a_shape):
+        n = len(k)
+        if n > 1:
+            if n == np.prod(a_shape):
                 tmptxt = 'All integrals did not converge'
             else:
-                tmptxt = '%d integrals did not converge' % (nk, )
+                tmptxt = '%d integrals did not converge' % (n, )
             tmptxt = tmptxt + '--singularities likely!'
         else:
             tmptxt = 'Integral did not converge--singularity likely!'
@@ -1000,42 +994,43 @@ class _Gaussq(object):
     def __call__(self, fun, a, b, releps=1e-3, abseps=1e-3, alpha=0, beta=0,
                  wfun=1, trace=False, args=(), max_iter=11):
         self.trace = trace
-        gn = 2
+        num_nodes = 2
 
-        aa, bb, args, a_shape = self._initialize(wfun, a, b, args)
+        a_0, b_0, args, a_shape = self._initialize(wfun, a, b, args)
 
-        jacob = self._get_jacob(wfun, aa, bb)
+        jacob = self._get_jacob(wfun, a_0, b_0)
         shift = int(wfun in [1, 4, 5, 6])
-        dx = self._get_dx(wfun, jacob, alpha, beta)
+        d_x = self._get_dx(wfun, jacob, alpha, beta)
 
         self._initialize_trace(max_iter)
 
         # Break out of the iteration loop for three reasons:
         #  1) the last update is very small (compared to int and to releps)
         #  2) There are more than 11 iterations. This should NEVER happen.
-        dtype = np.result_type(fun((aa+bb)*0.5, *args))
-        nk = np.prod(a_shape)  # # of integrals we have to compute
-        k = np.arange(nk)
-        opts = (nk, dtype)
+        dtype = np.result_type(fun((a_0+b_0)*0.5, *args))
+        n_k = np.prod(a_shape)  # # of integrals we have to compute
+        k = np.arange(n_k)
+        opts = (n_k, dtype)
         val, val_old, abserr = zeros(*opts), ones(*opts), zeros(*opts)
+        nodes_and_weights = self._nodes_and_weights
         for i in range(max_iter):
-            xn, w = self._points_and_weights(gn, wfun, alpha, beta)
-            x = (xn + shift) * jacob[k, :] + aa[k, :]
+            x_n, weights = nodes_and_weights(num_nodes, wfun, alpha, beta)
+            x = (x_n + shift) * jacob[k, :] + a_0[k, :]
 
-            pari = [xi[k, :] for xi in args]
-            y = fun(x, *pari)
+            params = [xi[k, :] for xi in args]
+            y = fun(x, *params)
             self._plot_trace(x, y)
-            val[k] = np.sum(w * y, axis=1) * dx[k]  # do the integration
+            val[k] = np.sum(weights * y, axis=1) * d_x[k]  # do the integration
             if any(np.isnan(val)):
                 val[np.isnan(val)] = val_old[np.isnan(val)]
             if 1 < i:
                 abserr[k] = abs(val_old[k] - val[k])  # absolute tolerance
                 k, = np.where(abserr > np.maximum(abs(releps * val), abseps))
-                nk = len(k)  # of integrals we have to compute again
-                if nk == 0:
+                n_k = len(k)  # of integrals we have to compute again
+                if n_k == 0:
                     break
             val_old[k] = val[k]
-            gn *= 2  # double the # of basepoints and weights
+            num_nodes *= 2  # double the # of basepoints and weights
         else:
             self._warn(k, a_shape)
 
@@ -1048,14 +1043,13 @@ class _Gaussq(object):
 gaussq = _Gaussq()
 
 
-def richardson(Q, k):
+def richardson(q_val, k):
     # license BSD
     # Richardson extrapolation with parameter estimation
-    c = np.real((Q[k - 1] - Q[k - 2]) / (Q[k] - Q[k - 1])) - 1.
+    c = np.real((q_val[k - 1] - q_val[k - 2]) / (q_val[k] - q_val[k - 1])) - 1.
     # The lower bound 0.07 admits the singularity x.^-0.9
     c = max(c, 0.07)
-    R = Q[k] + (Q[k] - Q[k - 1]) / c
-    return R
+    return q_val[k] + (q_val[k] - q_val[k - 1]) / c
 
 
 class _Quadgr(object):
@@ -1083,15 +1077,15 @@ class _Quadgr(object):
     @staticmethod
     def _nodes_and_weights():
         # Gauss-Legendre quadrature (12-point)
-        xq = np.asarray(
+        x = np.asarray(
             [0.12523340851146894, 0.36783149899818018, 0.58731795428661748,
              0.76990267419430469, 0.9041172563704748, 0.98156063424671924])
-        wq = np.asarray(
+        w = np.asarray(
             [0.24914704581340288, 0.23349253653835478, 0.20316742672306584,
              0.16007832854334636, 0.10693932599531818, 0.047175336386511842])
-        xq = np.hstack((xq, -xq))
-        wq = np.hstack((wq, wq))
-        return xq, wq
+        nodes = np.hstack((x, -x))
+        weights = np.hstack((w, w))
+        return nodes, weights
 
     @staticmethod
     def _get_best_estimate(vals0, vals1, vals2, k):
@@ -1121,24 +1115,23 @@ class _Quadgr(object):
         val1 = zeros(max_iter, dtype=dtype)  # First Richardson extrapolation
         val2 = zeros(max_iter, dtype=dtype)  # Second Richardson extrapolation
 
-        x0, w0 = self._nodes_and_weights()
-        nx0 = len(x0)
+        x_n, weights = self._nodes_and_weights()
+        n = len(x_n)
         # One interval
-        hh = (b - a) / 2             # Half interval length
-        x = (a + b) / 2 + hh * x0      # Nodes
+        d_x = (b - a) / 2                # Half interval length
+        x = (a + b) / 2 + d_x * x_n      # Nodes
         # Quadrature
-        val0[0] = hh * np.sum(w0 * fun(x), axis=0)
+        val0[0] = d_x * np.sum(weights * fun(x), axis=0)
 
         # Successive bisection of intervals
         for k in range(1, max_iter):
 
             # Interval bisection
-            hh = hh / 2
+            d_x = d_x / 2
             x = np.hstack([x + a, x + b]) / 2
             # Quadrature
-            val0[k] = hh * np.sum(w0 * np.sum(np.reshape(fun(x), (-1, nx0)),
-                                              axis=0),
-                                  axis=0)
+            val0[k] = np.sum(np.sum(np.reshape(fun(x), (-1, n)), axis=0) *
+                             weights, axis=0) * d_x
 
             # Richardson extrapolation
             if k >= 5:
@@ -1147,22 +1140,20 @@ class _Quadgr(object):
             elif k >= 3:
                 val1[k] = richardson(val0, k)
 
-            Q, err = self._get_best_estimate(val0, val1, val2, k)
-
-            # Convergence
+            q_val, err = self._get_best_estimate(val0, val1, val2, k)
 
-            converged = (err < abseps) | ~np.isfinite(Q)
+            converged = (err < abseps) | ~np.isfinite(q_val)
             if converged:
                 break
 
         _assert_warn(converged, 'Max number of iterations reached without '
                      'convergence.')
-        _assert_warn(np.isfinite(Q),
+        _assert_warn(np.isfinite(q_val),
                      'Integral approximation is Infinite or NaN.')
 
         # The error estimate should not be zero
-        err = err + 2 * np.finfo(Q).eps
-        return Q, err
+        err = err + 2 * np.finfo(q_val).eps
+        return q_val, err
 
     @staticmethod
     def _order_limits(a, b):
@@ -1174,9 +1165,9 @@ class _Quadgr(object):
         """
         Gauss-Legendre quadrature with Richardson extrapolation.
 
-        [Q,ERR] = QUADGR(FUN,A,B,TOL) approximates the integral of a function
+        [q_val,ERR] = QUADGR(FUN,A,B,TOL) approximates the integral of a function
         FUN from A to B with an absolute error tolerance TOL. FUN is a function
-        handle and must accept vector arguments. TOL is 1e-6 by default. Q is
+        handle and must accept vector arguments. TOL is 1e-6 by default. q_val is
         the integral approximation and ERR is an estimate of the absolute
         error.
 
@@ -1188,7 +1179,7 @@ class _Quadgr(object):
         Examples
         --------
         >>> import numpy as np
-        >>> Q, err = quadgr(np.log,0,1)
+        >>> q_val, err = quadgr(np.log,0,1)
         >>> q, err = quadgr(np.exp,0,9999*1j*np.pi)
         >>> np.allclose(q, -2.0000000000122662), err < 1.0e-08
         (True, True)
@@ -1223,25 +1214,24 @@ class _Quadgr(object):
         a = np.asarray(a)
         b = np.asarray(b)
         if a == b:
-            Q = b - a
+            q_val = b - a
             err = b - a
-            return Q, err
+            return q_val, err
 
         a, b, reverse = self._order_limits(a, b)
 
-        # Infinite limits
         improper_integral = np.isinf(a) | np.isinf(b)
-        if improper_integral:
-            Q, err = self._change_variable_and_integrate(fun, a, b, abseps,
-                                                         max_iter)
+        if improper_integral:  # Infinite limits
+            q_val, err = self._change_variable_and_integrate(fun, a, b, abseps,
+                                                             max_iter)
         else:
-            Q, err = self._integrate(fun, a, b, abseps, max_iter)
+            q_val, err = self._integrate(fun, a, b, abseps, max_iter)
 
         # Reverse direction
         if reverse:
-            Q = -Q
+            q_val = -q_val
 
-        return Q, err
+        return q_val, err
 
 quadgr = _Quadgr()
 
@@ -1380,13 +1370,13 @@ def qdemo(f, a, b, kmax=9, plot_error=False):
         err_dic.setdefault(name, []).append(abs(q[0] - true_val))
 
         name = 'Chebychev'
-        ck = np.polynomial.chebyshev.chebfit(x, y, deg=min(n-1, 36))
-        cki = np.polynomial.chebyshev.chebint(ck)
-        q = np.polynomial.chebyshev.chebval(x[-1], cki)
+        c_k = np.polynomial.chebyshev.chebfit(x, y, deg=min(n-1, 36))
+        c_ki = np.polynomial.chebyshev.chebint(c_k)
+        q = np.polynomial.chebyshev.chebval(x[-1], c_ki)
         vals_dic.setdefault(name, []).append(q)
         err_dic.setdefault(name, []).append(abs(q - true_val))
-        # ck = chebfit(f,n,a,b)
-        # q  = chebval(b,chebint(ck,a,b),a,b)
+        # c_k = chebfit(f,n,a,b)
+        # q  = chebval(b,chebint(c_k,a,b),a,b)
         # qc2[k] = q; ec2[k] = abs(q - true)
 
         name = 'Gauss-Legendre'  # quadrature