diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index f1bc491b3a737d3c867917abf2dccaf516ba28b1..6477b9ab4d9b7de6ba5b3a096ee57de42e29b51d 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -58,6 +58,9 @@ class ExpansionBase(object):
     def get_args(self):
         return self.kernel.get_args()
 
+    def get_source_args(self):
+        return self.kernel.get_source_args()
+
     # }}}
 
     def coefficients_from_source(self, expr, avec, bvec):
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index 84a19e9a78c51b0d4a99a85992056af6bf188c68..924055dc1d8fc670ff45a55eb4c2884b1e9a47b0 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -49,8 +49,23 @@ class VolumeTaylorMultipoleExpansion(
             from sumpy.symbolic import make_sym_vector
             dir_vec = make_sym_vector(kernel.dir_vec_name, kernel.dim)
 
-            result = [0]
-            1/0
+            coeff_identifiers = self.get_coefficient_identifiers()
+            result = [0] * len(coeff_identifiers)
+
+            for idim in xrange(kernel.dim):
+                for i, mi in enumerate(coeff_identifiers):
+                    if mi[idim] == 0:
+                        continue
+
+                    derivative_mi = tuple(mi_i - 1 if iaxis == idim else mi_i
+                            for iaxis, mi_i in enumerate(mi))
+
+                    result[i] += (
+                        - 1/mi_factorial(derivative_mi)
+                        * mi_power(avec, derivative_mi)
+                        * dir_vec[idim])
+
+            return result
         else:
             return [
                     mi_power(avec, mi) / mi_factorial(mi)
diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index d9133300f8c55cea0b50d4d9777d0c888bad4e36..0f9a9a474251584b04885d76e2ac34e56550656a 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -127,6 +127,13 @@ class Kernel(object):
         """
         return []
 
+    def get_source_args(self):
+        """Return list of :cls:`loopy.Argument` instances describing
+        extra arguments used by kernel in picking up contributions
+        from point sources.
+        """
+        return []
+
     def __sub__(self, other):
         return DifferenceKernel(self, other)
 
@@ -252,6 +259,9 @@ class DifferenceKernel(Kernel):
 
         Kernel.__init__(self, self.kernel_plus.dim)
 
+        # FIXME
+        raise NotImplementedError("difference kernel mostly unimplemented")
+
     def __getinitargs__(self):
         return (self.kernel_plus, self.kernel_minus)
 
@@ -265,8 +275,6 @@ class DifferenceKernel(Kernel):
 
     mapper_method = "map_difference_kernel"
 
-    # FIXME mostly unimplemented
-
 # }}}
 
 
@@ -305,6 +313,9 @@ class KernelWrapper(Kernel):
     def get_args(self):
         return self.inner_kernel.get_args()
 
+    def get_source_args(self):
+        return self.inner_kernel.get_source_args()
+
 # }}}
 
 
@@ -376,10 +387,10 @@ class DirectionalDerivative(DerivativeBase):
         return r"%s . \/_%s %s" % (
                 self.dir_vec_name, self.directional_kind[0], self.inner_kernel)
 
-    def get_args(self):
-        return self.inner_kernel.get_args() + [
-            lp.GlobalArg(self.dir_vec_name, None, shape=lp.auto,
-                dim_tags="sep,C")]
+    def get_source_args(self):
+        return [
+            lp.GlobalArg(self.dir_vec_name, None, shape=(self.dim, "nsrc"),
+                dim_tags="sep,C")] + self.inner_kernel.get_source_args()
 
 
 class DirectionalTargetDerivative(DirectionalDerivative):
diff --git a/sumpy/p2e.py b/sumpy/p2e.py
index 8114845ab4d78a781fbbe0f0d55a39170ed7da29..cd828d8d5d9a094edf4bfe0b721bc6eca4590785 100644
--- a/sumpy/p2e.py
+++ b/sumpy/p2e.py
@@ -77,19 +77,20 @@ class P2E(object):
                 complex_dtype=np.complex128  # FIXME
                 )
 
+        from sumpy.tools import gather_source_arguments
         arguments = (
                 [
-                    lp.GlobalArg("sources", None, shape=(self.dim, "nsources")),
-                    lp.GlobalArg("strengths", None, shape="nsources"),
+                    lp.GlobalArg("sources", None, shape=(self.dim, "nsrc")),
+                    lp.GlobalArg("strengths", None, shape="nsrc"),
                     lp.GlobalArg("box_source_starts,box_source_counts_nonchild",
                         None, shape=None),
                     lp.GlobalArg("centers", None, shape="dim, nboxes"),
                     lp.GlobalArg("expansions", None,
                         shape=("nboxes", len(coeff_names))),
                     lp.ValueArg("nboxes", np.int32),
-                    lp.ValueArg("nsources", np.int32),
+                    lp.ValueArg("nsrc", np.int32),
                     "..."
-                ] + self.expansion.get_args())
+                ] + gather_source_arguments([self.expansion]))
 
         loopy_knl = lp.make_kernel(self.device,
                 [
diff --git a/sumpy/p2p.py b/sumpy/p2p.py
index a581a45ee0559de7d87ca10522373371bc999d13..e7129684e3336fdbe268b4e52fd98acacf082d07 100644
--- a/sumpy/p2p.py
+++ b/sumpy/p2p.py
@@ -83,7 +83,7 @@ class P2P(KernelComputation):
                     * var("strength_%d" % self.strength_usage[i])[var("isrc")]
                 for i, name in enumerate(result_names)]
 
-        from sumpy.tools import gather_arguments
+        from sumpy.tools import gather_source_arguments
         arguments = (
                 [
                     lp.GlobalArg("src", None,
@@ -98,7 +98,7 @@ class P2P(KernelComputation):
                 ]+[
                     lp.GlobalArg("result_%d" % i, dtype, shape="ntgt", order="C")
                     for i, dtype in enumerate(self.value_dtypes)
-                ] + gather_arguments(self.kernels))
+                ] + gather_source_arguments(self.kernels))
 
         if self.exclude_self:
             from pymbolic.primitives import If, ComparisonOperator, Variable
diff --git a/sumpy/qbx.py b/sumpy/qbx.py
index 36a3cb02bc56c58dab943af6d6947820101dce3f..8ed6382711fda7e1db472c27daac86f247c61edd 100644
--- a/sumpy/qbx.py
+++ b/sumpy/qbx.py
@@ -118,7 +118,7 @@ class LayerPotentialBase(KernelComputation):
                 * self.get_strength_or_not(isrc_sym, i)
                 for i, name in enumerate(result_names)]
 
-        from sumpy.tools import gather_arguments
+        from sumpy.tools import gather_source_arguments
         arguments = (
                 [
                     lp.GlobalArg("src", None,
@@ -130,7 +130,7 @@ class LayerPotentialBase(KernelComputation):
                     lp.ValueArg("nsrc", None),
                     lp.ValueArg("ntgt", None),
                 ] + self.get_input_and_output_arguments()
-                + gather_arguments(self.kernels))
+                + gather_source_arguments(self.kernels))
 
         loopy_knl = lp.make_kernel(self.device,
                 "{[isrc,itgt,idim]: 0<=itgt<ntgt and 0<=isrc<nsrc "
diff --git a/sumpy/tools.py b/sumpy/tools.py
index 4dc1f80118a223a3e58f5e57fba137c0a3fcb7f6..446795091f68cc37c1ea058b93111a788ea736a4 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -94,6 +94,16 @@ def gather_arguments(kernel_likes):
     return sorted(result.itervalues(), key=lambda arg: arg.name)
 
 
+def gather_source_arguments(kernel_likes):
+    result = {}
+    for knl in kernel_likes:
+        for arg in knl.get_args() + knl.get_source_args():
+            result[arg.name] = arg
+            # FIXME: possibly check that arguments match before overwriting
+
+    return sorted(result.itervalues(), key=lambda arg: arg.name)
+
+
 # {{{  KernelComputation
 
 class KernelComputation(object):
diff --git a/test/test_kernels.py b/test/test_kernels.py
index 67bdd4526aa0f5594ca3c031749a47b038a5a740..00324bb019b9d6be5e8bb830e7d64f7d9cc8bfca 100644
--- a/test/test_kernels.py
+++ b/test/test_kernels.py
@@ -34,7 +34,8 @@ from pyopencl.tools import (  # noqa
 
 from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion
 from sumpy.expansion.local import VolumeTaylorLocalExpansion, H2DLocalExpansion
-from sumpy.kernel import LaplaceKernel, HelmholtzKernel, AxisTargetDerivative
+from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative,
+        DirectionalSourceDerivative)
 
 import logging
 logger = logging.getLogger(__name__)
@@ -90,7 +91,11 @@ def test_p2p(ctx_getter):
     (HelmholtzKernel(2), VolumeTaylorLocalExpansion),
     (HelmholtzKernel(2), H2DLocalExpansion),
     ])
-def test_p2e2p(ctx_getter, knl, expn_class, order):
+@pytest.mark.parametrize("with_source_derivative", [
+    False,
+    True
+    ])
+def test_p2e2p(ctx_getter, knl, expn_class, order, with_source_derivative):
     #logging.basicConfig(level=logging.INFO)
 
     ctx = ctx_getter()
@@ -105,6 +110,9 @@ def test_p2e2p(ctx_getter, knl, expn_class, order):
     if isinstance(knl, HelmholtzKernel):
         extra_kwargs["k"] = 0.05
 
+    if with_source_derivative:
+        knl = DirectionalSourceDerivative(knl, "dir_vec")
+
     out_kernels = [
             knl,
             AxisTargetDerivative(0, knl),
@@ -136,6 +144,12 @@ def test_p2e2p(ctx_getter, knl, expn_class, order):
     box_source_starts = np.array([0], dtype=np.int32)
     box_source_counts_nonchild = np.array([nsources], dtype=np.int32)
 
+    extra_source_kwargs = extra_kwargs.copy()
+    if with_source_derivative:
+        alpha = np.linspace(0, 2*np.pi, nsources, np.float64)
+        dir_vec = np.vstack([np.cos(alpha), np.sin(alpha)])
+        extra_source_kwargs["dir_vec"] = dir_vec
+
     for h in h_values:
         from sumpy.visualization import FieldPlotter
         if issubclass(expn_class, LocalExpansionBase):
@@ -161,7 +175,7 @@ def test_p2e2p(ctx_getter, knl, expn_class, order):
                 sources=sources,
                 strengths=strengths,
                 #flags="print_hl_cl",
-                out_host=True, **extra_kwargs)
+                out_host=True, **extra_source_kwargs)
 
         # }}}
 
@@ -191,7 +205,7 @@ def test_p2e2p(ctx_getter, knl, expn_class, order):
                 queue,
                 targets, sources, (strengths,),
                 out_host=True,
-                **extra_kwargs)
+                **extra_source_kwargs)
 
         err_pot = la.norm((pot - pot_direct)/res**2)
         err_grad_x = la.norm((grad_x - grad_x_direct)/res**2)
@@ -246,6 +260,10 @@ def test_p2e2p(ctx_getter, knl, expn_class, order):
             slack += 1
             grad_slack += 1
 
+    if isinstance(knl, DirectionalSourceDerivative):
+        slack += 1
+        grad_slack += 2
+
     assert eoc_rec_pot.order_estimate() > tgt_order - slack
     assert eoc_rec_grad_x.order_estimate() > tgt_order_grad - grad_slack