diff --git a/package/lib/utils.py b/package/lib/utils.py
index 204352a..e28bbff 100755
--- a/package/lib/utils.py
+++ b/package/lib/utils.py
@@ -39,23 +39,38 @@ 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):
+
+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)]
+        itheta = np.full(PA.shape, np.nan)
+        itheta[(xx + 0.5) != c[0]] = np.degrees(np.arctan((yy[(xx + 0.5) != c[0]] + 0.5 - c[1]) / (xx[(xx + 0.5) != c[0]] + 0.5 - c[0])))
+        itheta[(xx + 0.5) == c[0]] = PA[(xx + 0.5) == c[0]]
         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
+        return np.sum((princ_angle(PA[mask]) - ideal((c[0], c[1]))[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))
+
+    from scipy.optimize import minimize
+    from scipy.special import gammainc
+
+    conf[np.isfinite(PA)] = 1.0 - gammainc(0.5, 0.5 * chi2[np.isfinite(PA)])
+    result = minimize(chisq, np.array(PA.shape) / 2.0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])])
+    if result.success:
+        print("Center of emission found")
+    else:
+        print("Center of emission not found")
+    return conf, result.x
+
 
 def sci_not(v, err, rnd=1, out=str):
     """
diff --git a/package/test_center.py b/package/test_center.py
new file mode 100644
index 0000000..54c21d4
--- /dev/null
+++ b/package/test_center.py
@@ -0,0 +1,99 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from astropy.io.fits import open as fits_open
+from astropy.wcs import WCS
+from lib.utils import CenterConf, PCconf
+from matplotlib.colors import LogNorm
+from matplotlib.patches import Rectangle
+
+levelssnr = np.array([3.0, 4.0])
+levelsconf = np.array([0.99])
+
+NGC1068 = fits_open("./data/NGC1068/5144/NGC1068_FOC_b0.05arcsec_c0.07arcsec.fits")
+NGC1068conf = PCconf(
+    NGC1068["Q_STOKES"].data / NGC1068["I_STOKES"].data,
+    NGC1068["U_STOKES"].data / NGC1068["I_STOKES"].data,
+    np.sqrt(NGC1068["IQU_COV_MATRIX"].data[1, 1]) / NGC1068["I_STOKES"].data,
+    np.sqrt(NGC1068["IQU_COV_MATRIX"].data[2, 2]) / NGC1068["I_STOKES"].data,
+)
+NGC1068mask = NGC1068["DATA_MASK"].data.astype(bool)
+NGC1068snr = np.full(NGC1068mask.shape, np.nan)
+NGC1068snr[NGC1068["POL_DEG_ERR"].data > 0.0] = (
+    NGC1068["POL_DEG_DEBIASED"].data[NGC1068["POL_DEG_ERR"].data > 0.0] / NGC1068["POL_DEG_ERR"].data[NGC1068["POL_DEG_ERR"].data > 0.0]
+)
+
+NGC1068centconf, NGC1068center = CenterConf(NGC1068conf > 0.99, NGC1068["POL_ANG"].data, NGC1068["POL_ANG_ERR"].data)
+
+figngc, axngc = plt.subplots(1, 2, layout="tight", figsize=(18,9), subplot_kw=dict(projection=WCS(NGC1068[0].header)))
+
+axngc[0].set(xlabel="RA", ylabel="DEC", title="NGC1069 intensity map with SNR and confidence contours")
+imngc = axngc[0].imshow(NGC1068["I_STOKES"].data * NGC1068["I_STOKES"].header["PHOTFLAM"], norm=LogNorm(), cmap="inferno")
+ngcsnrcont = axngc[0].contour(NGC1068snr, levelssnr, colors="b")
+ngcconfcont = axngc[0].contour(NGC1068conf, levelsconf, colors="r")
+ngcconfcenter = axngc[0].plot(*np.unravel_index(np.argmax(NGC1068centconf), NGC1068centconf.shape)[::-1], "k+", label="Best confidence for center")
+ngcconfcentcont = axngc[0].contour(NGC1068centconf, 1.-levelsconf, colors="k")
+handles, labels = axngc[0].get_legend_handles_labels()
+labels.append("SNR contours")
+handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=ngcsnrcont.collections[0].get_edgecolor()[0]))
+labels.append("CONF99 contour")
+handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=ngcconfcont.collections[0].get_edgecolor()[0]))
+labels.append("Center CONF99 contour")
+handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=ngcconfcentcont.collections[0].get_edgecolor()[0]))
+axngc[0].legend(handles=handles, labels=labels)
+
+axngc[1].set(xlabel="RA", ylabel="DEC", title="Location of the nucleus confidence map")
+ngccent = axngc[1].imshow(NGC1068centconf, vmin=0.0, cmap="inferno")
+ngccentcont = axngc[1].contour(NGC1068centconf, 1.-levelsconf, colors="grey")
+ngccentcenter = axngc[1].plot(*np.unravel_index(np.argmax(NGC1068centconf), NGC1068centconf.shape)[::-1], "k+", label="Best confidence for center")
+handles, labels = axngc[1].get_legend_handles_labels()
+labels.append("CONF99 contour")
+handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=ngccentcont.collections[0].get_edgecolor()[0]))
+axngc[1].legend(handles=handles, labels=labels)
+
+figngc.savefig("NGC1068_center.pdf",dpi=150,facecolor="None")
+
+###################################################################################################
+
+MRK463E = fits_open("./data/MRK463E/5960/MRK463E_FOC_b0.05arcsec_c0.07arcsec.fits")
+MRK463Econf = PCconf(
+    MRK463E["Q_STOKES"].data / MRK463E["I_STOKES"].data,
+    MRK463E["U_STOKES"].data / MRK463E["I_STOKES"].data,
+    np.sqrt(MRK463E["IQU_COV_MATRIX"].data[1, 1]) / MRK463E["I_STOKES"].data,
+    np.sqrt(MRK463E["IQU_COV_MATRIX"].data[2, 2]) / MRK463E["I_STOKES"].data,
+)
+MRK463Emask = MRK463E["DATA_MASK"].data.astype(bool)
+MRK463Esnr = np.full(MRK463Emask.shape, np.nan)
+MRK463Esnr[MRK463E["POL_DEG_ERR"].data > 0.0] = (
+    MRK463E["POL_DEG_DEBIASED"].data[MRK463E["POL_DEG_ERR"].data > 0.0] / MRK463E["POL_DEG_ERR"].data[MRK463E["POL_DEG_ERR"].data > 0.0]
+)
+
+MRK463Ecentconf, MRK463Ecenter = CenterConf(MRK463Econf > 0.99, MRK463E["POL_ANG"].data, MRK463E["POL_ANG_ERR"].data)
+
+figmrk, axmrk = plt.subplots(1, 2, layout="tight", figsize=(18,9), subplot_kw=dict(projection=WCS(MRK463E[0].header)))
+
+axmrk[0].set(xlabel="RA", ylabel="DEC", title="NGC1069 intensity map with SNR and confidence contours")
+immrk = axmrk[0].imshow(MRK463E["I_STOKES"].data * MRK463E["I_STOKES"].header["PHOTFLAM"], norm=LogNorm(), cmap="inferno")
+mrksnrcont = axmrk[0].contour(MRK463Esnr, levelssnr, colors="b")
+mrkconfcont = axmrk[0].contour(MRK463Econf, levelsconf, colors="r")
+mrkconfcenter = axmrk[0].plot(*np.unravel_index(np.argmax(MRK463Ecentconf), MRK463Ecentconf.shape)[::-1], "k+", label="Best confidence for center")
+mrkconfcentcont = axmrk[0].contour(MRK463Ecentconf, 1.-levelsconf, colors="k")
+handles, labels = axmrk[1].get_legend_handles_labels()
+labels.append("SNR contours")
+handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=mrksnrcont.collections[0].get_edgecolor()[0]))
+labels.append("CONF99 contour")
+handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=mrkconfcont.collections[0].get_edgecolor()[0]))
+labels.append("Center CONF99 contour")
+handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=mrkconfcentcont.collections[0].get_edgecolor()[0]))
+axmrk[0].legend(handles=handles, labels=labels)
+
+axmrk[1].set(xlabel="RA", ylabel="DEC", title="Location of the nucleus confidence map")
+mrkcent = axmrk[1].imshow(MRK463Ecentconf, vmin=0.0, cmap="inferno")
+mrkcentcont = axmrk[1].contour(MRK463Ecentconf, 1.-levelsconf, colors="grey")
+mrkcentcenter = axmrk[1].plot(*np.unravel_index(np.argmax(MRK463Ecentconf), MRK463Ecentconf.shape)[::-1], "k+", label="Best confidence for center")
+handles, labels = axmrk[1].get_legend_handles_labels()
+labels.append("CONF99 contour")
+handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=mrkcentcont.collections[0].get_edgecolor()[0]))
+axmrk[1].legend(handles=handles, labels=labels)
+
+figmrk.savefig("MRK463E_center.pdf",dpi=150,facecolor="None")
+plt.show()