--- a/pymedcouple
+++ b/pymedcouple
@@ -112,17 +112,11 @@
     if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)):
         return 1.0
 
-    if debug:
-        print "Zmed = ", Zmed
-        
     # Centre Z wrt median, so that median(Z) = 0.
     Z = [z - Zmed for z in Z]
 
     # Scale inside [-0.5, 0.5], for greater numerical stability.
     Zden = 2*max(Z[0], -Z[-1])
-    if debug:
-        print "Zden =", Zden
-    
     Z = [z/Zden for z in Z]
     Zmed /= Zden
     
@@ -132,19 +126,9 @@
     Zplus   = [z for z in Z if z >= -Zeps]
     Zminus  = [z for z in Z if Zeps >= z]
 
-    if debug:
-        for zp in Zplus:
-            print "Zplus : %g" % zp
-        for zm in Zminus:
-            print "Zminus: %g" % zm
-
     n_plus = len(Zplus)
     n_minus = len(Zminus)
 
-    if debug:
-        print "n_plus =", n_plus
-        print "n_minus = ", n_minus
-
     def h_kern(i, j, debug=False):
         """Kernel function h for the medcouple, closing over the values of
         Zplus and Zminus just defined above.
@@ -161,9 +145,6 @@
         else:
             h = (a+b)/(a-b)
 
-        if False:
-            print "a = {a}, b = {b}, h({i},{j}) = {h}".format(**locals())
-
         return h
 
     # Init left and right borders
@@ -174,15 +155,9 @@
     Rtot = n_minus*n_plus
     mid_idx = (Rtot-1)//2
 
-    if debug:
-        print "mid_idx = ", mid_idx
-
     # kth pair algorithm (Johnson & Mizoguchi)
     while Rtot - Ltot > n_plus:
-        if debug:
-            print "L = ", L
-            print "R = ", R
-        
+       
         # First, compute the median inside the given bounds
         # (Be stingy, reuse same generator)
         [I1, I2] = tee(i for i in xrange(0, n_plus) if L[i] <= R[i])
@@ -190,8 +165,6 @@
         A = [h_kern(i, (L[i] + R[i])//2) for i in I1]
         W = [R[i] - L[i] + 1 for i in I2]
         Am = wmedian(A,W)
-        if debug:
-            print "Am = ", Am
             
         Am_eps = eps1*(eps1 + abs(Am))
 
@@ -218,10 +191,6 @@
         # the whole matrix may be.
         sumP = sum(P) + len(P)
         sumQ = sum(Q)
-        if debug:
-            print "P: ", sumP, P
-            print "Q: ", sumQ, Q
-            print
         
         if mid_idx <= sumP - 1:
             R = P
@@ -236,9 +205,6 @@
     # Didn't find the median, but now we have a very small search
     # space to find it in, just between the left and right boundaries.
     # This space is of size Rtot - Ltot which is <= n_plus
-    if debug:
-        print "L =", L
-        print "R =", R
     A = []
     for (i, (l, r)) in enumerate(izip(L, R)):
         for j in xrange(l, r + 1):
@@ -246,8 +212,6 @@
     A.sort()
     A.reverse()
 
-    if debug:
-        print "len(A) = ", len(A)
     Am =  A[mid_idx - Ltot]
     return Am
     
@@ -260,11 +224,13 @@
     else:
         return 0
 
+
 def main():
     import sys
     fname = sys.argv[1]
     with open(fname) as f:
         data = [float(x) for x in f.readlines() if x.strip() != ""]
+
     print "%.16g" % medcouple_1d(data)
 
 if __name__ == "__main__":