small improvments to ConfCenter et emission_center
This commit is contained in:
@@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA):
|
|||||||
|
|
||||||
conf[np.isfinite(PA)] = gammaincc(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])], method="trust-constr")
|
result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="TNC")
|
||||||
if result.success:
|
if result.success:
|
||||||
print("Center of emission found")
|
print("Center of emission found")
|
||||||
else:
|
else:
|
||||||
|
|||||||
@@ -23,7 +23,12 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None):
|
|||||||
stkI = Stokes["I_STOKES"].data
|
stkI = Stokes["I_STOKES"].data
|
||||||
QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan)
|
QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan)
|
||||||
for sflux, nflux in zip(
|
for sflux, nflux in zip(
|
||||||
[Stokes["Q_STOKES"].data, Stokes["U_STOKES"].data, np.sqrt(Stokes["IQU_COV_MATRIX"].data[1, 1]), np.sqrt(Stokes["IQU_COV_MATRIX"].data[2, 2])],
|
[
|
||||||
|
Stokes["Q_STOKES"].data,
|
||||||
|
Stokes["U_STOKES"].data,
|
||||||
|
np.sqrt(Stokes["IQU_COV_MATRIX"].data[1, 1]),
|
||||||
|
np.sqrt(Stokes["IQU_COV_MATRIX"].data[2, 2]),
|
||||||
|
],
|
||||||
[QN, UN, QN_ERR, UN_ERR],
|
[QN, UN, QN_ERR, UN_ERR],
|
||||||
):
|
):
|
||||||
nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0]
|
nflux[stkI > 0.0] = sflux[stkI > 0.0] / stkI[stkI > 0.0]
|
||||||
@@ -44,10 +49,25 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=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=P_cut, step_vec=1, scale_vec=5, display=display, fig=fig, width=0.33, linewidth=0.5)
|
fig, ax = polarization_map(
|
||||||
|
Stokes,
|
||||||
|
P_cut=P_cut,
|
||||||
|
step_vec=1,
|
||||||
|
scale_vec=5,
|
||||||
|
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="red", 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="red")
|
confcentcont = ax.contour(Stokescentconf, [0.01], colors="red")
|
||||||
# confcont = ax.contour(Stokesconf, [0.9905], colors="r")
|
# confcont = ax.contour(Stokesconf, [0.9905], colors="r")
|
||||||
@@ -62,11 +82,18 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None):
|
|||||||
# handles.append(Rectangle((0, 0), 1, 1, fill=False, ls="--", ec=snr3cont.get_edgecolor()[0]))
|
# handles.append(Rectangle((0, 0), 1, 1, fill=False, ls="--", ec=snr3cont.get_edgecolor()[0]))
|
||||||
# labels.append(r"$SNR_P \geq$ 4 contour")
|
# labels.append(r"$SNR_P \geq$ 4 contour")
|
||||||
# handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=snr4cont.get_edgecolor()[0]))
|
# handles.append(Rectangle((0, 0), 1, 1, fill=False, ec=snr4cont.get_edgecolor()[0]))
|
||||||
ax.legend(handles=handles, labels=labels, bbox_to_anchor=(-0.05, -0.02, 1.10, 0.01), loc="upper left", mode="expand", borderaxespad=0.0)
|
ax.legend(
|
||||||
|
handles=handles,
|
||||||
|
labels=labels,
|
||||||
|
bbox_to_anchor=(-0.05, -0.02, 1.10, 0.01),
|
||||||
|
loc="upper left",
|
||||||
|
mode="expand",
|
||||||
|
borderaxespad=0.0,
|
||||||
|
)
|
||||||
|
|
||||||
if output_dir is not None:
|
if output_dir is not None:
|
||||||
filename = pathjoin(output_dir, "%s_center.pdf" % target)
|
filename = pathjoin(output_dir, "%s_center.pdf" % target)
|
||||||
fig.savefig(filename, bbox_inches='tight', dpi=300, facecolor="None")
|
fig.savefig(filename, bbox_inches="tight", dpi=300, facecolor="None")
|
||||||
output.append(filename)
|
output.append(filename)
|
||||||
show()
|
show()
|
||||||
return output
|
return output
|
||||||
@@ -76,11 +103,57 @@ if __name__ == "__main__":
|
|||||||
import argparse
|
import argparse
|
||||||
|
|
||||||
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(
|
||||||
parser.add_argument("-f", "--file", metavar="path", required=False, help="The full or relative path to the data product", type=str, default=None)
|
"-t",
|
||||||
parser.add_argument("-c", "--pcut", metavar="pcut", required=False, help="The polarization cut for the data mask", type=float, default=0.99)
|
"--target",
|
||||||
parser.add_argument("-d", "--display", metavar="display", required=False, help="The map on which to display info", type=str, default="pf")
|
metavar="targetname",
|
||||||
parser.add_argument("-o", "--output_dir", metavar="directory_path", required=False, help="output directory path for the plots", type=str, default="./data")
|
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(
|
||||||
|
"-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",
|
||||||
|
)
|
||||||
args = parser.parse_args()
|
args = parser.parse_args()
|
||||||
exitcode = main(infile=args.file, P_cut=args.pcut, target=args.target, display=args.display, 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)
|
||||||
|
|||||||
Reference in New Issue
Block a user