From 703ad44ebd98bc2d6ece5af52905a90f6014d67f Mon Sep 17 00:00:00 2001
From: Michele Maris <michele.maris@inaf.it>
Date: Thu, 13 Jul 2023 19:18:46 +0200
Subject: [PATCH] u

---
 src/yapsut/__init__.py       |  1 +
 src/yapsut/periodic_stats.py | 29 +++++++++++++++++++++++++++++
 2 files changed, 30 insertions(+)
 create mode 100644 src/yapsut/periodic_stats.py

diff --git a/src/yapsut/__init__.py b/src/yapsut/__init__.py
index 2728e81..8c11f93 100644
--- a/src/yapsut/__init__.py
+++ b/src/yapsut/__init__.py
@@ -13,6 +13,7 @@ from .DummyCLI import DummyCLI
 from .nearly_even_split import nearly_even_split
 from .geometry import ElipsoidTriangulation 
 from .template_subs import discoverKeys, templateFiller
+from .periodic_stats import periodic_stats
 
 #from import .grep
 #import .intervalls
diff --git a/src/yapsut/periodic_stats.py b/src/yapsut/periodic_stats.py
new file mode 100644
index 0000000..9697b4b
--- /dev/null
+++ b/src/yapsut/periodic_stats.py
@@ -0,0 +1,29 @@
+""" periodic stats is a collection of tools for statistics over periodic quantities """
+
+import numpy as np
+
+def periodic_stats(alpha,deg=False) :
+    """ 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
+    
+    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
+
+if __name__=='__main__' :
+   alpha_deg=np.arange(350,371)
+   print(np.deg2rad(alpha-360).mean(),np.deg2rad(alpha-360).std())
+   print(periodic_stats(np.deg2rad(alpha)))    
-- 
GitLab