diff --git a/package/lib/background.py b/package/lib/background.py index ed9638c..d120ae6 100755 --- a/package/lib/background.py +++ b/package/lib/background.py @@ -86,20 +86,19 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non ) ax_h.plot([background[i] * convert_flux[i], background[i] * convert_flux[i]], [hist.min(), hist.max()], "x--", color="C{0:d}".format(i), alpha=0.8) if i == 0: - xmin, xmax = np.min(np.array(bins)[np.array(hist) > 1e2]) * convert_flux[0], np.max(np.array(bins)[np.array(hist) > 1e2]) * convert_flux[0] + xmin, xmax = np.min(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0], np.max(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0] else: xmin, xmax = ( - min(xmin, np.min(np.array(bins)[np.array(hist) > 1e2]) * convert_flux[0]), - max(xmax, np.max(np.array(bins)[np.array(hist) > 1e2]) * convert_flux[0]), + min(xmin, np.min(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0]), + max(xmax, np.max(np.array(bins)[np.array(hist) > 5e1]) * convert_flux[0]), ) if coeff is not None: # ax_h.plot(bins*convert_flux[i], gausspol(bins, *coeff[i]), '--', color="C{0:d}".format(i), alpha=0.8) ax_h.plot(bins * convert_flux[i], gauss(bins, *coeff[i]), "--", color="C{0:d}".format(i), alpha=0.8) ax_h.set_xscale("log") ax_h.set_yscale("log") - ax_h.set_ylim([100.0, np.max([hist.max() for hist in histograms])]) - # ax_h.set_xlim([np.min(background * convert_flux) * 1e-2, np.max(background * convert_flux) * 1e2]) - ax_h.set_xlim([xmin, xmax]) + ax_h.set_ylim([5e1, np.max([hist.max() for hist in histograms])]) + ax_h.set_xlim([max(xmin, np.min(background * convert_flux) * 1e-2), min(xmax, np.max(background * convert_flux) * 1e2)]) ax_h.set_xlabel(r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") ax_h.set_ylabel(r"Number of pixels in bin") ax_h.set_title("Histogram for each observation") @@ -135,7 +134,7 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non ax2.set(xlabel="pixel offset", ylabel="pixel offset", aspect="equal") fig2.subplots_adjust(hspace=0, wspace=0, right=1.0) - fig2.colorbar(im2, ax=ax2, location="right", aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + fig2.colorbar(im2, ax=ax2, location="right", shrink=0.60, aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") if savename is not None: this_savename = deepcopy(savename) diff --git a/package/lib/plots.py b/package/lib/plots.py index 604d58f..a2278a9 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -428,7 +428,7 @@ def polarization_map( P_cut = 0.99 else: # MaskP on the confidence value - SNRp_cut = 5.0 + SNRp_cut = 3.0 maskP = confP > P_cut mask = np.logical_and(maskI, maskP) @@ -504,72 +504,32 @@ def polarization_map( if display.lower() in ["i", "intensity"]: # If no display selected, show intensity map display = "i" - if flux_lim is None: - if mask.sum() > 0.0: - vmin, vmax = ( - 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) - else: - vmin, vmax = ( - 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) - else: + if flux_lim is not None: vmin, vmax = flux_lim - im = ax.imshow( - stkI * convert_flux, - norm=LogNorm(vmin, vmax), - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar( - im, - ax=ax, - aspect=30, - shrink=0.60, - pad=0.025, - label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]", - ) + elif mask.sum() > 0.0: + vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) + else: + vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) + im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") levelsI = np.array([0.8, 2.0, 5.0, 10.0, 20.0, 50.0]) / 100.0 * vmax print("Flux density contour levels : ", levelsI) ax.contour(stkI * convert_flux, levels=levelsI, colors="grey", linewidths=0.5) elif display.lower() in ["pf", "pol_flux"]: # Display polarization flux display = "pf" - if flux_lim is None: - if mask.sum() > 0.0: - vmin, vmax = ( - 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) - pfmax = (stkI[mask] * pol[mask] * convert_flux).max() - else: - vmin, vmax = ( - 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) - pfmax = (stkI[stkI > 0.0] * pol[stkI > 0.0] * convert_flux).max() - else: + if flux_lim is not None: vmin, vmax = flux_lim pfmax = (stkI[stkI > 0.0] * pol[stkI > 0.0] * convert_flux).max() - im = ax.imshow( - stkI * convert_flux * pol, - norm=LogNorm(vmin, vmax), - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar( - im, - ax=ax, - aspect=30, - shrink=0.60, - pad=0.025, - label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]", - ) - # levelsPf = np.linspace(0.0175, 0.50, 5) * pfmax + elif mask.sum() > 0.0: + vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) + pfmax = (stkI[mask] * pol[mask] * convert_flux).max() + else: + vmin, vmax = 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) + pfmax = (stkI[stkI > 0.0] * pol[stkI > 0.0] * convert_flux).max() + im = ax.imshow(stkI * convert_flux * pol, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + # levelsPf = np.linspace(0.0.60, 0.50, 5) * pfmax levelsPf = np.array([1.73, 13.0, 33.0, 66.0]) / 100.0 * pfmax print("Polarized flux density contour levels : ", levelsPf) ax.contour(stkI * convert_flux * pol, levels=levelsPf, colors="grey", linewidths=0.5) @@ -577,52 +537,24 @@ def polarization_map( # Display polarization degree map display = "p" vmin, vmax = 0.0, 100.0 - im = ax.imshow( - pol * 100.0, - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar(im, ax=ax, aspect=30, shrink=0.60, pad=0.025, label=r"$P$ [%]") + im = ax.imshow(pol * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$P$ [%]") elif display.lower() in ["pa", "pang", "pol_ang"]: # Display polarization degree map display = "pa" vmin, vmax = 0.0, 180.0 - im = ax.imshow( - princ_angle(pang), - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar(im, ax=ax, aspect=30, shrink=0.60, pad=0.025, label=r"$\theta_P$ [°]") + im = ax.imshow(princ_angle(pang), vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\theta_P$ [°]") elif display.lower() in ["s_p", "pol_err", "pol_deg_err"]: # Display polarization degree error map display = "s_p" if (SNRp > P_cut).any(): vmin, vmax = 0.0, np.max([pol_err[SNRp > P_cut].max(), 1.0]) * 100.0 - im = ax.imshow( - pol_err * 100.0, - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap="inferno_r", - alpha=1.0, - ) + im = ax.imshow(pol_err * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno_r", alpha=1.0) else: vmin, vmax = 0.0, 100.0 - im = ax.imshow( - pol_err * 100.0, - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap="inferno_r", - alpha=1.0, - ) - fig.colorbar(im, ax=ax, aspect=30, shrink=0.60, pad=0.025, label=r"$\sigma_P$ [%]") + im = ax.imshow(pol_err * 100.0, vmin=vmin, vmax=vmax, aspect="equal", cmap="inferno_r", alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\sigma_P$ [%]") elif display.lower() in ["s_i", "i_err"]: # Display intensity error map display = "s_i" @@ -631,73 +563,34 @@ def polarization_map( 1.0 / 2.0 * np.median(np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.0]) * convert_flux), np.max(np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.0]) * convert_flux), ) - im = ax.imshow( - np.sqrt(stk_cov[0, 0]) * convert_flux, - norm=LogNorm(vmin, vmax), - aspect="equal", - cmap="inferno_r", - alpha=1.0, - ) + im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap="inferno_r", alpha=1.0) else: - im = ax.imshow( - np.sqrt(stk_cov[0, 0]) * convert_flux, - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar( - im, - ax=ax, - aspect=30, - shrink=0.60, - pad=0.025, - label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]", - ) + im = ax.imshow(np.sqrt(stk_cov[0, 0]) * convert_flux, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") elif display.lower() in ["snri"]: # Display I_stokes signal-to-noise map display = "snri" vmin, vmax = 0.0, np.max(SNRi[np.isfinite(SNRi)]) if vmax * 0.99 > SNRi_cut: - im = ax.imshow( - SNRi, - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) + im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) levelsSNRi = np.linspace(SNRi_cut, vmax * 0.99, 5).astype(int) print("SNRi contour levels : ", levelsSNRi) ax.contour(SNRi, levels=levelsSNRi, colors="grey", linewidths=0.5) else: im = ax.imshow(SNRi, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar( - im, - ax=ax, - aspect=30, - shrink=0.60, - pad=0.025, - label=r"$I_{Stokes}/\sigma_{I}$", - ) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$I_{Stokes}/\sigma_{I}$") elif display.lower() in ["snr", "snrp"]: # Display polarization degree signal-to-noise map display = "snrp" vmin, vmax = 0.0, np.max(SNRp[np.isfinite(SNRp)]) if vmax * 0.99 > SNRp_cut: - im = ax.imshow( - SNRp, - vmin=vmin, - vmax=vmax, - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) + im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) levelsSNRp = np.linspace(P_cut, vmax * 0.99, 5).astype(int) print("SNRp contour levels : ", levelsSNRp) ax.contour(SNRp, levels=levelsSNRp, colors="grey", linewidths=0.5) else: im = ax.imshow(SNRp, aspect="equal", cmap=kwargs["cmap"], alpha=1.0) - fig.colorbar(im, ax=ax, aspect=30, shrink=0.60, pad=0.025, label=r"$P/\sigma_{P}$") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$P/\sigma_{P}$") elif display.lower() in ["conf", "confp"]: # Display polarization degree signal-to-noise map display = "confp" @@ -706,34 +599,17 @@ def polarization_map( levelsconfp = np.array([0.500, 0.900, 0.990, 0.999]) print("confp contour levels : ", levelsconfp) ax.contour(confP, levels=levelsconfp, colors="grey", linewidths=0.5) - fig.colorbar(im, ax=ax, aspect=30, shrink=0.60, pad=0.025, label=r"$Conf_{P}$") + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$Conf_{P}$") else: # Defaults to intensity map - if mask.sum() > 0.0: - vmin, vmax = ( - 1.0 * np.mean(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) + if flux_lim is not None: + vmin, vmax = flux_lim + elif mask.sum() > 0.0: + vmin, vmax = 1.0 * np.mean(np.sqrt(stk_cov[0, 0][mask]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) else: - vmin, vmax = ( - 1.0 * np.mean(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), - np.max(stkI[stkI > 0.0] * convert_flux), - ) - im = ax.imshow( - stkI * convert_flux, - norm=LogNorm(vmin, vmax), - aspect="equal", - cmap=kwargs["cmap"], - alpha=1.0, - ) - fig.colorbar( - im, - ax=ax, - aspect=30, - shrink=0.60, - pad=0.025, - label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]", - ) + vmin, vmax = 1.0 * np.mean(np.sqrt(stk_cov[0, 0][stkI > 0.0]) * convert_flux), np.max(stkI[stkI > 0.0] * convert_flux) + im = ax.imshow(stkI * convert_flux, norm=LogNorm(vmin, vmax), aspect="equal", cmap=kwargs["cmap"], alpha=1.0) + fig.colorbar(im, ax=ax, aspect=50, shrink=0.60, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") # Get integrated flux values from sum I_diluted = stkI[data_mask].sum() @@ -3148,12 +3024,7 @@ class pol_map(object): self.snr_conf = 0 b_snr_conf.label.set_text("Conf") self.s_P_cut = Slider( - self.ax_P_cut, - r"$SNR^{P}_{cut}$", - 1.0, - int(SNRp_max * 0.95), - valstep=1, - valinit=self.P_cut if P_cut >= 1.0 else 3, + self.ax_P_cut, r"$SNR^{P}_{cut}$", 1.0, max(int(SNRp_max * 0.95), 3), valstep=1, valinit=self.P_cut if P_cut >= 1.0 else 3 ) else: self.snr_conf = 1 diff --git a/package/lib/utils.py b/package/lib/utils.py index 214f0f2..86dcf8c 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="TNC") + result = minimize(chisq, c0, bounds=[(0, PA.shape[1]), (0.0, PA.shape[0])], method="trust-constr") if result.success: print("Center of emission found") else: diff --git a/package/src/emission_center.py b/package/src/emission_center.py index c793206..c9fe314 100755 --- a/package/src/emission_center.py +++ b/package/src/emission_center.py @@ -44,7 +44,7 @@ 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=3, 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"))) @@ -62,11 +62,11 @@ 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.0, -0.02, 1.0, 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, dpi=300, facecolor="None") + fig.savefig(filename, bbox_inches='tight', dpi=300, facecolor="None") output.append(filename) show() return output