diff --git a/snippets/circstats.py b/snippets/circstats.py
index 3c6deb1..f2be520 100644
--- a/snippets/circstats.py
+++ b/snippets/circstats.py
@@ -54,7 +54,7 @@
 ]
 
 import numpy as np
-from scipy.stats import chi2
+from scipy.stats import chi2, norm
 from scipy.special import iv
 
 
@@ -742,14 +742,20 @@ def _corrcc(alpha1, alpha2, nanrobust, axis=None):
     else:
         sumfunc = lambda x: np.sum(x, axis=axis)
 
-    num = 4 * (sumfunc(c1*c2) * sumfunc(s1*s2) -
-               sumfunc(c1*s2) * sumfunc(s1*c2))
-    den = np.sqrt((n**2 - sumfunc(c1_2)**2 - sumfunc(s1_2)**2) *
-                  (n**2 - sumfunc(c2_2)**2 - sumfunc(s2_2)**2))
+    num = 4 * (sumfunc(c1*c2) * sumfunc(s1*s2) - sumfunc(c1*s2) * sumfunc(s1*c2))
+    
+    den = np.sqrt((n**2 - sumfunc(c1_2)**2 - sumfunc(s1_2)**2) * (n**2 - sumfunc(c2_2)**2 - sumfunc(s2_2)**2))
 
     rho = num / den
 
-    return rho
+    l20 = np.mean(np.power(np.sin(alpha1 - mean(alpha1)),2))
+    l02 = np.mean(np.power(np.sin(alpha2 - mean(alpha2)),2))
+    l22 = np.mean(np.power(np.sin(alpha1 - mean(alpha1)),2) * np.power(np.sin(alpha2 - mean(alpha2)), 2));
+
+    ts = np.sqrt((n * l20 * l02)/l22) * rho;
+    pval = 2.0 * (1.0 - norm().cdf(np.abs(ts)));
+    
+    return (rho, pval)
 
 
 def corrcc(alpha1, alpha2, axis=None, nancheck=True):