diff --git a/package/lib/fits.py b/package/lib/fits.py index 6ea842a..528e602 100755 --- a/package/lib/fits.py +++ b/package/lib/fits.py @@ -46,6 +46,13 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): data_array.append(f[0].data) wcs_array.append(WCS(header=f[0].header, fobj=f).celestial) f.flush() + # Save pixel area for flux density computation + if headers[i]["PXFORMT"] == "NORMAL": + headers[i]["PXAREA"] = 6.25e-6 # 25x25 micron squared pixel area in cm^2 + elif headers[i]["PXFORMT"] == "ZOOM": + headers[i]["PXAREA"] = 1.25e-5 # 50x25 micron squared pixel area in cm^2 + else: + headers[i]["PXAREA"] = 1.0 # unknown default to 1 cm^2 data_array = np.array(data_array, dtype=np.double) # Prevent negative count value in imported data @@ -143,6 +150,7 @@ def save_Stokes( header["PHOTPLAM"] = (header_stokes["PHOTPLAM"], "Pivot Wavelength") header["PHOTBW"] = (header_stokes["PHOTBW"], "RMS Bandwidth of the Filter and Detector") header["PHOTFLAM"] = (header_stokes["PHOTFLAM"], "Inverse Sensitivity in DN/sec/cm**2/Angst") + header["PXAREA"] = (header_stokes["PXAREA"], "Pixel area in cm**2") header["EXPTIME"] = (header_stokes["EXPTIME"], "Total exposure time in sec") header["PROPOSID"] = (header_stokes["PROPOSID"], "PEP proposal identifier for observation") header["TARGNAME"] = (header_stokes["TARGNAME"], "Target name") diff --git a/package/lib/plots.py b/package/lib/plots.py index 8120809..0c9a201 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -114,7 +114,7 @@ def plot_obs(data_array, headers, rectangle=None, shifts=None, savename=None, pl used[r_ax][c_ax] = True else: ax_curr = ax[r_ax] - ax_curr[r_ax] = True + used[r_ax] = True # plots if "vmin" in kwargs.keys() or "vmax" in kwargs.keys(): vmin, vmax = kwargs["vmin"], kwargs["vmax"] @@ -338,9 +338,12 @@ def polarization_map( if P_cut >= 1.0: # MaskP on the signal-to-noise value - maskP[pol_err > 0.0] = SNRp[pol_err > 0.0] > P_cut + SNRp_cut = deepcopy(P_cut) + maskP[pol_err > 0.0] = SNRp[pol_err > 0.0] > SNRp_cut + P_cut = 0.99 else: # MaskP on the confidence value + SNRp_cut = 5.0 maskP = confP > P_cut mask = np.logical_and(maskI, maskP) @@ -363,7 +366,7 @@ def polarization_map( fig = plt.figure(figsize=(7 * ratiox, 7 * ratioy), layout="constrained") if ax is None: ax = fig.add_subplot(111, projection=wcs) - ax.set(aspect="equal", fc="k", ylim=[-stkI.shape[0] * 0.01, stkI.shape[0] * 1.01]) + ax.set(aspect="equal", fc="k", ylim=[-stkI.shape[0] * 0.00, stkI.shape[0] * 1.00]) # fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02) # ax.coords.grid(True, color='white', ls='dotted', alpha=0.5) @@ -375,8 +378,8 @@ def polarization_map( vmin, vmax = 0.0, stkI.max() * convert_flux for key, value in [ ["cmap", [["cmap", "inferno"]]], - ["width", [["width", 0.5]]], - ["linewidth", [["linewidth", 0.3]]], + ["width", [["width", 0.4]]], + ["linewidth", [["linewidth", 0.6]]], ]: try: _ = kwargs[key] @@ -424,7 +427,7 @@ def polarization_map( else: 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.75, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + 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}$]") 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) @@ -440,8 +443,9 @@ def polarization_map( pfmax = (stkI[stkI > 0.0] * pol[stkI > 0.0] * convert_flux).max() else: 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=50, shrink=0.75, pad=0.025, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + 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 levelsPf = np.array([1.73, 13.0, 33.0, 66.0]) / 100.0 * pfmax print("Polarized flux density contour levels : ", levelsPf) @@ -451,13 +455,13 @@ def polarization_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=50, shrink=0.75, pad=0.025, label=r"$P$ [%]") + fig.colorbar(im, ax=ax, aspect=30, 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=50, shrink=0.75, pad=0.025, label=r"$\theta_P$ [°]") + fig.colorbar(im, ax=ax, aspect=30, 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" @@ -467,7 +471,7 @@ def polarization_map( 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=50, shrink=0.75, pad=0.025, label=r"$\sigma_P$ [%]") + fig.colorbar(im, ax=ax, aspect=30, 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" @@ -479,7 +483,7 @@ def polarization_map( 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=50, shrink=0.75, pad=0.025, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") + 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}$]") elif display.lower() in ["snri"]: # Display I_stokes signal-to-noise map display = "snri" @@ -491,19 +495,19 @@ def polarization_map( 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=50, shrink=0.75, pad=0.025, label=r"$I_{Stokes}/\sigma_{I}$") + fig.colorbar(im, ax=ax, aspect=30, 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 > P_cut: + if vmax * 0.99 > SNRp_cut: 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=50, shrink=0.75, pad=0.025, label=r"$P/\sigma_{P}$") + fig.colorbar(im, ax=ax, aspect=30, 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" @@ -512,7 +516,7 @@ 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=50, shrink=0.75, pad=0.025, label=r"$Conf_{P}$") + fig.colorbar(im, ax=ax, aspect=30, shrink=0.60, pad=0.025, label=r"$Conf_{P}$") else: # Defaults to intensity map if mask.sum() > 0.0: @@ -520,18 +524,20 @@ def polarization_map( 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=50, shrink=0.75, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") + 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$]") - # Get integrated values from header - I_diluted = stkI[data_mask].sum() - I_diluted_err = np.sqrt(np.sum(stk_cov[0, 0][data_mask])) + # Get integrated flux values from sum + N_pix = data_mask.sum() + I_diluted = stkI[data_mask].sum() / N_pix + I_diluted_err = np.sqrt(np.sum(stk_cov[0, 0][data_mask])) / N_pix + # Get integrated polarization values from header P_diluted = Stokes[0].header["P_int"] P_diluted_err = Stokes[0].header["sP_int"] PA_diluted = Stokes[0].header["PA_int"] PA_diluted_err = Stokes[0].header["sPA_int"] - plt.rcParams.update({"font.size": 10}) + plt.rcParams.update({"font.size": 11}) px_size = wcs.wcs.get_cdelt()[0] * 3600.0 px_sc = AnchoredSizeBar(ax.transData, 1.0 / px_size, "1 arcsec", 3, pad=0.25, sep=5, borderpad=0.25, frameon=False, size_vertical=0.005, color=font_color) north_dir = AnchoredDirectionArrows( @@ -548,7 +554,7 @@ def polarization_map( head_length=10.0, head_width=10.0, angle=-Stokes[0].header["orientat"], - text_props={"ec": "k", "fc": font_color, "alpha": 1, "lw": 0.4}, + text_props={"ec": "k", "fc": font_color, "alpha": 1, "lw": 0.5}, arrow_props={"ec": "k", "fc": font_color, "alpha": 1, "lw": 1}, ) @@ -1930,7 +1936,7 @@ class crop_map(object): self.im.remove() self.im = self.ax.imshow(data * convert_flux, **kwargs) if hasattr(self, "cr"): - self.cr[0].set_data(*wcs.wcs.crpix) + self.cr[0].set_data([wcs.wcs.crpix[0]], [wcs.wcs.crpix[1]]) else: self.cr = self.ax.plot(*wcs.wcs.crpix, "r+") self.fig.canvas.draw_idle() @@ -3275,7 +3281,7 @@ class pol_map(object): str_conf = "" if self.region is None: s_I = np.sqrt(self.IQU_cov[0, 0]) - N_pix = self.I.size + N_pix = self.I.size() I_reg = self.I.sum() I_reg_err = np.sqrt(np.sum(s_I**2)) P_reg = self.Stokes[0].header["P_int"] @@ -3398,7 +3404,7 @@ class pol_map(object): self.an_int.remove() self.str_int = ( r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( - self.pivot_wav, sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2) + self.pivot_wav, sci_not(I_reg * self.map_convert / N_pix, I_reg_err * self.map_convert / N_pix, 2) ) + "\n" + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg * 100.0, np.ceil(P_reg_err * 1000.0) / 10.0) @@ -3410,7 +3416,7 @@ class pol_map(object): # self.str_cut = ( # "\n" # + r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( - # self.pivot_wav, sci_not(I_cut * self.map_convert, I_cut_err * self.map_convert, 2) + # self.pivot_wav, sci_not(I_cut * self.map_convert/N_cut, I_cut_err * self.map_convert/N_cut, 2) # ) # + "\n" # + r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut * 100.0, np.ceil(P_cut_err * 1000.0) / 10.0) @@ -3434,7 +3440,7 @@ class pol_map(object): else: str_int = ( r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( - self.pivot_wav, sci_not(I_reg * self.map_convert, I_reg_err * self.map_convert, 2) + self.pivot_wav, sci_not(I_reg * self.map_convert / N_pix, I_reg_err * self.map_convert / N_pix, 2) ) + "\n" + r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg * 100.0, np.ceil(P_reg_err * 1000.0) / 10.0) @@ -3446,7 +3452,7 @@ class pol_map(object): # str_cut = ( # "\n" # + r"$F_{{\lambda}}^{{cut}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format( - # self.pivot_wav, sci_not(I_cut * self.map_convert, I_cut_err * self.map_convert, 2) + # self.pivot_wav, sci_not(I_cut * self.map_convert/N_cut, I_cut_err * self.map_convert/N_cut, 2) # ) # + "\n" # + r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut * 100.0, np.ceil(P_cut_err * 1000.0) / 10.0) diff --git a/package/lib/reduction.py b/package/lib/reduction.py index 613e0d7..c62168a 100755 --- a/package/lib/reduction.py +++ b/package/lib/reduction.py @@ -676,7 +676,9 @@ def rebin_array(data_array, error_array, headers, pxsize=2, scale="px", operatio nw.wcs.crpix /= Dxy nw.array_shape = new_shape new_header["NAXIS1"], new_header["NAXIS2"] = nw.array_shape - new_header["PHOTFLAM"] = header["PHOTFLAM"] * (Dxy[0] * Dxy[1]) + # Compute new values of pixel area and flux density conversion factor + new_header["PXAREA"] = header["PXAREA"] * (Dxy[0] * Dxy[1]) + new_header["PHOTFLAM"] = header["PHOTFLAM"] / (Dxy[0] * Dxy[1]) for key, val in nw.to_header().items(): new_header.set(key, val) new_header["SAMPLING"] = (str(pxsize) + scale, "Resampling performed during reduction")