diff --git a/doc/Makefile b/doc/Makefile
index cae96bdb67392c43ff359b262b3c218ece5947ef..cd9f910938e3358daaf585921710b3b3c4081df8 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -3,7 +3,7 @@
 
 # You can set these variables from the command line.
 SPHINXOPTS    =
-SPHINXBUILD   = sphinx-build
+SPHINXBUILD   = python `which sphinx-build`
 PAPER         =
 BUILDDIR      = _build
 
diff --git a/requirements.txt b/requirements.txt
index 6e23ae957e1531c3a3dba6762f93010b69f94b07..6690e3d9651bad06a8fd28e444167fe03bbb119f 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,7 +2,7 @@ numpy
 sympy==1.0
 git+https://github.com/inducer/pymbolic
 git+https://github.com/inducer/islpy
-git+https://github.com/pyopencl/pyopencl
+git+https://github.com/inducer/pyopencl
 git+https://github.com/inducer/boxtree
 git+https://github.com/inducer/loopy
 git+https://github.com/inducer/pyfmmlib
diff --git a/sumpy/__init__.py b/sumpy/__init__.py
index b7d4293e5ac65e2af04a91552105f30455615bd6..c770e802e8bab11241b3f9593e5227ab7cb7d8ba 100644
--- a/sumpy/__init__.py
+++ b/sumpy/__init__.py
@@ -1,5 +1,4 @@
-from __future__ import division
-from __future__ import absolute_import
+from __future__ import division, absolute_import
 
 __copyright__ = "Copyright (C) 2013 Andreas Kloeckner"
 
diff --git a/sumpy/codegen.py b/sumpy/codegen.py
index ca6086907165deb6b72b48f55b9ccc8350055777..5d27708f220fcf2b93b804da44ef6eeb1ef36565 100644
--- a/sumpy/codegen.py
+++ b/sumpy/codegen.py
@@ -477,7 +477,7 @@ class PowerRewriter(CSECachingMapperMixin, IdentityMapper):
             if exp > 1 and exp % 2 == 0:
                 square = prim.wrap_in_cse(new_base*new_base)
                 return self.rec(prim.wrap_in_cse(square**(exp//2)))
-            if exp > 1 and exp % 2 == 1:
+            elif exp > 1 and exp % 2 == 1:
                 square = prim.wrap_in_cse(new_base*new_base)
                 return self.rec(prim.wrap_in_cse(square**((exp-1)//2))*new_base)
             elif exp == 1:
diff --git a/sumpy/e2e.py b/sumpy/e2e.py
index bff0f508c05ea378fb5c7cd31ce85f82c0ffdfcd..eca056a7695f8e2d0b1b7c18931a5866cb6b3395 100644
--- a/sumpy/e2e.py
+++ b/sumpy/e2e.py
@@ -98,6 +98,9 @@ class E2EBase(KernelCacheWrapper):
 
         src_coeff_exprs = [sym.Symbol("src_coeff%d" % i)
                 for i in range(len(self.src_expansion))]
+        src_rscale = sym.Symbol("src_rscale")
+
+        tgt_rscale = sym.Symbol("tgt_rscale")
 
         from sumpy.assignment_collection import SymbolicAssignmentCollection
         sac = SymbolicAssignmentCollection()
@@ -105,7 +108,8 @@ class E2EBase(KernelCacheWrapper):
                 sac.assign_unique("coeff%d" % i, coeff_i)
                 for i, coeff_i in enumerate(
                     self.tgt_expansion.translate_from(
-                        self.src_expansion, src_coeff_exprs, dvec))]
+                        self.src_expansion, src_coeff_exprs, src_rscale,
+                        dvec, tgt_rscale))]
 
         sac.run_global_cse()
 
@@ -196,6 +200,7 @@ class E2EFromCSR(E2EBase):
                 """],
                 [
                     lp.GlobalArg("centers", None, shape="dim, aligned_nboxes"),
+                    lp.ValueArg("src_rscale,tgt_rscale", None),
                     lp.GlobalArg("src_box_starts, src_box_lists",
                         None, shape=None, strides=(1,), offset=lp.auto),
                     lp.ValueArg("aligned_nboxes,tgt_base_ibox,src_base_ibox",
@@ -228,11 +233,22 @@ class E2EFromCSR(E2EBase):
         :arg src_expansions:
         :arg src_box_starts:
         :arg src_box_lists:
+        :arg src_rscale:
+        :arg tgt_rscale:
         :arg centers:
         """
         knl = self.get_cached_optimized_kernel()
 
-        return knl(queue, **kwargs)
+        centers = kwargs.pop("centers")
+        # "1" may be passed for rscale, which won't have its type
+        # meaningfully inferred. Make the type of rscale explicit.
+        src_rscale = centers.dtype.type(kwargs.pop("src_rscale"))
+        tgt_rscale = centers.dtype.type(kwargs.pop("tgt_rscale"))
+
+        return knl(queue,
+                centers=centers,
+                src_rscale=src_rscale, tgt_rscale=tgt_rscale,
+                **kwargs)
 
 # }}}
 
@@ -306,6 +322,7 @@ class E2EFromChildren(E2EBase):
                     lp.GlobalArg("target_boxes", None, shape=lp.auto,
                         offset=lp.auto),
                     lp.GlobalArg("centers", None, shape="dim, aligned_nboxes"),
+                    lp.ValueArg("src_rscale,tgt_rscale", None),
                     lp.GlobalArg("box_child_ids", None,
                         shape="nchildren, aligned_nboxes"),
                     lp.GlobalArg("tgt_expansions", None,
@@ -337,11 +354,22 @@ class E2EFromChildren(E2EBase):
         :arg src_expansions:
         :arg src_box_starts:
         :arg src_box_lists:
+        :arg src_rscale:
+        :arg tgt_rscale:
         :arg centers:
         """
         knl = self.get_cached_optimized_kernel()
 
-        return knl(queue, **kwargs)
+        centers = kwargs.pop("centers")
+        # "1" may be passed for rscale, which won't have its type
+        # meaningfully inferred. Make the type of rscale explicit.
+        src_rscale = centers.dtype.type(kwargs.pop("src_rscale"))
+        tgt_rscale = centers.dtype.type(kwargs.pop("tgt_rscale"))
+
+        return knl(queue,
+                centers=centers,
+                src_rscale=src_rscale, tgt_rscale=tgt_rscale,
+                **kwargs)
 
 # }}}
 
@@ -402,6 +430,7 @@ class E2EFromParent(E2EBase):
                     lp.GlobalArg("target_boxes", None, shape=lp.auto,
                         offset=lp.auto),
                     lp.GlobalArg("centers", None, shape="dim, naligned_boxes"),
+                    lp.ValueArg("src_rscale,tgt_rscale", None),
                     lp.ValueArg("naligned_boxes,nboxes", np.int32),
                     lp.ValueArg("tgt_base_ibox,src_base_ibox", np.int32),
                     lp.ValueArg("ntgt_level_boxes,nsrc_level_boxes", np.int32),
@@ -430,11 +459,22 @@ class E2EFromParent(E2EBase):
         :arg src_expansions:
         :arg src_box_starts:
         :arg src_box_lists:
+        :arg src_rscale:
+        :arg tgt_rscale:
         :arg centers:
         """
         knl = self.get_cached_optimized_kernel()
 
-        return knl(queue, **kwargs)
+        centers = kwargs.pop("centers")
+        # "1" may be passed for rscale, which won't have its type
+        # meaningfully inferred. Make the type of rscale explicit.
+        src_rscale = centers.dtype.type(kwargs.pop("src_rscale"))
+        tgt_rscale = centers.dtype.type(kwargs.pop("tgt_rscale"))
+
+        return knl(queue,
+                centers=centers,
+                src_rscale=src_rscale, tgt_rscale=tgt_rscale,
+                **kwargs)
 
 # }}}
 
diff --git a/sumpy/e2p.py b/sumpy/e2p.py
index 9a2a8bc1ea564e37cf4e2d248de6952cfd386067..5a5c11549c98e89229801ebce9e83619c25a2176 100644
--- a/sumpy/e2p.py
+++ b/sumpy/e2p.py
@@ -81,12 +81,15 @@ class E2PBase(KernelCacheWrapper):
         from sumpy.symbolic import make_sym_vector
         bvec = make_sym_vector("b", self.dim)
 
+        import sumpy.symbolic as sp
+        rscale = sp.Symbol("rscale")
+
         from sumpy.assignment_collection import SymbolicAssignmentCollection
         sac = SymbolicAssignmentCollection()
 
         coeff_exprs = [sym.Symbol("coeff%d" % i)
                 for i in range(len(self.expansion.get_coefficient_identifiers()))]
-        value = self.expansion.evaluate(coeff_exprs, bvec)
+        value = self.expansion.evaluate(coeff_exprs, bvec, rscale)
 
         result_names = [
             sac.assign_unique("result_%d_p" % i,
@@ -112,7 +115,8 @@ class E2PBase(KernelCacheWrapper):
         sympy_conv = SympyToPymbolicMapper()
         return [lp.Assignment(id=None,
                     assignee="kernel_scaling",
-                    expression=sympy_conv(self.expansion.kernel.get_scaling()),
+                    expression=sympy_conv(
+                        self.expansion.kernel.get_global_scaling_const()),
                     temp_var_type=lp.auto)]
 
     def get_cache_key(self):
@@ -169,6 +173,7 @@ class E2PFromSingleBox(E2PBase):
                     lp.GlobalArg("box_target_starts,box_target_counts_nonchild",
                         None, shape=None),
                     lp.GlobalArg("centers", None, shape="dim, naligned_boxes"),
+                    lp.ValueArg("rscale", None),
                     lp.GlobalArg("result", None, shape="nresults, ntargets",
                         dim_tags="sep,C"),
                     lp.GlobalArg("src_expansions", None,
@@ -209,7 +214,12 @@ class E2PFromSingleBox(E2PBase):
         """
         knl = self.get_cached_optimized_kernel()
 
-        return knl(queue, **kwargs)
+        centers = kwargs.pop("centers")
+        # "1" may be passed for rscale, which won't have its type
+        # meaningfully inferred. Make the type of rscale explicit.
+        rscale = centers.dtype.type(kwargs.pop("rscale"))
+
+        return knl(queue, centers=centers, rscale=rscale, **kwargs)
 
 # }}}
 
@@ -304,7 +314,13 @@ class E2PFromCSR(E2PBase):
 
     def __call__(self, queue, **kwargs):
         knl = self.get_cached_optimized_kernel()
-        return knl(queue, **kwargs)
+
+        centers = kwargs.pop("centers")
+        # "1" may be passed for rscale, which won't have its type
+        # meaningfully inferred. Make the type of rscale explicit.
+        rscale = centers.dtype.type(kwargs.pop("rscale"))
+
+        return knl(queue, centers=centers, rscale=rscale, **kwargs)
 
 # }}}
 
diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index b5096e0c572af3bc9e779976d3368d6c377cf5d0..6b9862331448b9d3bcb7c4ae8a4e1f21179cf513 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -1,6 +1,4 @@
-from __future__ import division
-from __future__ import absolute_import
-from six.moves import range
+from __future__ import division, absolute_import
 
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
@@ -24,11 +22,16 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
+from six.moves import range
+import six
 import numpy as np
 import logging
 from pytools import memoize_method
 import sumpy.symbolic as sym
 from sumpy.tools import MiDerivativeTaker
+import sumpy.symbolic as sp
+from collections import defaultdict
+
 
 __doc__ = """
 .. autoclass:: ExpansionBase
@@ -47,8 +50,17 @@ logger = logging.getLogger(__name__)
 # {{{ base class
 
 class ExpansionBase(object):
+    """
+    .. automethod:: with_kernel
+    .. automethod:: __len__
+    .. automethod:: get_coefficient_identifiers
+    .. automethod:: coefficients_from_source
+    .. automethod:: translate_from
+    .. automethod:: __eq__
+    .. automethod:: __ne__
+    """
 
-    def __init__(self, kernel, order):
+    def __init__(self, kernel, order, use_rscale=None):
         # Don't be tempted to remove target derivatives here.
         # Line Taylor QBX can't do without them, because it can't
         # apply those derivatives to the expanded quantity.
@@ -56,8 +68,16 @@ class ExpansionBase(object):
         self.kernel = kernel
         self.order = order
 
+        if use_rscale is None:
+            use_rscale = True
+
+        self.use_rscale = use_rscale
+
     # {{{ propagate kernel interface
 
+    # This is here to conform this to enough of the kernel interface
+    # to make it fit into sumpy.qbx.LayerPotential.
+
     @property
     def dim(self):
         return self.kernel.dim
@@ -72,8 +92,8 @@ class ExpansionBase(object):
     def get_code_transformer(self):
         return self.kernel.get_code_transformer()
 
-    def get_scaling(self):
-        return self.kernel.get_scaling()
+    def get_global_scaling_const(self):
+        return self.kernel.get_global_scaling_const()
 
     def get_args(self):
         return self.kernel.get_args()
@@ -84,13 +104,18 @@ class ExpansionBase(object):
     # }}}
 
     def with_kernel(self, kernel):
-        return type(self)(kernel, self.order)
-        # FIXME: add self.use_rscale once the rscale MR is in
+        return type(self)(kernel, self.order, self.use_rscale)
 
     def __len__(self):
         return len(self.get_coefficient_identifiers())
 
-    def coefficients_from_source(self, avec, bvec):
+    def get_coefficient_identifiers(self):
+        """
+        Returns the identifiers of the coefficients that actually get stored.
+        """
+        raise NotImplementedError
+
+    def coefficients_from_source(self, avec, bvec, rscale):
         """Form an expansion from a source point.
 
         :arg avec: vector from source to center.
@@ -102,7 +127,7 @@ class ExpansionBase(object):
         """
         raise NotImplementedError
 
-    def evaluate(self, coeffs, bvec):
+    def evaluate(self, coeffs, bvec, rscale):
         """
         :return: a :mod:`sympy` expression corresponding
             to the evaluated expansion with the coefficients
@@ -111,16 +136,22 @@ class ExpansionBase(object):
 
         raise NotImplementedError
 
+    def translate_from(self, src_expansion, src_coeff_exprs, src_rscale,
+            dvec, tgt_rscale):
+        raise NotImplementedError
+
     def update_persistent_hash(self, key_hash, key_builder):
         key_hash.update(type(self).__name__.encode("utf8"))
         key_builder.rec(key_hash, self.kernel)
         key_builder.rec(key_hash, self.order)
+        key_builder.rec(key_hash, self.use_rscale)
 
     def __eq__(self, other):
         return (
                 type(self) == type(other)
                 and self.kernel == other.kernel
-                and self.order == other.order)
+                and self.order == other.order
+                and self.use_rscale == other.use_rscale)
 
     def __ne__(self, other):
         return not self.__eq__(other)
@@ -139,10 +170,12 @@ class DerivativeWrangler(object):
     def get_coefficient_identifiers(self):
         raise NotImplementedError
 
-    def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives):
+    def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives,
+            rscale):
         raise NotImplementedError
 
-    def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients):
+    def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients,
+            rscale):
         raise NotImplementedError
 
     def get_derivative_taker(self, expr, var_list):
@@ -169,7 +202,8 @@ class FullDerivativeWrangler(DerivativeWrangler):
     get_coefficient_identifiers = (
             DerivativeWrangler.get_full_coefficient_identifiers)
 
-    def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives):
+    def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives,
+            rscale):
         return stored_kernel_derivatives
 
     get_stored_mpole_coefficients_from_full = (
@@ -212,28 +246,37 @@ def _spmv(spmat, x, sparse_vectors):
 
 
 class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
-
-    def __init__(self, order, dim):
-        DerivativeWrangler.__init__(self, order, dim)
+    _rscale_symbol = sp.Symbol("_sumpy_rscale_placeholder")
 
     def get_coefficient_identifiers(self):
         return self.stored_identifiers
 
-    def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives):
-        return _spmv(self.coeff_matrix, stored_kernel_derivatives,
+    def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives,
+            rscale):
+        coeff_matrix = self.get_coefficient_matrix(rscale)
+        return _spmv(coeff_matrix, stored_kernel_derivatives,
                      sparse_vectors=False)
 
-    def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients):
+    def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients,
+            rscale):
         # = M^T x, where M = coeff matrix
+
+        coeff_matrix = self.get_coefficient_matrix(rscale)
         result = [0] * len(self.stored_identifiers)
         for row, coeff in enumerate(full_mpole_coefficients):
-            for col, val in self.coeff_matrix[row]:
+            for col, val in coeff_matrix[row]:
                 result[col] += coeff * val
         return result
 
-    def precompute_recurrences(self):
+    @property
+    def stored_identifiers(self):
+        stored_identifiers, coeff_matrix = self._get_stored_ids_and_coeff_mat()
+        return stored_identifiers
+
+    @memoize_method
+    def get_coefficient_matrix(self, rscale):
         """
-        Build up a matrix that expresses every derivative in terms of a
+        Return a matrix that expresses every derivative in terms of a
         set of "stored" derivatives.
 
         For example, for the recurrence
@@ -250,6 +293,20 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
 
            ^ rows = one for every derivative
         """
+        stored_identifiers, coeff_matrix = self._get_stored_ids_and_coeff_mat()
+
+        # substitute actual rscale for internal placeholder
+        return defaultdict(lambda: [],
+                ((irow, [
+                    (icol,
+                        coeff.subs(self._rscale_symbol, rscale)
+                        if isinstance(coeff, sp.Basic)
+                        else coeff)
+                    for icol, coeff in row])
+                for irow, row in six.iteritems(coeff_matrix)))
+
+    @memoize_method
+    def _get_stored_ids_and_coeff_mat(self):
         stored_identifiers = []
         identifiers_so_far = {}
 
@@ -258,14 +315,13 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
         logger.debug("computing recurrence for Taylor coefficients: start")
 
         # Sparse matrix, indexed by row
-        from collections import defaultdict
         coeff_matrix_transpose = defaultdict(lambda: [])
 
         # Build up the matrix transpose by row.
         from six import iteritems
         for i, identifier in enumerate(self.get_full_coefficient_identifiers()):
             expr = self.try_get_recurrence_for_derivative(
-                    identifier, identifiers_so_far)
+                    identifier, identifiers_so_far, rscale=self._rscale_symbol)
 
             if expr is None:
                 # Identifier should be stored
@@ -280,7 +336,7 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
 
             identifiers_so_far[identifier] = i
 
-        self.stored_identifiers = stored_identifiers
+        stored_identifiers = stored_identifiers
 
         coeff_matrix = defaultdict(lambda: [])
         for i, row in iteritems(coeff_matrix_transpose):
@@ -293,11 +349,18 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
 
         logger.debug("number of Taylor coefficients was reduced from {orig} to {red}"
                      .format(orig=len(self.get_full_coefficient_identifiers()),
-                             red=len(self.stored_identifiers)))
+                             red=len(stored_identifiers)))
 
-        self.coeff_matrix = coeff_matrix
+        return stored_identifiers, coeff_matrix
+
+    def try_get_recurrence_for_derivative(self, deriv, in_terms_of, rscale):
+        """
+        :arg deriv: a tuple of integers identifying a derivative for which
+            a recurrence is sought
+        :arg in_terms_of: a container supporting efficient containment testing
+            indicating availability of derivatives for use in recurrences
+        """
 
-    def try_get_recurrence_for_derivative(self, deriv, in_terms_of):
         raise NotImplementedError
 
     def get_derivative_taker(self, expr, var_list):
@@ -307,14 +370,10 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
 
 class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
 
-    def __init__(self, order, dim):
-        LinearRecurrenceBasedDerivativeWrangler.__init__(self, order, dim)
-        self.precompute_recurrences()
-
-    def try_get_recurrence_for_derivative(self, deriv, in_terms_of):
+    def try_get_recurrence_for_derivative(self, deriv, in_terms_of, rscale):
         deriv = np.array(deriv, dtype=int)
 
-        for dim in np.nonzero(2 <= deriv)[0]:
+        for dim in np.where(2 <= deriv)[0]:
             # Check if we can reduce this dimension in terms of the other
             # dimensions.
 
@@ -340,14 +399,13 @@ class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
 class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
 
     def __init__(self, order, dim, helmholtz_k_name):
-        LinearRecurrenceBasedDerivativeWrangler.__init__(self, order, dim)
+        super(HelmholtzDerivativeWrangler, self).__init__(order, dim)
         self.helmholtz_k_name = helmholtz_k_name
-        self.precompute_recurrences()
 
-    def try_get_recurrence_for_derivative(self, deriv, in_terms_of):
+    def try_get_recurrence_for_derivative(self, deriv, in_terms_of, rscale):
         deriv = np.array(deriv, dtype=int)
 
-        for dim in np.nonzero(2 <= deriv)[0]:
+        for dim in np.where(2 <= deriv)[0]:
             # Check if we can reduce this dimension in terms of the other
             # dimensions.
 
@@ -368,7 +426,7 @@ class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
                 coeffs[needed_deriv] = -1
             else:
                 k = sym.Symbol(self.helmholtz_k_name)
-                coeffs[tuple(reduced_deriv)] = -k*k
+                coeffs[tuple(reduced_deriv)] = -k*k*rscale*rscale
                 return coeffs
 
 # }}}
@@ -420,7 +478,8 @@ class VolumeTaylorExpansion(VolumeTaylorExpansionBase):
     derivative_wrangler_class = FullDerivativeWrangler
     derivative_wrangler_cache = {}
 
-    def __init__(self, kernel, order):
+    # not user-facing, be strict about having to pass use_rscale
+    def __init__(self, kernel, order, use_rscale):
         self.derivative_wrangler_key = (order, kernel.dim)
 
 
@@ -429,7 +488,8 @@ class LaplaceConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase):
     derivative_wrangler_class = LaplaceDerivativeWrangler
     derivative_wrangler_cache = {}
 
-    def __init__(self, kernel, order):
+    # not user-facing, be strict about having to pass use_rscale
+    def __init__(self, kernel, order, use_rscale):
         self.derivative_wrangler_key = (order, kernel.dim)
 
 
@@ -438,7 +498,8 @@ class HelmholtzConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase):
     derivative_wrangler_class = HelmholtzDerivativeWrangler
     derivative_wrangler_cache = {}
 
-    def __init__(self, kernel, order):
+    # not user-facing, be strict about having to pass use_rscale
+    def __init__(self, kernel, order, use_rscale):
         helmholtz_k_name = kernel.get_base_kernel().helmholtz_k_name
         self.derivative_wrangler_key = (order, kernel.dim, helmholtz_k_name)
 
@@ -499,8 +560,7 @@ class DefaultExpansionFactory(ExpansionFactoryBase):
                 and base_kernel.dim == 2):
             from sumpy.expansion.local import Y2DLocalExpansion
             return Y2DLocalExpansion
-        elif (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel)
-                and base_kernel.dim == 3):
+        elif isinstance(base_kernel.get_base_kernel(), HelmholtzKernel):
             from sumpy.expansion.local import \
                     HelmholtzConformingVolumeTaylorLocalExpansion
             return HelmholtzConformingVolumeTaylorLocalExpansion
@@ -528,8 +588,7 @@ class DefaultExpansionFactory(ExpansionFactoryBase):
             from sumpy.expansion.multipole import (
                     LaplaceConformingVolumeTaylorMultipoleExpansion)
             return LaplaceConformingVolumeTaylorMultipoleExpansion
-        elif (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel)
-                and base_kernel.dim == 3):
+        elif isinstance(base_kernel.get_base_kernel(), HelmholtzKernel):
             from sumpy.expansion.multipole import (
                     HelmholtzConformingVolumeTaylorMultipoleExpansion)
             return HelmholtzConformingVolumeTaylorMultipoleExpansion
diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py
index 76f3d1ebaaa6a12037d333559deb49ffc8c19e48..9e891a2dc48399bfded9e7cb0fc6ed63f38c40e4 100644
--- a/sumpy/expansion/local.py
+++ b/sumpy/expansion/local.py
@@ -1,5 +1,4 @@
 from __future__ import division, absolute_import
-from six.moves import range, zip
 
 __copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
 
@@ -23,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
+from six.moves import range, zip
 import sumpy.symbolic as sym
 
 from sumpy.expansion import (
@@ -57,7 +57,8 @@ class LineTaylorLocalExpansion(LocalExpansionBase):
     def get_coefficient_identifiers(self):
         return list(range(self.order+1))
 
-    def coefficients_from_source(self, avec, bvec):
+    def coefficients_from_source(self, avec, bvec, rscale):
+        # no point in heeding rscale here--just ignore it
         if bvec is None:
             raise RuntimeError("cannot use line-Taylor expansions in a setting "
                     "where the center-target vector is not known at coefficient "
@@ -97,7 +98,8 @@ class LineTaylorLocalExpansion(LocalExpansionBase):
                     .subs("tau", 0)
                     for i in self.get_coefficient_identifiers()]
 
-    def evaluate(self, coeffs, bvec):
+    def evaluate(self, coeffs, bvec, rscale):
+        # no point in heeding rscale here--just ignore it
         from pytools import factorial
         return sym.Add(*(
                 coeffs[self.get_storage_index(i)] / factorial(i)
@@ -113,33 +115,42 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
     Coefficients represent derivative values of the kernel.
     """
 
-    def coefficients_from_source(self, avec, bvec):
+    def coefficients_from_source(self, avec, bvec, rscale):
         from sumpy.tools import MiDerivativeTaker
         ppkernel = self.kernel.postprocess_at_source(
                 self.kernel.get_expression(avec), avec)
 
         taker = MiDerivativeTaker(ppkernel, avec)
-        return [taker.diff(mi) for mi in self.get_coefficient_identifiers()]
+        return [
+                taker.diff(mi) * rscale ** sum(mi)
+                for mi in self.get_coefficient_identifiers()]
 
-    def evaluate(self, coeffs, bvec):
+    def evaluate(self, coeffs, bvec, rscale):
         from sumpy.tools import mi_power, mi_factorial
         evaluated_coeffs = (
-            self.derivative_wrangler.get_full_kernel_derivatives_from_stored(coeffs))
+            self.derivative_wrangler.get_full_kernel_derivatives_from_stored(
+                coeffs, rscale))
+        bvec = bvec / rscale
         result = sum(
                 coeff
-                * self.kernel.postprocess_at_target(mi_power(bvec, mi), bvec)
+                * mi_power(bvec, mi)
                 / mi_factorial(mi)
                 for coeff, mi in zip(
                         evaluated_coeffs, self.get_full_coefficient_identifiers()))
         return result
 
-    def translate_from(self, src_expansion, src_coeff_exprs, dvec):
+    def translate_from(self, src_expansion, src_coeff_exprs, src_rscale,
+            dvec, tgt_rscale):
         logger.info("building translation operator: %s(%d) -> %s(%d): start"
                 % (type(src_expansion).__name__,
                     src_expansion.order,
                     type(self).__name__,
                     self.order))
 
+        if not self.use_rscale:
+            src_rscale = 1
+            tgt_rscale = 1
+
         from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansionBase
         if isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase):
             # We know the general form of the multipole expansion is:
@@ -151,6 +162,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
             #
             # This code speeds up derivative taking by caching all kernel
             # derivatives.
+
             taker = src_expansion.get_kernel_derivative_taker(dvec)
 
             from sumpy.tools import add_mi
@@ -161,44 +173,62 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
                 for coeff, term in zip(
                         src_coeff_exprs,
                         src_expansion.get_coefficient_identifiers()):
-                    kernel_deriv = taker.diff(add_mi(deriv, term))
-                    local_result.append(coeff * kernel_deriv)
+
+                    kernel_deriv = (
+                            src_expansion.get_scaled_multipole(
+                                taker.diff(add_mi(deriv, term)),
+                                dvec, src_rscale,
+                                nderivatives=sum(deriv) + sum(term),
+                                nderivatives_for_scaling=sum(term)))
+
+                    local_result.append(
+                            coeff * kernel_deriv * tgt_rscale**sum(deriv))
                 result.append(sym.Add(*local_result))
         else:
             from sumpy.tools import MiDerivativeTaker
-            expr = src_expansion.evaluate(src_coeff_exprs, dvec)
+            expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=src_rscale)
             taker = MiDerivativeTaker(expr, dvec)
-            result = [taker.diff(mi) for mi in self.get_coefficient_identifiers()]
+
+            # Rscale/operand magnitude is fairly sensitive to the order of
+            # operations--which is something we don't have fantastic control
+            # over at the symbolic level. The '.expand()' below moves the two
+            # canceling "rscales" closer to each other in the hope of helping
+            # with that.
+            result = [
+                    (taker.diff(mi) * tgt_rscale**sum(mi)).expand(deep=False)
+                    for mi in self.get_coefficient_identifiers()]
 
         logger.info("building translation operator: done")
         return result
 
 
 class VolumeTaylorLocalExpansion(
-        VolumeTaylorLocalExpansionBase,
-        VolumeTaylorExpansion):
+        VolumeTaylorExpansion,
+        VolumeTaylorLocalExpansionBase):
 
-    def __init__(self, kernel, order):
-        VolumeTaylorLocalExpansionBase.__init__(self, kernel, order)
-        VolumeTaylorExpansion.__init__(self, kernel, order)
+    def __init__(self, kernel, order, use_rscale=None):
+        VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale)
+        VolumeTaylorExpansion.__init__(self, kernel, order, use_rscale)
 
 
 class LaplaceConformingVolumeTaylorLocalExpansion(
-        VolumeTaylorLocalExpansionBase,
-        LaplaceConformingVolumeTaylorExpansion):
+        LaplaceConformingVolumeTaylorExpansion,
+        VolumeTaylorLocalExpansionBase):
 
-    def __init__(self, kernel, order):
-        VolumeTaylorLocalExpansionBase.__init__(self, kernel, order)
-        LaplaceConformingVolumeTaylorExpansion.__init__(self, kernel, order)
+    def __init__(self, kernel, order, use_rscale=None):
+        VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale)
+        LaplaceConformingVolumeTaylorExpansion.__init__(
+                self, kernel, order, use_rscale)
 
 
 class HelmholtzConformingVolumeTaylorLocalExpansion(
-        VolumeTaylorLocalExpansionBase,
-        HelmholtzConformingVolumeTaylorExpansion):
+        HelmholtzConformingVolumeTaylorExpansion,
+        VolumeTaylorLocalExpansionBase):
 
-    def __init__(self, kernel, order):
-        VolumeTaylorLocalExpansionBase.__init__(self, kernel, order)
-        HelmholtzConformingVolumeTaylorExpansion.__init__(self, kernel, order)
+    def __init__(self, kernel, order, use_rscale=None):
+        VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale)
+        HelmholtzConformingVolumeTaylorExpansion.__init__(
+                self, kernel, order, use_rscale)
 
 # }}}
 
@@ -212,7 +242,10 @@ class _FourierBesselLocalExpansion(LocalExpansionBase):
     def get_coefficient_identifiers(self):
         return list(range(-self.order, self.order+1))
 
-    def coefficients_from_source(self, avec, bvec):
+    def coefficients_from_source(self, avec, bvec, rscale):
+        if not self.use_rscale:
+            rscale = 1
+
         from sumpy.symbolic import sym_real_norm_2
         hankel_1 = sym.Function("hankel_1")
 
@@ -223,10 +256,14 @@ class _FourierBesselLocalExpansion(LocalExpansionBase):
         avec_len = sym_real_norm_2(avec)
         return [self.kernel.postprocess_at_source(
                     hankel_1(l, arg_scale * avec_len)
+                    * rscale ** abs(l)
                     * sym.exp(sym.I * l * source_angle_rel_center), avec)
                     for l in self.get_coefficient_identifiers()]
 
-    def evaluate(self, coeffs, bvec):
+    def evaluate(self, coeffs, bvec, rscale):
+        if not self.use_rscale:
+            rscale = 1
+
         from sumpy.symbolic import sym_real_norm_2
         bessel_j = sym.Function("bessel_j")
         bvec_len = sym_real_norm_2(bvec)
@@ -237,12 +274,18 @@ class _FourierBesselLocalExpansion(LocalExpansionBase):
         return sum(coeffs[self.get_storage_index(l)]
                    * self.kernel.postprocess_at_target(
                        bessel_j(l, arg_scale * bvec_len)
+                       / rscale ** abs(l)
                        * sym.exp(sym.I * l * -target_angle_rel_center), bvec)
                 for l in self.get_coefficient_identifiers())
 
-    def translate_from(self, src_expansion, src_coeff_exprs, dvec):
+    def translate_from(self, src_expansion, src_coeff_exprs, src_rscale,
+            dvec, tgt_rscale):
         from sumpy.symbolic import sym_real_norm_2
 
+        if not self.use_rscale:
+            src_rscale = 1
+            tgt_rscale = 1
+
         arg_scale = self.get_bessel_arg_scaling()
 
         if isinstance(src_expansion, type(self)):
@@ -254,6 +297,8 @@ class _FourierBesselLocalExpansion(LocalExpansionBase):
                 translated_coeffs.append(
                     sum(src_coeff_exprs[src_expansion.get_storage_index(m)]
                         * bessel_j(m - l, arg_scale * dvec_len)
+                        / src_rscale ** abs(m)
+                        * tgt_rscale ** abs(l)
                         * sym.exp(sym.I * (m - l) * -new_center_angle_rel_old_center)
                     for m in src_expansion.get_coefficient_identifiers()))
             return translated_coeffs
@@ -265,10 +310,14 @@ class _FourierBesselLocalExpansion(LocalExpansionBase):
             translated_coeffs = []
             for l in self.get_coefficient_identifiers():
                 translated_coeffs.append(
-                    sum((-1) ** l * hankel_1(m + l, arg_scale * dvec_len)
+                    sum(
+                        (-1) ** l
+                        * hankel_1(m + l, arg_scale * dvec_len)
+                        * src_rscale ** abs(m)
+                        * tgt_rscale ** abs(l)
                         * sym.exp(sym.I * (m + l) * new_center_angle_rel_old_center)
                         * src_coeff_exprs[src_expansion.get_storage_index(m)]
-                    for m in src_expansion.get_coefficient_identifiers()))
+                        for m in src_expansion.get_coefficient_identifiers()))
             return translated_coeffs
 
         raise RuntimeError("do not know how to translate %s to %s"
@@ -277,12 +326,12 @@ class _FourierBesselLocalExpansion(LocalExpansionBase):
 
 
 class H2DLocalExpansion(_FourierBesselLocalExpansion):
-    def __init__(self, kernel, order):
+    def __init__(self, kernel, order, use_rscale=None):
         from sumpy.kernel import HelmholtzKernel
         assert (isinstance(kernel.get_base_kernel(), HelmholtzKernel)
                 and kernel.dim == 2)
 
-        super(H2DLocalExpansion, self).__init__(kernel, order)
+        super(H2DLocalExpansion, self).__init__(kernel, order, use_rscale)
 
         from sumpy.expansion.multipole import H2DMultipoleExpansion
         self.mpole_expn_class = H2DMultipoleExpansion
@@ -292,12 +341,12 @@ class H2DLocalExpansion(_FourierBesselLocalExpansion):
 
 
 class Y2DLocalExpansion(_FourierBesselLocalExpansion):
-    def __init__(self, kernel, order):
+    def __init__(self, kernel, order, use_rscale=None):
         from sumpy.kernel import YukawaKernel
         assert (isinstance(kernel.get_base_kernel(), YukawaKernel)
                 and kernel.dim == 2)
 
-        super(Y2DLocalExpansion, self).__init__(kernel, order)
+        super(Y2DLocalExpansion, self).__init__(kernel, order, use_rscale)
 
         from sumpy.expansion.multipole import Y2DMultipoleExpansion
         self.mpole_expn_class = Y2DMultipoleExpansion
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index 958bca6dbec6b8fdcbbb77545f039ae92a65d3d1..4db9250b329602ead272e122c0c128059829277d 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -25,6 +25,7 @@ THE SOFTWARE.
 from six.moves import range, zip
 import sumpy.symbolic as sym  # noqa
 
+from sumpy.symbolic import vector_subs
 from sumpy.expansion import (
     ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion,
     HelmholtzConformingVolumeTaylorExpansion)
@@ -53,12 +54,15 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
     Coefficients represent the terms in front of the kernel derivatives.
     """
 
-    def coefficients_from_source(self, avec, bvec):
+    def coefficients_from_source(self, avec, bvec, rscale):
         from sumpy.kernel import DirectionalSourceDerivative
         kernel = self.kernel
 
         from sumpy.tools import mi_power, mi_factorial
 
+        if not self.use_rscale:
+            rscale = 1
+
         if isinstance(kernel, DirectionalSourceDerivative):
             if kernel.get_base_kernel() is not kernel.inner_kernel:
                 raise NotImplementedError("more than one source derivative "
@@ -82,33 +86,61 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
                         - mi_power(avec, derivative_mi) * mi[idim]
                         * dir_vec[idim])
             for i, mi in enumerate(coeff_identifiers):
-                result[i] /= mi_factorial(mi)
+                result[i] /= (mi_factorial(mi) * rscale ** sum(mi))
         else:
+            avec = avec/rscale
+
             result = [
                     mi_power(avec, mi) / mi_factorial(mi)
                     for mi in self.get_full_coefficient_identifiers()]
         return (
-            self.derivative_wrangler.get_stored_mpole_coefficients_from_full(result))
+            self.derivative_wrangler.get_stored_mpole_coefficients_from_full(
+                result, rscale))
+
+    def get_scaled_multipole(self, expr, bvec, rscale, nderivatives,
+            nderivatives_for_scaling=None):
+        if nderivatives_for_scaling is None:
+            nderivatives_for_scaling = nderivatives
+
+        if self.kernel.has_efficient_scale_adjustment:
+            return (
+                    self.kernel.adjust_for_kernel_scaling(
+                        vector_subs(
+                            expr,
+                            bvec, bvec/rscale),
+                        rscale, nderivatives)
+                    / rscale ** (nderivatives - nderivatives_for_scaling))
+        else:
+            return (rscale**nderivatives_for_scaling * expr)
+
+    def evaluate(self, coeffs, bvec, rscale):
+        if not self.use_rscale:
+            rscale = 1
 
-    def evaluate(self, coeffs, bvec):
         taker = self.get_kernel_derivative_taker(bvec)
+
         result = sym.Add(*tuple(
-                coeff * taker.diff(mi)
+                coeff
+                * self.get_scaled_multipole(taker.diff(mi), bvec, rscale, sum(mi))
                 for coeff, mi in zip(coeffs, self.get_coefficient_identifiers())))
+
         return result
 
     def get_kernel_derivative_taker(self, bvec):
-        ppkernel = self.kernel.postprocess_at_target(
-                self.kernel.get_expression(bvec), bvec)
+        return (self.derivative_wrangler.get_derivative_taker(
+            self.kernel.get_expression(bvec), bvec))
 
-        return self.derivative_wrangler.get_derivative_taker(ppkernel, bvec)
-
-    def translate_from(self, src_expansion, src_coeff_exprs, dvec):
+    def translate_from(self, src_expansion, src_coeff_exprs, src_rscale,
+            dvec, tgt_rscale):
         if not isinstance(src_expansion, type(self)):
             raise RuntimeError("do not know how to translate %s to "
                     "Taylor multipole expansion"
                                % type(src_expansion).__name__)
 
+        if not self.use_rscale:
+            src_rscale = 1
+            tgt_rscale = 1
+
         logger.info("building translation operator: %s(%d) -> %s(%d): start"
                 % (type(src_expansion).__name__,
                     src_expansion.order,
@@ -148,40 +180,45 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
                     contrib *= (binomial(n, k)
                             * dvec[idim]**(n-k))
 
-                result[i] += contrib
+                result[i] += (
+                        contrib
+                        * (src_rscale**sum(src_mi) / tgt_rscale**sum(tgt_mi)))
 
             result[i] /= mi_factorial(tgt_mi)
 
         logger.info("building translation operator: done")
         return (
-            self.derivative_wrangler.get_stored_mpole_coefficients_from_full(result))
+            self.derivative_wrangler.get_stored_mpole_coefficients_from_full(
+                result, tgt_rscale))
 
 
 class VolumeTaylorMultipoleExpansion(
-        VolumeTaylorMultipoleExpansionBase,
-        VolumeTaylorExpansion):
+        VolumeTaylorExpansion,
+        VolumeTaylorMultipoleExpansionBase):
 
-    def __init__(self, kernel, order):
-        VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order)
-        VolumeTaylorExpansion.__init__(self, kernel, order)
+    def __init__(self, kernel, order, use_rscale=None):
+        VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale)
+        VolumeTaylorExpansion.__init__(self, kernel, order, use_rscale)
 
 
 class LaplaceConformingVolumeTaylorMultipoleExpansion(
-        VolumeTaylorMultipoleExpansionBase,
-        LaplaceConformingVolumeTaylorExpansion):
+        LaplaceConformingVolumeTaylorExpansion,
+        VolumeTaylorMultipoleExpansionBase):
 
-    def __init__(self, kernel, order):
-        VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order)
-        LaplaceConformingVolumeTaylorExpansion.__init__(self, kernel, order)
+    def __init__(self, kernel, order, use_rscale=None):
+        VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale)
+        LaplaceConformingVolumeTaylorExpansion.__init__(
+                self, kernel, order, use_rscale)
 
 
 class HelmholtzConformingVolumeTaylorMultipoleExpansion(
-        VolumeTaylorMultipoleExpansionBase,
-        HelmholtzConformingVolumeTaylorExpansion):
+        HelmholtzConformingVolumeTaylorExpansion,
+        VolumeTaylorMultipoleExpansionBase):
 
-    def __init__(self, kernel, order):
-        VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order)
-        HelmholtzConformingVolumeTaylorExpansion.__init__(self, kernel, order)
+    def __init__(self, kernel, order, use_rscale=None):
+        VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale)
+        HelmholtzConformingVolumeTaylorExpansion.__init__(
+                self, kernel, order, use_rscale)
 
 # }}}
 
@@ -195,7 +232,10 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase):
     def get_coefficient_identifiers(self):
         return list(range(-self.order, self.order+1))
 
-    def coefficients_from_source(self, avec, bvec):
+    def coefficients_from_source(self, avec, bvec, rscale):
+        if not self.use_rscale:
+            rscale = 1
+
         from sumpy.symbolic import sym_real_norm_2
         bessel_j = sym.Function("bessel_j")
         avec_len = sym_real_norm_2(avec)
@@ -204,12 +244,18 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase):
 
         # The coordinates are negated since avec points from source to center.
         source_angle_rel_center = sym.atan2(-avec[1], -avec[0])
-        return [self.kernel.postprocess_at_source(
-                    bessel_j(l, arg_scale * avec_len) *
-                    sym.exp(sym.I * l * -source_angle_rel_center), avec)
+        return [
+                self.kernel.postprocess_at_source(
+                    bessel_j(l, arg_scale * avec_len)
+                    / rscale ** abs(l)
+                    * sym.exp(sym.I * l * -source_angle_rel_center),
+                    avec)
                 for l in self.get_coefficient_identifiers()]
 
-    def evaluate(self, coeffs, bvec):
+    def evaluate(self, coeffs, bvec, rscale):
+        if not self.use_rscale:
+            rscale = 1
+
         from sumpy.symbolic import sym_real_norm_2
         hankel_1 = sym.Function("hankel_1")
         bvec_len = sym_real_norm_2(bvec)
@@ -220,14 +266,21 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase):
         return sum(coeffs[self.get_storage_index(l)]
                    * self.kernel.postprocess_at_target(
                        hankel_1(l, arg_scale * bvec_len)
+                       * rscale ** abs(l)
                        * sym.exp(sym.I * l * target_angle_rel_center), bvec)
                 for l in self.get_coefficient_identifiers())
 
-    def translate_from(self, src_expansion, src_coeff_exprs, dvec):
+    def translate_from(self, src_expansion, src_coeff_exprs, src_rscale,
+            dvec, tgt_rscale):
         if not isinstance(src_expansion, type(self)):
             raise RuntimeError("do not know how to translate %s to %s"
                                % (type(src_expansion).__name__,
                                    type(self).__name__))
+
+        if not self.use_rscale:
+            src_rscale = 1
+            tgt_rscale = 1
+
         from sumpy.symbolic import sym_real_norm_2
         dvec_len = sym_real_norm_2(dvec)
         bessel_j = sym.Function("bessel_j")
@@ -240,30 +293,34 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase):
             translated_coeffs.append(
                 sum(src_coeff_exprs[src_expansion.get_storage_index(m)]
                     * bessel_j(m - l, arg_scale * dvec_len)
+                    * src_rscale ** abs(m)
+                    / tgt_rscale ** abs(l)
                     * sym.exp(sym.I * (m - l) * new_center_angle_rel_old_center)
                 for m in src_expansion.get_coefficient_identifiers()))
         return translated_coeffs
 
 
 class H2DMultipoleExpansion(_HankelBased2DMultipoleExpansion):
-    def __init__(self, kernel, order):
+    def __init__(self, kernel, order, use_rscale=None):
         from sumpy.kernel import HelmholtzKernel
         assert (isinstance(kernel.get_base_kernel(), HelmholtzKernel)
                 and kernel.dim == 2)
 
-        super(H2DMultipoleExpansion, self).__init__(kernel, order)
+        super(H2DMultipoleExpansion, self).__init__(
+                kernel, order, use_rscale=use_rscale)
 
     def get_bessel_arg_scaling(self):
         return sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name)
 
 
 class Y2DMultipoleExpansion(_HankelBased2DMultipoleExpansion):
-    def __init__(self, kernel, order):
+    def __init__(self, kernel, order, use_rscale=None):
         from sumpy.kernel import YukawaKernel
         assert (isinstance(kernel.get_base_kernel(), YukawaKernel)
                 and kernel.dim == 2)
 
-        super(Y2DMultipoleExpansion, self).__init__(kernel, order)
+        super(Y2DMultipoleExpansion, self).__init__(
+                kernel, order, use_rscale=use_rscale)
 
     def get_bessel_arg_scaling(self):
         return sym.I * sym.Symbol(self.kernel.get_base_kernel().yukawa_lambda_name)
diff --git a/sumpy/fmm.py b/sumpy/fmm.py
index 77226a9a4df94d72f6c95cad39ecd8e103bd486c..7c1894c43eb9af122e085f6e8c2d15ca5df7a249 100644
--- a/sumpy/fmm.py
+++ b/sumpy/fmm.py
@@ -43,6 +43,10 @@ from sumpy import (
         E2EFromCSR, E2EFromChildren, E2EFromParent)
 
 
+def level_to_rscale(tree, level):
+    return tree.root_extent * (2**-level)
+
+
 # {{{ expansion wrangler code container
 
 class SumpyExpansionWranglerCodeContainer(object):
@@ -56,7 +60,7 @@ class SumpyExpansionWranglerCodeContainer(object):
     def __init__(self, cl_context,
             multipole_expansion_factory,
             local_expansion_factory,
-            out_kernels, exclude_self=False):
+            out_kernels, exclude_self=False, use_rscale=None):
         """
         :arg multipole_expansion_factory: a callable of a single argument (order)
             that returns a multipole expansion.
@@ -69,16 +73,17 @@ class SumpyExpansionWranglerCodeContainer(object):
         self.local_expansion_factory = local_expansion_factory
         self.out_kernels = out_kernels
         self.exclude_self = exclude_self
+        self.use_rscale = use_rscale
 
         self.cl_context = cl_context
 
     @memoize_method
     def multipole_expansion(self, order):
-        return self.multipole_expansion_factory(order)
+        return self.multipole_expansion_factory(order, self.use_rscale)
 
     @memoize_method
     def local_expansion(self, order):
-        return self.local_expansion_factory(order)
+        return self.local_expansion_factory(order, self.use_rscale)
 
     @memoize_method
     def p2m(self, tgt_order):
@@ -324,6 +329,7 @@ class SumpyExpansionWrangler(object):
                     tgt_base_ibox=level_start_ibox,
 
                     wait_for=src_weights.events,
+                    rscale=level_to_rscale(self.tree, lev),
 
                     **kwargs)
 
@@ -343,18 +349,23 @@ class SumpyExpansionWrangler(object):
         tree = self.tree
 
         evt = None
-        # 2 is the last relevant source_level.
-        # 1 is the last relevant target_level.
-        # (Nobody needs a multipole on level 0, i.e. for the root box.)
-        for source_level in range(tree.nlevels-1, 1, -1):
-            start, stop = level_start_source_parent_box_nrs[
-                            source_level:source_level+2]
-            if start == stop:
-                continue
 
+        # nlevels-1 is the last valid level index
+        # nlevels-2 is the last valid level that could have children
+        #
+        # 3 is the last relevant source_level.
+        # 2 is the last relevant target_level.
+        # (because no level 1 box will be well-separated from another)
+        for source_level in range(tree.nlevels-1, 2, -1):
             target_level = source_level - 1
             assert target_level > 0
 
+            start, stop = level_start_source_parent_box_nrs[
+                            target_level:target_level+2]
+            if start == stop:
+                print("source", source_level, "empty")
+                continue
+
             m2m = self.code.m2m(
                     self.level_orders[source_level],
                     self.level_orders[target_level])
@@ -375,13 +386,17 @@ class SumpyExpansionWrangler(object):
                     box_child_ids=self.tree.box_child_ids,
                     centers=self.tree.box_centers,
 
-                    **self.kernel_extra_kwargs)
+                    src_rscale=level_to_rscale(self.tree, source_level),
+                    tgt_rscale=level_to_rscale(self.tree, target_level),
 
+                    **self.kernel_extra_kwargs)
             assert mpoles_res is target_mpoles_view
 
         if evt is not None:
             mpoles.add_event(evt)
 
+        # }}}
+
         return mpoles
 
     def eval_direct(self, target_boxes, source_box_starts,
@@ -441,6 +456,9 @@ class SumpyExpansionWrangler(object):
                     src_box_lists=src_box_lists,
                     centers=self.tree.box_centers,
 
+                    src_rscale=level_to_rscale(self.tree, lev),
+                    tgt_rscale=level_to_rscale(self.tree, lev),
+
                     wait_for=mpole_exps.events,
 
                     **self.kernel_extra_kwargs)
@@ -484,6 +502,8 @@ class SumpyExpansionWrangler(object):
                     centers=self.tree.box_centers,
                     result=pot,
 
+                    rscale=level_to_rscale(self.tree, isrc_level),
+
                     wait_for=wait_for,
 
                     **kwargs)
@@ -530,6 +550,8 @@ class SumpyExpansionWrangler(object):
                     tgt_expansions=target_local_exps_view,
                     tgt_base_ibox=target_level_start_ibox,
 
+                    rscale=level_to_rscale(self.tree, lev),
+
                     wait_for=src_weights.events,
 
                     **kwargs)
@@ -572,6 +594,9 @@ class SumpyExpansionWrangler(object):
                     box_parent_ids=self.tree.box_parent_ids,
                     centers=self.tree.box_centers,
 
+                    src_rscale=level_to_rscale(self.tree, source_lev),
+                    tgt_rscale=level_to_rscale(self.tree, target_lev),
+
                     **self.kernel_extra_kwargs)
 
             assert local_exps_res is target_local_exps_view
@@ -607,6 +632,8 @@ class SumpyExpansionWrangler(object):
                     centers=self.tree.box_centers,
                     result=pot,
 
+                    rscale=level_to_rscale(self.tree, lev),
+
                     wait_for=local_exps.events,
 
                     **kwargs)
diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index 780cf8eb4dd2eb482ea0492779232bbff92c75c3..ddfba6829ca6c1fdb4bfc6930cc9352e7c1d89e0 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -31,8 +31,8 @@ from sumpy.symbolic import pymbolic_real_norm_2
 from pymbolic.primitives import make_sym_vector
 from pymbolic import var
 
-__doc__ = """
 
+__doc__ = """
 Kernel interface
 ----------------
 
@@ -48,6 +48,7 @@ PDE kernels
 -----------
 
 .. autoclass:: LaplaceKernel
+.. autoclass:: BiharmonicKernel
 .. autoclass:: HelmholtzKernel
 .. autoclass:: YukawaKernel
 .. autoclass:: StokesletKernel
@@ -101,12 +102,21 @@ class Kernel(object):
     .. attribute:: is_complex_valued
     .. attribute:: dim
 
-        *dim* is allowed to be *None* if the dimensionality is not yet
-        known.
+    .. automethod:: get_base_kernel
+    .. automethod:: prepare_loopy_kernel
+    .. automethod:: get_code_transformer
+    .. automethod:: get_expression
+    .. attribute:: has_efficient_scale_adjustment
+    .. automethod:: adjust_for_kernel_scaling
+    .. automethod:: postprocess_at_source
+    .. automethod:: postprocess_at_target
+    .. automethod:: get_global_scaling_const
+    .. automethod:: get_args
+    .. automethod:: get_source_args
     """
 
     def __init__(self, dim=None):
-        self._dim = dim
+        self.dim = dim
 
     # {{{ hashing/pickling/equality
 
@@ -140,21 +150,10 @@ class Kernel(object):
         # Can't use trivial pickling: hash_value cache must stay unset
         assert len(self.init_arg_names) == len(state)
         for name, value in zip(self.init_arg_names, state):
-            if name == "dim":
-                name = "_dim"
-
             setattr(self, name, value)
 
     # }}}
 
-    @property
-    def dim(self):
-        if self._dim is None:
-            raise RuntimeError("the number of dimensions for this kernel "
-                    "has not yet been set")
-
-        return self._dim
-
     def get_base_kernel(self):
         """Return the kernel being wrapped by this one, or else
         *self*.
@@ -175,18 +174,87 @@ class Kernel(object):
         return lambda expr: expr
 
     def get_expression(self, dist_vec):
-        """Return a :mod:`sympy` expression for the kernel.
+        r"""Return a :mod:`sympy` expression for the kernel."""
+        raise NotImplementedError
+
+    has_efficient_scale_adjustment = False
+
+    def adjust_for_kernel_scaling(self, expr, rscale, nderivatives):
+        r"""
+        Consider a Taylor multipole expansion:
+
+        .. math::
+
+            f (x - y) = \sum_{i = 0}^{\infty} (\partial_y^i f) (x - y) \big|_{y = c}
+           \frac{(y - c)^i}{i!} .
+
+        Now suppose we would like to use a scaled version :math:`g` of the
+        kernel :math:`f`:
+
+        .. math::
+
+            \begin{eqnarray*}
+              f (x) & = & g (x / \alpha),\\
+              f^{(i)} (x) & = & \frac{1}{\alpha^i} g^{(i)} (x / \alpha) .
+            \end{eqnarray*}
+
+        where :math:`\alpha` is chosen to be on a length scale similar to
+        :math:`x` (for example by choosing :math:`\alpha` proporitional to the
+        size of the box for which the expansion is intended) so that :math:`x /
+        \alpha` is roughly of unit magnitude, to avoid arithmetic issues with
+        small arguments. This yields
+
+        .. math::
+
+            f (x - y) = \sum_{i = 0}^{\infty} (\partial_y^i g)
+            \left( \frac{x - y}{\alpha} \right) \Bigg|_{y = c}
+            \cdot
+            \frac{(y - c)^i}{\alpha^i \cdot i!}.
+
+        Observe that the :math:`(y - c)` term is now scaled to unit magnitude,
+        as is the argument of :math:`g`.
+
+        With :math:`\xi = x / \alpha`, we find
+
+        .. math::
+
+            \begin{eqnarray*}
+              g (\xi) & = & f (\alpha \xi),\\
+              g^{(i)} (\xi) & = & \alpha^i f^{(i)} (\alpha \xi) .
+            \end{eqnarray*}
+
+        Generically for all kernels, :math:`f^{(i)} (\alpha \xi)` is computable
+        by taking a sufficient number of symbolic derivatives of :math:`f` and
+        providing :math:`\alpha \xi = x` as the argument.
+
+        Now, for some kernels, like :math:`f (x) = C \log x`, the powers of
+        :math:`\alpha^i` from the chain rule cancel with the ones from the
+        argument substituted into the kernel derivative:
 
-        :arg dist_vec: target - source
+        .. math::
 
-        (Assumes translation invariance of the kernel.)
+            g^{(i)} (\xi) = \alpha^i f^{(i)} (\alpha \xi) = C' \cdot \alpha^i \cdot
+            \frac{1}{(\alpha x)^i} \quad (i > 0),
+
+        making them what you might call *scale-invariant*. In this case, one
+        may set :attr:`has_efficient_scale_adjustment`. For these kernels only,
+        :meth:`adjust_for_kernel_scaling` provides a shortcut for scaled kernel
+        derivative computation. Given :math:`f^{(i)}` as *expr*, it directly
+        returns an expression for :math:`g^{(i)}`, where :math:`i` is given
+        as *nderivatives*.
+
+        :arg rscale: The scaling parameter :math:`\alpha` above.
         """
+
         raise NotImplementedError
 
     def postprocess_at_source(self, expr, avec):
         """Transform a kernel evaluation or expansion expression in a place
         where the vector a (something - source) is known. ("something" may be
         an expansion center or a target.)
+
+        The typical use of this function is to apply source-variable
+        derivatives to the kernel.
         """
         return expr
 
@@ -194,11 +262,20 @@ class Kernel(object):
         """Transform a kernel evaluation or expansion expression in a place
         where the vector b (target - something) is known. ("something" may be
         an expansion center or a target.)
+
+        The typical use of this function is to apply target-variable
+        derivatives to the kernel.
         """
         return expr
 
-    def get_scaling(self):
-        """Return a global scaling of the kernel."""
+    def get_global_scaling_const(self):
+        """Return a global scaling constant of the kernel.
+         Typically, this ensures that the kernel is scaled so that
+         :math:`\mathcal L(G)(x)=C\delta(x)` with a constant of 1, where
+         :math:`\mathcal L` is the PDE operator associated with the kernel.
+         Not to be confused with *rscale*, which keeps expansion
+         coefficients benignly scaled.
+        """
         raise NotImplementedError
 
     def get_args(self):
@@ -220,68 +297,62 @@ class Kernel(object):
 class ExpressionKernel(Kernel):
     is_complex_valued = False
 
-    init_arg_names = ("dim", "expression", "scaling", "is_complex_valued")
+    init_arg_names = ("dim", "expression", "global_scaling_const",
+            "is_complex_valued")
 
-    def __init__(self, dim, expression, scaling, is_complex_valued):
+    def __init__(self, dim, expression, global_scaling_const,
+            is_complex_valued):
         """
         :arg expression: A :mod:`pymbolic` expression depending on
             variables *d_1* through *d_N* where *N* equals *dim*.
             (These variables match what is returned from
             :func:`pymbolic.primitives.make_sym_vector` with
             argument `"d"`.)
-        :arg scaling: A :mod:`pymbolic` expression for the scaling
-            of the kernel.
+        :arg global_scaling_const: A constant :mod:`pymbolic` expression for the
+            global scaling of the kernel. Typically, this ensures that
+            the kernel is scaled so that :math:`\mathcal L(G)(x)=C\delta(x)`
+            with a constant of 1, where :math:`\mathcal L` is the PDE
+            operator associated with the kernel. Not to be confused with
+            *rscale*, which keeps expansion coefficients benignly scaled.
         """
 
-        # expression and scaling are pymbolic objects because those pickle
-        # cleanly. D'oh, sympy!
+        # expression and global_scaling_const are pymbolic objects because
+        # those pickle cleanly. D'oh, sympy!
 
         Kernel.__init__(self, dim)
 
         self.expression = expression
-        self.scaling = scaling
+        self.global_scaling_const = global_scaling_const
         self.is_complex_valued = is_complex_valued
 
     def __getinitargs__(self):
-        return (self._dim, self.expression, self.scaling,
+        return (self.dim, self.expression, self.global_scaling_const,
                 self.is_complex_valued)
 
     def __repr__(self):
-        if self._dim is not None:
-            return "ExprKnl%dD" % self.dim
-        else:
-            return "ExprKnl"
-
-    def get_expression(self, dist_vec):
-        if self.expression is None:
-            raise RuntimeError("expression in ExpressionKernel has not "
-                    "been determined yet (this could be due to a PDE kernel "
-                    "not having learned its dimensionality yet)")
+        return "ExprKnl%dD" % self.dim
 
+    def get_expression(self, scaled_dist_vec):
         from sumpy.symbolic import PymbolicToSympyMapperWithSymbols
         expr = PymbolicToSympyMapperWithSymbols()(self.expression)
 
-        if self.dim != len(dist_vec):
+        if self.dim != len(scaled_dist_vec):
             raise ValueError("dist_vec length does not match expected dimension")
 
         from sumpy.symbolic import Symbol
         expr = expr.subs(dict(
             (Symbol("d%d" % i), dist_vec_i)
-            for i, dist_vec_i in enumerate(dist_vec)
+            for i, dist_vec_i in enumerate(scaled_dist_vec)
             ))
 
         return expr
 
-    def get_scaling(self):
+    def get_global_scaling_const(self):
         """Return a global scaling of the kernel."""
 
-        if self.scaling is None:
-            raise RuntimeError("scaling in ExpressionKernel has not "
-                    "been determined yet (this could be due to a PDE kernel "
-                    "not having learned its dimensionality yet)")
-
         from sumpy.symbolic import PymbolicToSympyMapperWithSymbols
-        return PymbolicToSympyMapperWithSymbols()(self.scaling)
+        return PymbolicToSympyMapperWithSymbols()(
+                self.global_scaling_const)
 
     def update_persistent_hash(self, key_hash, key_builder):
         key_hash.update(type(self).__name__.encode("utf8"))
@@ -299,12 +370,12 @@ class ExpressionKernel(Kernel):
 one_kernel_2d = ExpressionKernel(
         dim=2,
         expression=1,
-        scaling=1,
+        global_scaling_const=1,
         is_complex_valued=False)
 one_kernel_3d = ExpressionKernel(
         dim=3,
         expression=1,
-        scaling=1,
+        global_scaling_const=1,
         is_complex_valued=False)
 
 
@@ -324,23 +395,35 @@ class LaplaceKernel(ExpressionKernel):
             expr = 1/r
             scaling = 1/(4*var("pi"))
         else:
-            raise RuntimeError("unsupported dimensionality")
+            raise NotImplementedError("unsupported dimensionality")
 
-        ExpressionKernel.__init__(
-                self,
+        super(LaplaceKernel, self).__init__(
                 dim,
                 expression=expr,
-                scaling=scaling,
+                global_scaling_const=scaling,
                 is_complex_valued=False)
 
+    has_efficient_scale_adjustment = True
+
+    def adjust_for_kernel_scaling(self, expr, rscale, nderivatives):
+        if self.dim == 2:
+            if nderivatives == 0:
+                import sumpy.symbolic as sp
+                return (expr + sp.log(rscale))
+            else:
+                return expr
+
+        elif self.dim == 3:
+            return expr / rscale
+
+        else:
+            raise NotImplementedError("unsupported dimensionality")
+
     def __getinitargs__(self):
-        return (self._dim,)
+        return (self.dim,)
 
     def __repr__(self):
-        if self._dim is not None:
-            return "LapKnl%dD" % self.dim
-        else:
-            return "LapKnl"
+        return "LapKnl%dD" % self.dim
 
     mapper_method = "map_laplace_kernel"
 
@@ -359,21 +442,17 @@ class BiharmonicKernel(ExpressionKernel):
         else:
             raise RuntimeError("unsupported dimensionality")
 
-        ExpressionKernel.__init__(
-                self,
+        super(BiharmonicKernel, self).__init__(
                 dim,
                 expression=expr,
-                scaling=scaling,
+                global_scaling_const=scaling,
                 is_complex_valued=False)
 
     def __getinitargs__(self):
-        return (self._dim,)
+        return (self.dim,)
 
     def __repr__(self):
-        if self._dim is not None:
-            return "BiharmKnl%dD" % self.dim
-        else:
-            return "BiharmKnl"
+        return "BiharmKnl%dD" % self.dim
 
     mapper_method = "map_biharmonic_kernel"
 
@@ -403,18 +482,17 @@ class HelmholtzKernel(ExpressionKernel):
         else:
             raise RuntimeError("unsupported dimensionality")
 
-        ExpressionKernel.__init__(
-                self,
+        super(HelmholtzKernel, self).__init__(
                 dim,
                 expression=expr,
-                scaling=scaling,
+                global_scaling_const=scaling,
                 is_complex_valued=True)
 
         self.helmholtz_k_name = helmholtz_k_name
         self.allow_evanescent = allow_evanescent
 
     def __getinitargs__(self):
-        return (self._dim, self.helmholtz_k_name,
+        return (self.dim, self.helmholtz_k_name,
                 self.allow_evanescent)
 
     def update_persistent_hash(self, key_hash, key_builder):
@@ -470,17 +548,16 @@ class YukawaKernel(ExpressionKernel):
         else:
             raise RuntimeError("unsupported dimensionality")
 
-        ExpressionKernel.__init__(
-                self,
+        super(YukawaKernel, self).__init__(
                 dim,
                 expression=expr,
-                scaling=scaling,
+                global_scaling_const=scaling,
                 is_complex_valued=True)
 
         self.yukawa_lambda_name = yukawa_lambda_name
 
     def __getinitargs__(self):
-        return (self._dim, self.yukawa_lambda_name)
+        return (self.dim, self.yukawa_lambda_name)
 
     def update_persistent_hash(self, key_hash, key_builder):
         key_hash.update(type(self).__name__.encode("utf8"))
@@ -549,15 +626,14 @@ class StokesletKernel(ExpressionKernel):
         self.icomp = icomp
         self.jcomp = jcomp
 
-        ExpressionKernel.__init__(
-                self,
+        super(StokesletKernel, self).__init__(
                 dim,
                 expression=expr,
-                scaling=scaling,
+                global_scaling_const=scaling,
                 is_complex_valued=False)
 
     def __getinitargs__(self):
-        return (self._dim, self.icomp, self.jcomp, self.viscosity_mu_name)
+        return (self.dim, self.icomp, self.jcomp, self.viscosity_mu_name)
 
     def update_persistent_hash(self, key_hash, key_builder):
         key_hash.update(type(self).__name__.encode())
@@ -615,15 +691,14 @@ class StressletKernel(ExpressionKernel):
         self.jcomp = jcomp
         self.kcomp = kcomp
 
-        ExpressionKernel.__init__(
-                self,
+        super(StressletKernel, self).__init__(
                 dim,
                 expression=expr,
-                scaling=scaling,
+                global_scaling_const=scaling,
                 is_complex_valued=False)
 
     def __getinitargs__(self):
-        return (self._dim, self.icomp, self.jcomp, self.kcomp,
+        return (self.dim, self.icomp, self.jcomp, self.kcomp,
                       self.viscosity_mu_name)
 
     def update_persistent_hash(self, key_hash, key_builder):
@@ -665,8 +740,16 @@ class KernelWrapper(Kernel):
     def is_complex_valued(self):
         return self.inner_kernel.is_complex_valued
 
-    def get_expression(self, dist_vec):
-        return self.inner_kernel.get_expression(dist_vec)
+    def get_expression(self, scaled_dist_vec):
+        return self.inner_kernel.get_expression(scaled_dist_vec)
+
+    @property
+    def has_efficient_scale_adjustment(self):
+        return self.inner_kernel.has_efficient_scale_adjustment
+
+    def adjust_for_kernel_scaling(self, expr, rscale, nderivatives):
+        return self.inner_kernel.adjust_for_kernel_scaling(
+                expr, rscale, nderivatives)
 
     def postprocess_at_source(self, expr, avec):
         return self.inner_kernel.postprocess_at_source(expr, avec)
@@ -674,8 +757,8 @@ class KernelWrapper(Kernel):
     def postprocess_at_target(self, expr, avec):
         return self.inner_kernel.postprocess_at_target(expr, avec)
 
-    def get_scaling(self):
-        return self.inner_kernel.get_scaling()
+    def get_global_scaling_const(self):
+        return self.inner_kernel.get_global_scaling_const()
 
     def get_code_transformer(self):
         return self.inner_kernel.get_code_transformer()
diff --git a/sumpy/p2e.py b/sumpy/p2e.py
index b8c60a6423bc1dd704cf63cf915831c084480b7b..d9354da71aeb53b606fe3dc711c8da1262f8ada7 100644
--- a/sumpy/p2e.py
+++ b/sumpy/p2e.py
@@ -77,13 +77,16 @@ class P2EBase(KernelCacheWrapper):
         from sumpy.symbolic import make_sym_vector
         avec = make_sym_vector("a", self.dim)
 
+        import sumpy.symbolic as sp
+        rscale = sp.Symbol("rscale")
+
         from sumpy.assignment_collection import SymbolicAssignmentCollection
         sac = SymbolicAssignmentCollection()
 
         coeff_names = [
                 sac.assign_unique("coeff%d" % i, coeff_i)
                 for i, coeff_i in enumerate(
-                    self.expansion.coefficients_from_source(avec, None))]
+                    self.expansion.coefficients_from_source(avec, None, rscale))]
 
         sac.run_global_cse()
 
@@ -97,7 +100,7 @@ class P2EBase(KernelCacheWrapper):
                 )
 
     def get_cache_key(self):
-        return (type(self).__name__, self.expansion)
+        return (type(self).__name__, self.name, self.expansion)
 
 # }}}
 
@@ -144,6 +147,7 @@ class P2EFromSingleBox(P2EBase):
                     lp.GlobalArg("box_source_starts,box_source_counts_nonchild",
                         None, shape=None),
                     lp.GlobalArg("centers", None, shape="dim, aligned_nboxes"),
+                    lp.ValueArg("rscale", None),
                     lp.GlobalArg("tgt_expansions", None,
                         shape=("nboxes", ncoeffs), offset=lp.auto),
                     lp.ValueArg("nboxes,aligned_nboxes,tgt_base_ibox", np.int32),
@@ -178,10 +182,16 @@ class P2EFromSingleBox(P2EBase):
         :arg centers:
         :arg sources:
         :arg strengths:
+        :arg rscale:
         """
+        centers = kwargs.pop("centers")
+        # "1" may be passed for rscale, which won't have its type
+        # meaningfully inferred. Make the type of rscale explicit.
+        rscale = centers.dtype.type(kwargs.pop("rscale"))
+
         knl = self.get_cached_optimized_kernel()
 
-        return knl(queue, **kwargs)
+        return knl(queue, centers=centers, rscale=rscale, **kwargs)
 
 # }}}
 
@@ -277,10 +287,16 @@ class P2EFromCSR(P2EBase):
         :arg centers:
         :arg sources:
         :arg strengths:
+        :arg rscale:
         """
         knl = self.get_cached_optimized_kernel()
 
-        return knl(queue, **kwargs)
+        centers = kwargs.pop("centers")
+        # "1" may be passed for rscale, which won't have its type
+        # meaningfully inferred. Make the type of rscale explicit.
+        rscale = centers.dtype.type(kwargs.pop("rscale"))
+
+        return knl(queue, centers=centers, rscale=rscale, **kwargs)
 
 # }}}
 
diff --git a/sumpy/p2p.py b/sumpy/p2p.py
index 9e07fbe01adffbe9a0732c46a80e87ee75b86e98..fb0a1cf573231d112118a18fd1b1a1e0d3337013 100644
--- a/sumpy/p2p.py
+++ b/sumpy/p2p.py
@@ -79,8 +79,10 @@ class P2PBase(KernelComputation, KernelCacheWrapper):
                 sac.assign_unique("knl%d" % i,
                     knl.postprocess_at_target(
                         knl.postprocess_at_source(
-                            knl.get_expression(dvec), dvec),
-                        dvec))
+                            knl.get_expression(dvec),
+                            dvec),
+                        dvec)
+                    )
                 for i, knl in enumerate(self.kernels)]
 
         sac.run_global_cse()
diff --git a/sumpy/point_calculus.py b/sumpy/point_calculus.py
index 2706aa98f403ec3f00480469ab97dd2b399aaba4..e938a9d8ae34d5a510a72914c6db95376c627bdf 100644
--- a/sumpy/point_calculus.py
+++ b/sumpy/point_calculus.py
@@ -46,7 +46,7 @@ class CalculusPatch(object):
     .. automethod:: dy
     .. automethod:: dy
     .. automethod:: laplace
-    .. automethod:: eval_at_0
+    .. automethod:: eval_at_center
     .. autoattribute:: x
     .. autoattribute:: y
     .. autoattribute:: z
diff --git a/sumpy/qbx.py b/sumpy/qbx.py
index 8892739e3b9050d3d617e8f210b288b3da4b7e28..b5806ef2ddbc63c50f99c443c1ecda2a011ae84f 100644
--- a/sumpy/qbx.py
+++ b/sumpy/qbx.py
@@ -61,7 +61,9 @@ def stringify_expn_index(i):
 
 
 def expand(expansion_nr, sac, expansion, avec, bvec):
-    coefficients = expansion.coefficients_from_source(avec, bvec)
+    rscale = sym.Symbol("rscale")
+
+    coefficients = expansion.coefficients_from_source(avec, bvec, rscale)
 
     assigned_coeffs = [
             sym.Symbol(
@@ -71,7 +73,7 @@ def expand(expansion_nr, sac, expansion, avec, bvec):
             for i in expansion.get_coefficient_identifiers()]
 
     return sac.assign_unique("expn%d_result" % expansion_nr,
-            expansion.evaluate(assigned_coeffs, bvec))
+            expansion.evaluate(assigned_coeffs, bvec, rscale))
 
 
 # {{{ layer potential computation
@@ -101,6 +103,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper):
         return """
             <> a[idim] = center[idim,itgt] - src[idim,isrc] {dup=idim}
             <> b[idim] = tgt[idim,itgt] - center[idim,itgt] {dup=idim}
+            <> rscale = expansion_radii[itgt]
             """
 
     def get_src_tgt_arguments(self):
@@ -111,6 +114,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper):
                     shape=(self.dim, "ntargets"), order="C"),
                 lp.GlobalArg("center", None,
                     shape=(self.dim, "ntargets"), order="C"),
+                lp.GlobalArg("expansion_radii", None, shape="ntargets"),
                 lp.ValueArg("nsources", None),
                 lp.ValueArg("ntargets", None),
                 ]
@@ -242,7 +246,8 @@ class LayerPotential(LayerPotentialBase):
                 for iknl in range(len(self.expansions))
                 ]
 
-    def __call__(self, queue, targets, sources, centers, strengths, **kwargs):
+    def __call__(self, queue, targets, sources, centers, strengths, expansion_radii,
+            **kwargs):
         """
         :arg strengths: are required to have area elements and quadrature weights
             already multiplied in.
@@ -253,7 +258,8 @@ class LayerPotential(LayerPotentialBase):
         for i, dens in enumerate(strengths):
             kwargs["strength_%d" % i] = dens
 
-        return knl(queue, src=sources, tgt=targets, center=centers, **kwargs)
+        return knl(queue, src=sources, tgt=targets, center=centers,
+                expansion_radii=expansion_radii, **kwargs)
 
 # }}}
 
@@ -283,11 +289,11 @@ class LayerPotentialMatrixGenerator(LayerPotentialBase):
                 for iknl in range(len(self.expansions))
                 ]
 
-    def __call__(self, queue, targets, sources, centers, **kwargs):
+    def __call__(self, queue, targets, sources, centers, expansion_radii, **kwargs):
         knl = self.get_optimized_kernel()
 
         return knl(queue, src=sources, tgt=targets, center=centers,
-                **kwargs)
+                expansion_radii=expansion_radii, **kwargs)
 
 # }}}
 
diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py
index 93295bca64b2e5af09b793267b8825a6019f58a6..55f91284d62c536c290c2ccf8f2c07f4d829796f 100644
--- a/sumpy/symbolic.py
+++ b/sumpy/symbolic.py
@@ -25,7 +25,6 @@ THE SOFTWARE.
 
 import six
 from six.moves import range
-from six.moves import zip
 
 import numpy as np
 from pymbolic.mapper import IdentityMapper as IdentityMapperBase
@@ -77,7 +76,7 @@ _find_symbolic_backend()
 # Symbolic API common to SymEngine and sympy.
 # Before adding a function here, make sure it's present in both modules.
 SYMBOLIC_API = """
-Add Basic Mul Pow exp sqrt symbols sympify cos sin atan2 Function Symbol
+Add Basic Mul Pow exp sqrt log symbols sympify cos sin atan2 Function Symbol
 Derivative Integer Matrix Subs I pi functions""".split()
 
 if USE_SYMENGINE:
@@ -193,12 +192,14 @@ def make_sym_vector(name, components):
             [sym.Symbol("%s%d" % (name, i)) for i in range(components)])
 
 
-def vector_subs(expr, old, new):
-    result = expr
-    for old_i, new_i in zip(old, new):
-        result = result.subs(old_i, new_i)
+def vector_subs(expr, from_vec, to_vec):
+    subs_pairs = []
+    assert (from_vec.rows, from_vec.cols) == (to_vec.rows, to_vec.cols)
+    for irow in range(from_vec.rows):
+        for icol in range(from_vec.cols):
+            subs_pairs.append((from_vec[irow, icol], to_vec[irow, icol]))
 
-    return result
+    return expr.subs(subs_pairs)
 
 
 def find_power_of(base, prod):
diff --git a/sumpy/tools.py b/sumpy/tools.py
index 2aa4732c02b6181d1324b9aa1dedee209cf949f9..65cd711e794f5dd103d03f9f177374366cdb233d 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -95,14 +95,21 @@ class MiDerivativeTaker(object):
 
 class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker):
     """
-    See :class:`sumpy.expansion.LinearRecurrenceBasedDerivativeWrangler`
+    The derivative taker for expansions that use
+    :class:`sumpy.expansion.LinearRecurrenceBasedDerivativeWrangler`
     """
 
     def __init__(self, expr, var_list, wrangler):
-        MiDerivativeTaker.__init__(self, expr, var_list)
+        super(LinearRecurrenceBasedMiDerivativeTaker, self).__init__(
+                expr, var_list)
         self.wrangler = wrangler
 
+    @memoize_method
     def diff(self, mi):
+        """
+        :arg mi: a multi-index (tuple) indicating how many x/y derivatives are
+            to be taken.
+        """
         try:
             expr = self.cache_by_mi[mi]
         except KeyError:
@@ -119,7 +126,7 @@ class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker):
 
                 recurrence = (
                         self.wrangler.try_get_recurrence_for_derivative(
-                            next_mi, self.cache_by_mi))
+                            next_mi, self.cache_by_mi, rscale=1))
 
                 if recurrence is not None:
                     expr = Add(*tuple(
@@ -129,6 +136,7 @@ class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker):
                     expr = expr.diff(next_deriv)
 
                 self.cache_by_mi[next_mi] = expr
+
         return expr
 
 # }}}
@@ -271,7 +279,7 @@ class KernelComputation(object):
         return [
                 lp.Assignment(id=None,
                     assignee="knl_%d_scaling" % i,
-                    expression=sympy_conv(kernel.get_scaling()),
+                    expression=sympy_conv(kernel.get_global_scaling_const()),
                     temp_var_type=dtype)
                 for i, (kernel, dtype) in enumerate(
                     zip(self.kernels, self.value_dtypes))]
diff --git a/sumpy/toys.py b/sumpy/toys.py
index b32613d37e64103f0db0cc024aa351e3db6b39b0..c4d438ba7354a29be6119699849d78b83a4ca459 100644
--- a/sumpy/toys.py
+++ b/sumpy/toys.py
@@ -27,6 +27,7 @@ from six.moves import range  # noqa: F401
 
 from pytools import memoize_method
 from numbers import Number
+from sumpy.kernel import TargetDerivativeRemover
 
 import numpy as np  # noqa: F401
 import loopy as lp  # noqa: F401
@@ -43,11 +44,14 @@ class ToyContext(object):
             mpole_expn_class=None,
             local_expn_class=None,
             expansion_factory=None,
-            extra_source_kwargs=None):
+            extra_source_kwargs=None,
+            extra_kernel_kwargs=None):
         self.cl_context = cl_context
         self.queue = cl.CommandQueue(self.cl_context)
         self.kernel = kernel
 
+        self.no_target_deriv_kernel = TargetDerivativeRemover()(kernel)
+
         if expansion_factory is None:
             from sumpy.expansion import DefaultExpansionFactory
             expansion_factory = DefaultExpansionFactory()
@@ -58,12 +62,20 @@ class ToyContext(object):
             local_expn_class = \
                     expansion_factory.get_local_expansion_class(kernel)
 
+        self.mpole_expn_class = mpole_expn_class
+        self.local_expn_class = local_expn_class
+
         if extra_source_kwargs is None:
             extra_source_kwargs = {}
+        if extra_kernel_kwargs is None:
+            extra_kernel_kwargs = {}
 
-        self.mpole_expn_class = mpole_expn_class
-        self.local_expn_class = local_expn_class
         self.extra_source_kwargs = extra_source_kwargs
+        self.extra_kernel_kwargs = extra_kernel_kwargs
+
+        extra_source_and_kernel_kwargs = extra_source_kwargs.copy()
+        extra_source_and_kernel_kwargs.update(extra_kernel_kwargs)
+        self.extra_source_and_kernel_kwargs = extra_source_and_kernel_kwargs
 
     @memoize_method
     def get_p2p(self):
@@ -74,57 +86,57 @@ class ToyContext(object):
     def get_p2m(self, order):
         from sumpy import P2EFromSingleBox
         return P2EFromSingleBox(self.cl_context,
-                self.mpole_expn_class(self.kernel, order),
+                self.mpole_expn_class(self.no_target_deriv_kernel, order),
                 [self.kernel])
 
     @memoize_method
     def get_p2l(self, order):
         from sumpy import P2EFromSingleBox
         return P2EFromSingleBox(self.cl_context,
-                self.local_expn_class(self.kernel, order),
+                self.local_expn_class(self.no_target_deriv_kernel, order),
                 [self.kernel])
 
     @memoize_method
     def get_m2p(self, order):
         from sumpy import E2PFromSingleBox
         return E2PFromSingleBox(self.cl_context,
-                self.mpole_expn_class(self.kernel, order),
+                self.mpole_expn_class(self.no_target_deriv_kernel, order),
                 [self.kernel])
 
     @memoize_method
     def get_l2p(self, order):
         from sumpy import E2PFromSingleBox
         return E2PFromSingleBox(self.cl_context,
-                self.local_expn_class(self.kernel, order),
+                self.local_expn_class(self.no_target_deriv_kernel, order),
                 [self.kernel])
 
     @memoize_method
     def get_m2m(self, from_order, to_order):
         from sumpy import E2EFromCSR
         return E2EFromCSR(self.cl_context,
-                self.mpole_expn_class(self.kernel, from_order),
-                self.mpole_expn_class(self.kernel, to_order))
+                self.mpole_expn_class(self.no_target_deriv_kernel, from_order),
+                self.mpole_expn_class(self.no_target_deriv_kernel, to_order))
 
     @memoize_method
     def get_m2l(self, from_order, to_order):
         from sumpy import E2EFromCSR
         return E2EFromCSR(self.cl_context,
-                self.mpole_expn_class(self.kernel, from_order),
-                self.local_expn_class(self.kernel, to_order))
+                self.mpole_expn_class(self.no_target_deriv_kernel, from_order),
+                self.local_expn_class(self.no_target_deriv_kernel, to_order))
 
     @memoize_method
     def get_l2l(self, from_order, to_order):
         from sumpy import E2EFromCSR
         return E2EFromCSR(self.cl_context,
-                self.local_expn_class(self.kernel, from_order),
-                self.local_expn_class(self.kernel, to_order))
+                self.local_expn_class(self.no_target_deriv_kernel, from_order),
+                self.local_expn_class(self.no_target_deriv_kernel, to_order))
 
 # }}}
 
 
 # {{{ helpers
 
-def _p2e(psource, center, order, p2e, expn_class, expn_kwargs):
+def _p2e(psource, center, rscale, order, p2e, expn_class, expn_kwargs):
     source_boxes = np.array([0], dtype=np.int32)
     box_source_starts = np.array([0], dtype=np.int32)
     box_source_counts_nonchild = np.array(
@@ -142,14 +154,16 @@ def _p2e(psource, center, order, p2e, expn_class, expn_kwargs):
             centers=centers,
             sources=psource.points,
             strengths=psource.weights,
+            rscale=rscale,
             nboxes=1,
             tgt_base_ibox=0,
 
             #flags="print_hl_cl",
-            out_host=True, **toy_ctx.extra_source_kwargs)
+            out_host=True,
+            **toy_ctx.extra_source_and_kernel_kwargs)
 
-    return expn_class(toy_ctx, center, order, coeffs[0], derived_from=psource,
-            **expn_kwargs)
+    return expn_class(toy_ctx, center, rscale, order, coeffs[0],
+            derived_from=psource, **expn_kwargs)
 
 
 def _e2p(psource, targets, e2p):
@@ -173,14 +187,15 @@ def _e2p(psource, targets, e2p):
             box_target_starts=box_target_starts,
             box_target_counts_nonchild=box_target_counts_nonchild,
             centers=centers,
+            rscale=psource.rscale,
             targets=targets,
             #flags="print_hl_cl",
-            out_host=True, **toy_ctx.extra_source_kwargs)
+            out_host=True, **toy_ctx.extra_kernel_kwargs)
 
     return pot
 
 
-def _e2e(psource, to_center, to_order, e2e, expn_class, expn_kwargs):
+def _e2e(psource, to_center, to_rscale, to_order, e2e, expn_class, expn_kwargs):
     toy_ctx = psource.toy_ctx
 
     target_boxes = np.array([1], dtype=np.int32)
@@ -211,10 +226,14 @@ def _e2e(psource, to_center, to_order, e2e, expn_class, expn_kwargs):
             src_box_starts=src_box_starts,
             src_box_lists=src_box_lists,
             centers=centers,
+
+            src_rscale=psource.rscale,
+            tgt_rscale=to_rscale,
+
             #flags="print_hl_cl",
-            out_host=True, **toy_ctx.extra_source_kwargs)
+            out_host=True, **toy_ctx.extra_kernel_kwargs)
 
-    return expn_class(toy_ctx, to_center, to_order, to_coeffs[1],
+    return expn_class(toy_ctx, to_center, to_rscale, to_order, to_coeffs[1],
             derived_from=psource, **expn_kwargs)
 
 # }}}
@@ -298,7 +317,8 @@ class PointSources(PotentialSource):
     def eval(self, targets):
         evt, (potential,) = self.toy_ctx.get_p2p()(
                 self.toy_ctx.queue, targets, self.points, [self.weights],
-                out_host=True, **self.toy_ctx.extra_source_kwargs)
+                out_host=True,
+                **self.toy_ctx.extra_source_and_kernel_kwargs)
 
         return potential
 
@@ -321,10 +341,11 @@ class ExpansionPotentialSource(PotentialSource):
        Passed to matplotlib.text(). Used for customizing the expansion
        label. Just for visualization, purely advisory.
     """
-    def __init__(self, toy_ctx, center, order, coeffs, derived_from,
+    def __init__(self, toy_ctx, center, rscale, order, coeffs, derived_from,
             radius=None, expn_style=None, text_kwargs=None):
         super(ExpansionPotentialSource, self).__init__(toy_ctx)
         self.center = np.asarray(center)
+        self.rscale = rscale
         self.order = order
         self.coeffs = coeffs
 
@@ -383,19 +404,19 @@ class Product(PotentialExpressionNode):
 # }}}
 
 
-def multipole_expand(psource, center, order=None, **expn_kwargs):
+def multipole_expand(psource, center, order=None, rscale=1, **expn_kwargs):
     if isinstance(psource, PointSources):
         if order is None:
             raise ValueError("order may not be None")
 
-        return _p2e(psource, center, order, psource.toy_ctx.get_p2m(order),
+        return _p2e(psource, center, rscale, order, psource.toy_ctx.get_p2m(order),
                 MultipoleExpansion, expn_kwargs)
 
     elif isinstance(psource, MultipoleExpansion):
         if order is None:
             order = psource.order
 
-        return _e2e(psource, center, order,
+        return _e2e(psource, center, rscale, order,
                 psource.toy_ctx.get_m2m(psource.order, order),
                 MultipoleExpansion, expn_kwargs)
 
@@ -404,19 +425,19 @@ def multipole_expand(psource, center, order=None, **expn_kwargs):
                 % type(psource).__name__)
 
 
-def local_expand(psource, center, order=None, **expn_kwargs):
+def local_expand(psource, center, order=None, rscale=1, **expn_kwargs):
     if isinstance(psource, PointSources):
         if order is None:
             raise ValueError("order may not be None")
 
-        return _p2e(psource, center, order, psource.toy_ctx.get_p2l(order),
+        return _p2e(psource, center, rscale, order, psource.toy_ctx.get_p2l(order),
                 LocalExpansion, expn_kwargs)
 
     elif isinstance(psource, MultipoleExpansion):
         if order is None:
             order = psource.order
 
-        return _e2e(psource, center, order,
+        return _e2e(psource, center, rscale, order,
                 psource.toy_ctx.get_m2l(psource.order, order),
                 LocalExpansion, expn_kwargs)
 
@@ -424,7 +445,7 @@ def local_expand(psource, center, order=None, **expn_kwargs):
         if order is None:
             order = psource.order
 
-        return _e2e(psource, center, order,
+        return _e2e(psource, center, rscale, order,
                 psource.toy_ctx.get_l2l(psource.order, order),
                 LocalExpansion, expn_kwargs)
 
diff --git a/sumpy/version.py b/sumpy/version.py
index 6333fc0cf50d9373f1beb568f773eba9c6422f25..97d6e94e2ed589c3909cf397fb480ed28b834da6 100644
--- a/sumpy/version.py
+++ b/sumpy/version.py
@@ -25,4 +25,4 @@ VERSION = (2016, 1)
 VERSION_STATUS = "beta1"
 VERSION_TEXT = ".".join(str(x) for x in VERSION) + VERSION_STATUS
 
-KERNEL_VERSION = 15
+KERNEL_VERSION = 18
diff --git a/test/test_codegen.py b/test/test_codegen.py
index fb1655986bb2ffa6b69cd6c1bd86bd5de594e342..0a02b6b4c39d9223380555bb705b9227ed1caf17 100644
--- a/test/test_codegen.py
+++ b/test/test_codegen.py
@@ -91,7 +91,7 @@ def test_line_taylor_coeff_growth():
     expn = LineTaylorLocalExpansion(LaplaceKernel(2), order)
     avec = make_sym_vector("a", 2)
     bvec = make_sym_vector("b", 2)
-    coeffs = expn.coefficients_from_source(avec, bvec)
+    coeffs = expn.coefficients_from_source(avec, bvec, rscale=1)
 
     sym2pymbolic = SympyToPymbolicMapper()
     coeffs_pymbolic = [sym2pymbolic(c) for c in coeffs]
diff --git a/test/test_kernels.py b/test/test_kernels.py
index c32435d3d30cbf49e924dd7de4a89a7d0784fcee..65295e29451e039619ef65076926cc1b5b182d22 100644
--- a/test/test_kernels.py
+++ b/test/test_kernels.py
@@ -214,6 +214,8 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative):
 
         targets = fp.points
 
+        rscale = 0.5  # pick something non-1
+
         # {{{ apply p2e
 
         evt, (mpoles,) = p2e(queue,
@@ -225,6 +227,7 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative):
                 strengths=strengths,
                 nboxes=1,
                 tgt_base_ibox=0,
+                rscale=rscale,
 
                 #flags="print_hl_cl",
                 out_host=True, **extra_source_kwargs)
@@ -247,6 +250,8 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative):
                 box_target_counts_nonchild=box_target_counts_nonchild,
                 centers=centers,
                 targets=targets,
+                rscale=rscale,
+
                 #flags="print_hl_cl",
                 out_host=True, **extra_kwargs)
 
@@ -402,7 +407,7 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
         orders = [2, 3, 4]
     nboxes = centers.shape[-1]
 
-    def eval_at(e2p, source_box_nr):
+    def eval_at(e2p, source_box_nr, rscale):
         e2p_target_boxes = np.array([source_box_nr], dtype=np.int32)
 
         # These are indexed by global box numbers.
@@ -422,6 +427,9 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
                 box_target_counts_nonchild=e2p_box_target_counts_nonchild,
                 centers=centers,
                 targets=targets,
+
+                rscale=rscale,
+
                 out_host=True, **extra_kwargs
                 )
 
@@ -452,6 +460,11 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
 
         # }}}
 
+        m1_rscale = 0.5
+        m2_rscale = 0.25
+        l1_rscale = 0.5
+        l2_rscale = 0.25
+
         # {{{ apply P2M
 
         p2m_source_boxes = np.array([0], dtype=np.int32)
@@ -469,6 +482,7 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
                 sources=sources,
                 strengths=strengths,
                 nboxes=nboxes,
+                rscale=m1_rscale,
 
                 tgt_base_ibox=0,
 
@@ -479,7 +493,7 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
 
         ntargets = targets.shape[-1]
 
-        pot = eval_at(m2p, 0)
+        pot = eval_at(m2p, 0, m1_rscale)
 
         err = la.norm((pot - pot_direct)/res**2)
         err = err / (la.norm(pot_direct) / res**2)
@@ -503,12 +517,16 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
                 src_box_starts=m2m_src_box_starts,
                 src_box_lists=m2m_src_box_lists,
                 centers=centers,
+
+                src_rscale=m1_rscale,
+                tgt_rscale=m2_rscale,
+
                 #flags="print_hl_cl",
                 out_host=True, **extra_kwargs)
 
         # }}}
 
-        pot = eval_at(m2p, 1)
+        pot = eval_at(m2p, 1, m2_rscale)
 
         err = la.norm((pot - pot_direct)/res**2)
         err = err / (la.norm(pot_direct) / res**2)
@@ -531,12 +549,16 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
                 src_box_starts=m2l_src_box_starts,
                 src_box_lists=m2l_src_box_lists,
                 centers=centers,
+
+                src_rscale=m2_rscale,
+                tgt_rscale=l1_rscale,
+
                 #flags="print_hl_cl",
                 out_host=True, **extra_kwargs)
 
         # }}}
 
-        pot = eval_at(l2p, 2)
+        pot = eval_at(l2p, 2, l1_rscale)
 
         err = la.norm((pot - pot_direct)/res**2)
         err = err / (la.norm(pot_direct) / res**2)
@@ -559,12 +581,16 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
                 src_box_starts=l2l_src_box_starts,
                 src_box_lists=l2l_src_box_lists,
                 centers=centers,
+
+                src_rscale=l1_rscale,
+                tgt_rscale=l2_rscale,
+
                 #flags="print_hl_wrapper",
                 out_host=True, **extra_kwargs)
 
         # }}}
 
-        pot = eval_at(l2p, 3)
+        pot = eval_at(l2p, 3, l2_rscale)
 
         err = la.norm((pot - pot_direct)/res**2)
         err = err / (la.norm(pot_direct) / res**2)
diff --git a/test/test_qbx.py b/test/test_qbx.py
index aa86edf876fd9dc9aef4e26c2ce10992cfa73240..3758a5d7b5af4f031dabc9770da884860c585382 100644
--- a/test/test_qbx.py
+++ b/test/test_qbx.py
@@ -82,8 +82,11 @@ def test_direct(ctx_getter):
         radius = 7 * h
         centers = unit_circle * (1 - radius)
 
+        expansion_radii = np.ones(n) * radius
+
         strengths = (sigma * h,)
-        evt, (result_qbx,) = lpot(queue, targets, sources, centers, strengths)
+        evt, (result_qbx,) = lpot(queue, targets, sources, centers, strengths,
+                expansion_radii=expansion_radii)
 
         eocrec.add_data_point(h, np.max(np.abs(result_ref - result_qbx)))