diff --git a/package/lib/utils.py b/package/lib/utils.py index 86dcf8c..214f0f2 100755 --- a/package/lib/utils.py +++ b/package/lib/utils.py @@ -109,7 +109,7 @@ def CenterConf(mask, PA, sPA): conf[np.isfinite(PA)] = gammaincc(0.5, 0.5 * chi2[np.isfinite(PA)]) 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: print("Center of emission found") else: diff --git a/package/src/emission_center.py b/package/src/emission_center.py index c9fe314..097362c 100755 --- a/package/src/emission_center.py +++ b/package/src/emission_center.py @@ -23,7 +23,12 @@ def main(infile, P_cut=0.99, target=None, display="pf", output_dir=None): stkI = Stokes["I_STOKES"].data QN, UN, QN_ERR, UN_ERR = np.full((4, stkI.shape[0], stkI.shape[1]), np.nan) 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], ): 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"] 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="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) confcentcont = ax.contour(Stokescentconf, [0.01], colors="red") # 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])) # labels.append(r"$SNR_P \geq$ 4 contour") # 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: 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) show() return output @@ -76,11 +103,57 @@ if __name__ == "__main__": import argparse 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("-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( + "-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( + "-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() - 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)