diff --git a/src/yapsut/periodic_stats.py b/src/yapsut/periodic_stats.py
index 9697b4b3906e58458c6eb252f701b2da411ccb63..6f2972524a0cadf2741668cbe6ca0ee6f19fec4e 100644
--- a/src/yapsut/periodic_stats.py
+++ b/src/yapsut/periodic_stats.py
@@ -3,27 +3,38 @@
 import numpy as np
 
 def periodic_stats(alpha,deg=False) :
-    """ statistics of a periodic alpha (radiants)
+   """ statistics of a periodic alpha (radiants)
     
-    If alpha is a vector of angles (radiants) find the average number and the std deviation avoiding the problem of the switch between -180,180 or 360 and 0
+   If alpha is a vector of angles (radiants) find the average number and the std deviation avoiding the problem of the switch between -180,180 or 360 and 0
     
-    example: alpha=np.array([-180,180]) 
+   example: alpha=np.array([-180,180]) 
     
-    """
-    u=np.cos(alpha)
-    v=np.sin(alpha)
-    #
-    rr=np.array([np.nanmean(u),np.nanmean(v)])
-    rr=rr/rr.dot(rr)**0.5
-    #
-    up=rr[0]*u+rr[1]*v ;
-    vp=-rr[1]*u+rr[0]*v
-    #
-    ap=np.rad2deg(np.arctan2(vp,up).std())
-    #
-    return np.arctan2(rr[1],rr[0]), np.nanvar(np.arctan2(vp,up))**0.5
+   """
+   u=np.cos(alpha)
+   v=np.sin(alpha)
+   #
+   rr=np.array([np.nanmean(u),np.nanmean(v)])
+   rr=rr/rr.dot(rr)**0.5
+   #
+   up=rr[0]*u+rr[1]*v ;
+   vp=-rr[1]*u+rr[0]*v
+   #
+   ap=np.rad2deg(np.arctan2(vp,up).std())
+   #
+   return np.arctan2(rr[1],rr[0]), np.nanvar(np.arctan2(vp,up))**0.5
 
 if __name__=='__main__' :
-   alpha_deg=np.arange(350,371)
+   print("Test1")
+   alpha=np.arange(350,371)
+   print(alpha)
    print(np.deg2rad(alpha-360).mean(),np.deg2rad(alpha-360).std())
-   print(periodic_stats(np.deg2rad(alpha)))    
+   print(periodic_stats(np.deg2rad(alpha)))
+   print()
+
+   print("Test2")
+   alpha=np.array([-170,-171,-172,-173,-174,-175,-176,-177,-178,-179,-180,180,179,178,177,176,175,174,173,172,171,170])
+   print(alpha)
+   print(np.mod(alpha,360))
+   print(np.deg2rad(np.mod(alpha,360)).mean(),np.deg2rad(np.mod(alpha,360)).std())
+   print(periodic_stats(np.deg2rad(alpha)))
+   print()