Tibeuleu/FOC_Reduction
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

make use of existing scipy.special.gammaincc function for emission_center script

Browse Source
This commit is contained in:
Thibault Barnouin
2025-02-13 17:09:14 +01:00
parent 56f963ac7d
commit 450c9f3e59
2 changed files with 15 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
package/lib/utils.py
View File

@@ -88,8 +88,8 @@ def CenterConf(mask, PA, sPA):
2D matrix of same shape as input containing the confidence on the polarization 2D matrix of same shape as input containing the confidence on the polarization
computation between 0 and 1 for 2 parameters of interest (Q and U Stokes fluxes). computation between 0 and 1 for 2 parameters of interest (Q and U Stokes fluxes).
""" """
chi2 = np.full(PA.shape, np.nan) chi2 = np.full(PA.shape, np.nan, dtype=np.float64)
conf = np.full(PA.shape, -1.0) conf = np.full(PA.shape, -1.0, dtype=np.float64)
yy, xx = np.indices(PA.shape) yy, xx = np.indices(PA.shape)
def ideal(c): def ideal(c):
@@ -105,9 +105,9 @@ def CenterConf(mask, PA, sPA):
chi2[y, x] = chisq((x, y)) chi2[y, x] = chisq((x, y))
from scipy.optimize import minimize from scipy.optimize import minimize
from scipy.special import gammainc from scipy.special import gammaincc
conf[np.isfinite(PA)] = 1.0 - gammainc(0.5, 0.5 * chi2[np.isfinite(PA)]) conf[np.isfinite(PA)] = gammaincc(0.5, 0.5 * chi2[np.isfinite(PA)])
c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1] c0 = np.unravel_index(np.argmax(conf), conf.shape)[::-1]
result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])]) result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])])
if result.success: if result.success:

17
package/src/emission_center.py
View File

@@ -6,7 +6,7 @@ from sys import path as syspath
syspath.append(str(Path(__file__).parent.parent)) syspath.append(str(Path(__file__).parent.parent))
def main(infile, target=None, output_dir=None): def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None):
from os.path import join as pathjoin from os.path import join as pathjoin
import numpy as np import numpy as np
@@ -34,19 +34,22 @@ def main(infile, target=None, output_dir=None):
Stokes["POL_DEG_DEBIASED"].data[Stokes["POL_DEG_ERR"].data > 0.0] / Stokes["POL_DEG_ERR"].data[Stokes["POL_DEG_ERR"].data > 0.0] Stokes["POL_DEG_DEBIASED"].data[Stokes["POL_DEG_ERR"].data > 0.0] / Stokes["POL_DEG_ERR"].data[Stokes["POL_DEG_ERR"].data > 0.0]
) )
Stokescentconf, Stokescenter = CenterConf(Stokesconf > 0.99, Stokes["POL_ANG"].data, Stokes["POL_ANG_ERR"].data) if P_cut < 1.0:
Stokescentconf, Stokescenter = CenterConf(Stokesconf > P_cut, Stokes["POL_ANG"].data, Stokes["POL_ANG_ERR"].data)
else:
Stokescentconf, Stokescenter = CenterConf(Stokessnr > P_cut, Stokes["POL_ANG"].data, Stokes["POL_ANG_ERR"].data)
Stokespos = WCS(Stokes[0].header).pixel_to_world(*Stokescenter) Stokespos = WCS(Stokes[0].header).pixel_to_world(*Stokescenter)
if target is None: if target is None:
target = Stokes[0].header["TARGNAME"] target = Stokes[0].header["TARGNAME"]
fig = figure(figsize=(8, 8.5), layout="constrained") fig = figure(figsize=(8, 8.5), layout="constrained")
fig, ax = polarization_map(Stokes, P_cut=0.99, step_vec=1, scale_vec=3, display="pf", fig=fig, width=0.33, linewidth=0.5) fig, ax = polarization_map(Stokes, P_cut=P_cut, step_vec=1, scale_vec=3, display=display, fig=fig, width=0.33, linewidth=0.5)
ax.plot(*Stokescenter, marker="+", color="k", lw=3) ax.plot(*Stokescenter, marker="+", color="k", lw=3)
ax.plot(*Stokescenter, marker="+", color="gray", lw=1.5, label="Best confidence for center: {0}".format(Stokespos.to_string("hmsdms"))) ax.plot(*Stokescenter, marker="+", color="red", lw=1.5, label="Best confidence for center: {0}".format(Stokespos.to_string("hmsdms")))
ax.contour(Stokescentconf, [0.01], colors="k", linewidths=3) ax.contour(Stokescentconf, [0.01], colors="k", linewidths=3)
confcentcont = ax.contour(Stokescentconf, [0.01], colors="gray") confcentcont = ax.contour(Stokescentconf, [0.01], colors="red")
# confcont = ax.contour(Stokesconf, [0.9905], colors="r") # confcont = ax.contour(Stokesconf, [0.9905], colors="r")
# snr3cont = ax.contour(Stokessnr, [3.0], colors="b", linestyles="dashed") # snr3cont = ax.contour(Stokessnr, [3.0], colors="b", linestyles="dashed")
# snr4cont = ax.contour(Stokessnr, [4.0], colors="b") # snr4cont = ax.contour(Stokessnr, [4.0], colors="b")
@@ -75,7 +78,9 @@ if __name__ == "__main__":
parser = argparse.ArgumentParser(description="Look for the center of emission for a given reduced observation") parser = argparse.ArgumentParser(description="Look for the center of emission for a given reduced observation")
parser.add_argument("-t", "--target", metavar="targetname", required=False, help="the name of the target", type=str, default=None) parser.add_argument("-t", "--target", metavar="targetname", required=False, help="the name of the target", type=str, default=None)
parser.add_argument("-f", "--file", metavar="path", required=False, help="The full or relative path to the data product", type=str, default=None) parser.add_argument("-f", "--file", metavar="path", required=False, help="The full or relative path to the data product", type=str, default=None)
parser.add_argument("-c", "--pcut", metavar="pcut", required=False, help="The polarization cut for the data mask", type=float, default=0.99)
parser.add_argument("-d", "--display", metavar="display", required=False, help="The map on which to display info", type=str, default="pf")
parser.add_argument("-o", "--output_dir", metavar="directory_path", required=False, help="output directory path for the plots", type=str, default="./data") parser.add_argument("-o", "--output_dir", metavar="directory_path", required=False, help="output directory path for the plots", type=str, default="./data")
args = parser.parse_args() args = parser.parse_args()
exitcode = main(infile=args.file, target=args.target, output_dir=args.output_dir) exitcode = main(infile=args.file, P_cut=args.pcut, target=args.target, display=args.display, output_dir=args.output_dir)
print("Written to: ", exitcode) print("Written to: ", exitcode)