diff --git a/setup.py b/setup.py
index c793e73cb43393db60b5951e870b04df71b6c2fe..58cc30f016edd4430860f74f38395b17209de635 100644
--- a/setup.py
+++ b/setup.py
@@ -37,7 +37,7 @@ setup(name="sumpy",
       author="Rio Yokota, Andreas Kloeckner, Matthew Knepley",
       author_email="yokota@bu.edu",
       license = "MIT",
-      packages=["sumpy"],
+      packages=["sumpy", "sumpy.symbolic"],
 
       # 2to3 invocation
       cmdclass={'build_py': build_py})
diff --git a/sumpy/symbolic/__init__.py b/sumpy/symbolic/__init__.py
index 9c184b1770ea4a6fe2b6855abb6b44144c6f21c9..58d4f8d532eb8f80b1020bd8f2143cd41e9ca3d9 100644
--- a/sumpy/symbolic/__init__.py
+++ b/sumpy/symbolic/__init__.py
@@ -61,7 +61,7 @@ class IdentityMapper(SympyMapper):
 
     map_Symbol = map_Integer
     map_Real = map_Integer
-
+    map_ImaginaryUnit = map_Integer
 
 
 
@@ -129,6 +129,9 @@ class CSEMapper(IdentityMapper):
         else:
             return IdentityMapper.map_Rational(self, expr)
 
+    def map_Subs(self, expr):
+        return expr
+
 
 
 
@@ -154,7 +157,7 @@ def eliminate_common_subexpressions(exprs, sym_gen):
             else:
                 iterable = node
 
-            for arg in node.args:
+            for arg in iterable:
                 gather_use_counts(arg)
 
     for expr in exprs:
@@ -204,20 +207,58 @@ def diff_multi_index(expr, multi_index, var_name):
 
 
 
-def make_coulomb_kernel_in(var_name, dimensions):
-    from sumpy.symbolic import make_sym_vector
-    dist = make_sym_vector(var_name, dimensions)
+def make_coulomb_kernel(dist_vec):
+    dimensions = len(dist_vec)
+    r = sp.sqrt((dist_vec.T*dist_vec)[0,0])
 
     if dimensions == 2:
-        return sp.log(sp.sqrt((dist.T*dist)[0,0]))
+        return sp.log(r)
     elif dimensions == 3:
-        return 1/sp.sqrt((dist.T*dist)[0,0])
+        return 1/r
     else:
         raise RuntimeError("unsupported dimensionality")
 
 
 
 
+
+def make_helmholtz_kernel(dist_vec):
+    dimensions = len(dist_vec)
+    r = sp.sqrt((dist_vec.T*dist_vec)[0,0])
+
+    i = sp.sqrt(-1)
+    k = sp.Symbol("k")
+
+    if dimensions == 2:
+        return i/4 * sp.Function("H1_0")(k*r)
+    elif dimensions == 3:
+        return sp.exp(i*k*r)/r
+    else:
+        raise RuntimeError("unsupported dimensionality")
+
+
+
+
+def make_coulomb_kernel_in(var_name, dimensions):
+    from warnings import warn
+    warn("make_coulomb_kernel_in is deprecated", DeprecationWarning, stacklevel=2)
+
+    from sumpy.symbolic import make_sym_vector
+    return make_coulomb_kernel(make_sym_vector(var_name, dimensions))
+
+
+
+
+def make_helmholtz_kernel_in(var_name, dimensions):
+    from warnings import warn
+    warn("make_helmholtz_kernel_in is deprecated", DeprecationWarning, stacklevel=2)
+
+    from sumpy.symbolic import make_sym_vector
+    return make_helmholtz_kernel(make_sym_vector(var_name, dimensions))
+
+
+
+
 def make_coulomb_kernel_ts(dimensions):
     old = make_sym_vector("d", dimensions)
     new = (make_sym_vector("t", dimensions)
diff --git a/sumpy/symbolic/codegen.py b/sumpy/symbolic/codegen.py
index a30523ca2e19ec2c044de07391ce59ef6aea50fb..8b1523f3ec65e8239292d38dfa9b2127200c3e8b 100644
--- a/sumpy/symbolic/codegen.py
+++ b/sumpy/symbolic/codegen.py
@@ -88,6 +88,10 @@ class PowRewriter(IdentityMapper):
 
         return IdentityMapper.map_Pow(self, expr)
 
+    def map_Subs(self, expr):
+        return expr
+
+
 
 
 
diff --git a/sumpy/tools.py b/sumpy/tools.py
index b9d6c8baa304a550bccb006529b24c4b6e16f910..f9f102c2766901d5f16f3cb90d0fceb78548bd47 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -21,6 +21,12 @@ def mi_power(vector, mi):
         result *= vec_i**mi_i
     return result
 
+def mi_derivative(expr, vector, mi):
+    for mi_i, vec_i in zip(mi, vector):
+        expr = expr.diff(vec_i, mi_i)
+    return expr
+
+
 # }}}
 
 
diff --git a/test/test_kernels.py b/test/test_kernels.py
index eb466aea1214532fe273ac03ea4a37218e9dcc6c..cb5d585efefd6405d3d62fcb2aab9dad4c2bda89 100644
--- a/test/test_kernels.py
+++ b/test/test_kernels.py
@@ -156,6 +156,8 @@ def test_p2m2p(ctx_getter):
 
     assert la.norm((potential_dev-potential_dev_direct).get())/res**2 < 1e-3
 
+    # }}}
+
 
 
 
@@ -263,6 +265,8 @@ def test_p2m2m2p(ctx_getter):
 
     assert la.norm((potential_dev-potential_dev_direct).get())/res**2 < 1e-3
 
+    # }}}
+