From bf5373d4e03c29b4122b427256ce49da0927dc5e Mon Sep 17 00:00:00 2001
From: Thibault Barnouin <thibault.barnouin@astro.unistra.fr>
Date: Fri, 6 Sep 2024 17:28:01 +0200
Subject: [PATCH] add fonction to compute confidence map of center of emission
 location

---
 package/lib/utils.py | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)

diff --git a/package/lib/utils.py b/package/lib/utils.py
index 858b8aa..204352a 100755
--- a/package/lib/utils.py
+++ b/package/lib/utils.py
@@ -39,6 +39,23 @@ def PCconf(QN, UN, QN_ERR, UN_ERR):
     conf[mask] = 1.0 - np.exp(-0.5 * chi2[mask])
     return conf
 
+def Centerconf(mask, PA, sPA):
+    """
+    Compute the confidence map for the position of the center of emission.
+    """
+    chi2 = np.full(PA.shape, np.nan)
+    conf = np.full(PA.shape, -1.0)
+    yy, xx = np.indices(PA.shape)
+    def ideal(c):
+        itheta = np.degrees(np.arctan((yy+0.5-c[1])/(xx+0.5-c[0])))
+        itheta[np.isnan(itheta)] = PA[np.isnan(itheta)]
+        return princ_angle(itheta)
+    def chisq(c):
+        return np.sum((princ_angle(PA[mask])-ideal((x,y))[mask])**2/sPA[mask]**2)/np.sum(mask)
+    for x,y in zip(xx[np.isfinite(PA)],yy[np.isfinite(PA)]):
+        chi2[y,x] = chisq((x,y))
+    conf[mask] = 1.0 - np.exp(-0.5*chi2[mask])
+    return conf
 
 def sci_not(v, err, rnd=1, out=str):
     """