diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py
index ca349b59..98d88107 100644
--- a/GPy/kern/constructors.py
+++ b/GPy/kern/constructors.py
@@ -292,7 +292,8 @@ except ImportError:
 if sympy_available:
     from parts.sympykern import spkern
     from sympy.parsing.sympy_parser import parse_expr
-
+    from GPy.util.symbolic import sinc
+    
     def rbf_sympy(input_dim, ARD=False, variance=1., lengthscale=1.):
         """
         Radial Basis Function covariance.
@@ -312,6 +313,27 @@ if sympy_available:
             f =  variance*sp.exp(-dist/(2*lengthscale**2))
         return kern(input_dim, [spkern(input_dim, f, name='rbf_sympy')])
 
+    def sinc(input_dim, ARD=False, variance=1., lengthscale=1.):
+        """
+        TODO: Not clear why this isn't working, suggests argument of sinc is not a number.
+        sinc covariance funciton
+        """
+        X = [sp.var('x%i' % i) for i in range(input_dim)]
+        Z = [sp.var('z%i' % i) for i in range(input_dim)]
+        variance = sp.var('variance',positive=True)
+        if ARD:
+            lengthscales = [sp.var('lengthscale_%i' % i, positive=True) for i in range(input_dim)]
+            dist_string = ' + '.join(['(x%i-z%i)**2/lengthscale_%i**2' % (i, i, i) for i in range(input_dim)])
+            dist = parse_expr(dist_string)
+            f =  variance*sinc(sp.pi*sp.sqrt(dist))
+        else:
+            lengthscale = sp.var('lengthscale',positive=True)
+            dist_string = ' + '.join(['(x%i-z%i)**2' % (i, i) for i in range(input_dim)])
+            dist = parse_expr(dist_string)
+            f =  variance*sinc(sp.pi*sp.sqrt(dist)/lengthscale)
+            
+        return kern(input_dim, [spkern(input_dim, f, name='sinc')])
+
     def sympykern(input_dim, k,name=None):
         """
         A base kernel object, where all the hard work in done by sympy.
diff --git a/GPy/kern/parts/sympy_helpers.cpp b/GPy/kern/parts/sympy_helpers.cpp
index ec8ac87b..76dba4eb 100644
--- a/GPy/kern/parts/sympy_helpers.cpp
+++ b/GPy/kern/parts/sympy_helpers.cpp
@@ -9,3 +9,17 @@ double DiracDelta(double x){
 double DiracDelta(double x,int foo){
     return 0.0;
 };
+
+double sinc(double x){
+  if (x==0)
+    return 1.0;
+  else 
+    return sin(x)/x;
+}
+
+double sinc_grad(double x){
+  if (x==0)
+    return 0.0;
+  else 
+    return (x*cos(x) - sin(x))/(x*x);
+}
diff --git a/GPy/kern/parts/sympy_helpers.h b/GPy/kern/parts/sympy_helpers.h
index 29244eca..d5b495ca 100644
--- a/GPy/kern/parts/sympy_helpers.h
+++ b/GPy/kern/parts/sympy_helpers.h
@@ -1,3 +1,6 @@
 #include <math.h>
 double DiracDelta(double x);
 double DiracDelta(double x, int foo);
+
+double sinc(double x);
+double sinc_grad(double x);
diff --git a/GPy/util/symbolic.py b/GPy/util/symbolic.py
new file mode 100644
index 00000000..f4f5fda0
--- /dev/null
+++ b/GPy/util/symbolic.py
@@ -0,0 +1,32 @@
+from sympy import Function, S, oo, I, cos, sin
+
+
+class sinc_grad(Function):
+    nargs = 1
+    
+    def fdiff(self, argindex=1):
+        return ((2-x*x)*sin(self.args[0]) - 2*x*cos(x))/(x*x*x)
+    
+    @classmethod
+    def eval(cls, x):
+        if x is S.Zero:
+            return S.Zero
+        else:
+            return (x*cos(x) - sin(x))/(x*x)
+            
+class sinc(Function):
+    
+    nargs = 1
+    
+    def fdiff(self, argindex=1):
+        return sinc_grad(self.args[0])
+    
+    @classmethod
+    def eval(cls, x):
+        if x is S.Zero:
+            return S.One
+        else:
+            return sin(x)/x
+    
+    def _eval_is_real(self):
+        return self.args[0].is_real