diff --git a/README.md b/README.md index 1731062..3fe06c9 100755 --- a/README.md +++ b/README.md @@ -1,2 +1,6 @@ # FOC_Reduction FOC reduction pipeline + +TODO: + - Add all polarimetry capables instruments from HST (starting with FOS and ACS) + - Build science case for future UV polarimeters (AGN outflow geometry, dust scattering, torus outline ?) diff --git a/doc/pipeline.pdf b/doc/pipeline.pdf new file mode 100644 index 0000000..d6e7e05 Binary files /dev/null and b/doc/pipeline.pdf differ diff --git a/package/FOC_reduction.py b/package/FOC_reduction.py index e6920d0..42a7bab 100755 --- a/package/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -41,8 +41,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Background estimation error_sub_type = "freedman-diaconis" # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) - subtract_error = 0.5 - display_bkg = False + subtract_error = 1.0 + display_bkg = True # Data binning pxsize = 0.05 @@ -51,7 +51,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Alignement align_center = "center" # If None will not align the images - display_align = False + display_align = True display_data = False # Transmittance correction @@ -176,6 +176,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= proj_plots.plot_obs( data_array, headers, + shifts=shifts, savename="_".join([figname, str(align_center)]), plots_folder=plots_folder, norm=LogNorm(vmin=data_array[data_array > 0.0].min() * headers[0]["photflam"], vmax=data_array[data_array > 0.0].max() * headers[0]["photflam"]), @@ -305,114 +306,25 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= savename="_".join([figname]), plots_folder=plots_folder, ) - proj_plots.polarization_map( - deepcopy(Stokes_hdul), - data_mask, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - flux_lim=flux_lim, - step_vec=step_vec, - scale_vec=scale_vec, - savename="_".join([figname, "I"]), - plots_folder=plots_folder, - display="Intensity", - ) - proj_plots.polarization_map( - deepcopy(Stokes_hdul), - data_mask, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - flux_lim=flux_lim, - step_vec=step_vec, - scale_vec=scale_vec, - savename="_".join([figname, "P_flux"]), - plots_folder=plots_folder, - display="Pol_Flux", - ) - proj_plots.polarization_map( - deepcopy(Stokes_hdul), - data_mask, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - flux_lim=flux_lim, - step_vec=step_vec, - scale_vec=scale_vec, - savename="_".join([figname, "P"]), - plots_folder=plots_folder, - display="Pol_deg", - ) - proj_plots.polarization_map( - deepcopy(Stokes_hdul), - data_mask, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - flux_lim=flux_lim, - step_vec=step_vec, - scale_vec=scale_vec, - savename="_".join([figname, "PA"]), - plots_folder=plots_folder, - display="Pol_ang", - ) - proj_plots.polarization_map( - deepcopy(Stokes_hdul), - data_mask, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - flux_lim=flux_lim, - step_vec=step_vec, - scale_vec=scale_vec, - savename="_".join([figname, "I_err"]), - plots_folder=plots_folder, - display="I_err", - ) - proj_plots.polarization_map( - deepcopy(Stokes_hdul), - data_mask, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - flux_lim=flux_lim, - step_vec=step_vec, - scale_vec=scale_vec, - savename="_".join([figname, "P_err"]), - plots_folder=plots_folder, - display="Pol_deg_err", - ) - proj_plots.polarization_map( - deepcopy(Stokes_hdul), - data_mask, - P_cut=P_cut, - SNRi_cut=SNRi_cut, - flux_lim=flux_lim, - step_vec=step_vec, - scale_vec=scale_vec, - savename="_".join([figname, "SNRi"]), - plots_folder=plots_folder, - display="SNRi", - ) - proj_plots.polarization_map( - deepcopy(Stokes_hdul), - data_mask, - P_cut=P_cut if P_cut >= 1.0 else 3.0, - SNRi_cut=SNRi_cut, - flux_lim=flux_lim, - step_vec=step_vec, - scale_vec=scale_vec, - savename="_".join([figname, "SNRp"]), - plots_folder=plots_folder, - display="SNRp", - ) - proj_plots.polarization_map( - deepcopy(Stokes_hdul), - data_mask, - P_cut=P_cut if P_cut < 1.0 else 0.99, - SNRi_cut=SNRi_cut, - flux_lim=flux_lim, - step_vec=step_vec, - scale_vec=scale_vec, - savename="_".join([figname, "confP"]), - plots_folder=plots_folder, - display="confp", - ) + for figtype, figsuffix in zip( + ["Intensity", "Pol_flux", "Pol_deg", "Pol_ang", "I_err", "P_err", "SNRi", "SNRp", "confp"], + ["I", "P_flux", "P", "PA", "I_err", "P_err", "SNRi", "SNRp", "confP"], + ): + try: + proj_plots.polarization_map( + deepcopy(Stokes_hdul), + data_mask, + P_cut=P_cut, + SNRi_cut=SNRi_cut, + flux_lim=flux_lim, + step_vec=step_vec, + scale_vec=scale_vec, + savename="_".join([figname, figsuffix]), + plots_folder=plots_folder, + display=figtype, + ) + except ValueError: + pass elif not interactive: proj_plots.polarization_map( deepcopy(Stokes_hdul), data_mask, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=figname, plots_folder=plots_folder, display="integrate" diff --git a/package/lib/plots.py b/package/lib/plots.py index d4e49df..12e3b2e 100755 --- a/package/lib/plots.py +++ b/package/lib/plots.py @@ -66,7 +66,7 @@ except ImportError: from utils import PCconf, princ_angle, rot2D, sci_not -def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="", **kwargs): +def plot_obs(data_array, headers, rectangle=None, shifts=None, savename=None, plots_folder="", **kwargs): """ Plots raw observation imagery with some information on the instrument and filters. @@ -77,16 +77,14 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" single observation with multiple polarizers of an instrument headers : header list List of headers corresponding to the images in data_array - vmin : float, optional - Min pixel value that should be displayed. - Defaults to 0. - vmax : float, optional - Max pixel value that should be displayed. - Defaults to 6. rectangle : numpy.ndarray, optional Array of parameters for matplotlib.patches.Rectangle objects that will be displayed on each output image. If None, no rectangle displayed. Defaults to None. + shifts : numpy.ndarray, optional + Array of vector coordinates corresponding to images shifts with respect + to the others. If None, no shifts displayed. + Defaults to None. savename : str, optional Name of the figure the map should be saved to. If None, the map won't be saved (only displayed). @@ -100,6 +98,7 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" nb_obs = np.max([np.sum([head["filtnam1"] == curr_pol for head in headers]) for curr_pol in ["POL0", "POL60", "POL120"]]) shape = np.array((3, nb_obs)) fig, ax = plt.subplots(shape[0], shape[1], figsize=(3 * shape[1], 3 * shape[0]), layout="constrained", sharex=True, sharey=True) + used = np.zeros(shape, dtype=bool) r_pol = dict(pol0=0, pol60=1, pol120=2) c_pol = dict(pol0=0, pol60=0, pol120=0) for i, (data, head) in enumerate(zip(data_array, headers)): @@ -112,15 +111,17 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" c_pol[filt.lower()] += 1 if shape[1] != 1: ax_curr = ax[r_ax][c_ax] + used[r_ax][c_ax] = True else: ax_curr = ax[r_ax] + ax_curr[r_ax] = True # plots if "vmin" in kwargs.keys() or "vmax" in kwargs.keys(): vmin, vmax = kwargs["vmin"], kwargs["vmax"] del kwargs["vmin"], kwargs["vmax"] else: vmin, vmax = convert * data[data > 0.0].min() / 10.0, convert * data[data > 0.0].max() - for key, value in [["cmap", [["cmap", "gray"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: + for key, value in [["cmap", [["cmap", "inferno"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: try: _ = kwargs[key] except KeyError: @@ -129,17 +130,29 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" # im = ax[r_ax][c_ax].imshow(convert*data, origin='lower', **kwargs) data[data * convert < vmin * 10.0] = vmin * 10.0 / convert im = ax_curr.imshow(convert * data, origin="lower", **kwargs) + if shifts is not None: + x, y = np.array(data.shape[::-1]) / 2.0 - shifts[i] + dx, dy = shifts[i] + ax_curr.arrow(x, y, dx, dy, length_includes_head=True, width=0.1, head_width=0.3, color="g") + ax_curr.plot([x, x], [0, data.shape[0] - 1], "--", lw=2, color="g", alpha=0.85) + ax_curr.plot([0, data.shape[1] - 1], [y, y], "--", lw=2, color="g", alpha=0.85) if rectangle is not None: x, y, width, height, angle, color = rectangle[i] ax_curr.add_patch(Rectangle((x, y), width, height, angle=angle, edgecolor=color, fill=False)) # position of centroid - ax_curr.plot([data.shape[1] / 2, data.shape[1] / 2], [0, data.shape[0] - 1], "--", lw=1, color="grey", alpha=0.5) - ax_curr.plot([0, data.shape[1] - 1], [data.shape[1] / 2, data.shape[1] / 2], "--", lw=1, color="grey", alpha=0.5) + ax_curr.plot([data.shape[1] / 2, data.shape[1] / 2], [0, data.shape[0] - 1], "--", lw=2, color="b", alpha=0.85) + ax_curr.plot([0, data.shape[1] - 1], [data.shape[0] / 2, data.shape[0] / 2], "--", lw=2, color="b", alpha=0.85) ax_curr.annotate( - instr + ":" + rootname, color="white", fontsize=5, xy=(0.01, 1.00), xycoords="axes fraction", verticalalignment="top", horizontalalignment="left" + instr + ":" + rootname, color="white", fontsize=10, xy=(0.01, 1.00), xycoords="axes fraction", verticalalignment="top", horizontalalignment="left" ) - ax_curr.annotate(filt, color="white", fontsize=10, xy=(0.01, 0.01), xycoords="axes fraction", verticalalignment="bottom", horizontalalignment="left") - ax_curr.annotate(exptime, color="white", fontsize=5, xy=(1.00, 0.01), xycoords="axes fraction", verticalalignment="bottom", horizontalalignment="right") + ax_curr.annotate(filt, color="white", fontsize=15, xy=(0.01, 0.01), xycoords="axes fraction", verticalalignment="bottom", horizontalalignment="left") + ax_curr.annotate( + exptime, color="white", fontsize=10, xy=(1.00, 0.01), xycoords="axes fraction", verticalalignment="bottom", horizontalalignment="right" + ) + unused = np.logical_not(used) + ii, jj = np.indices(shape) + for i, j in zip(ii[unused], jj[unused]): + fig.delaxes(ax[i][j]) # fig.subplots_adjust(hspace=0.01, wspace=0.01, right=1.02) fig.colorbar(im, ax=ax, location="right", shrink=0.75, aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") @@ -349,7 +362,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.10, stkI.shape[0] * 1.15]) + ax.set(aspect="equal", fc="k", ylim=[-stkI.shape[0] * 0.01, stkI.shape[0] * 1.01]) # fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02) # ax.coords.grid(True, color='white', ls='dotted', alpha=0.5) @@ -677,7 +690,7 @@ class align_maps(object): except KeyError: for key_i, val_i in value: kwargs[key_i] = val_i - self.map_ax.imshow(self.map_data * self.map_convert, aspect="equal", **kwargs) + self.im = self.map_ax.imshow(self.map_data * self.map_convert, aspect="equal", **kwargs) if kwargs["cmap"] in [ "inferno", @@ -766,7 +779,7 @@ class align_maps(object): except KeyError: for key_i, val_i in value: other_kwargs[key_i] = val_i - self.other_ax.imshow(self.other_data * self.other_convert, aspect="equal", **other_kwargs) + self.other_im = self.other_ax.imshow(self.other_data * self.other_convert, aspect="equal", **other_kwargs) px_size2 = self.other_wcs.wcs.get_cdelt()[0] * 3600.0 px_sc2 = AnchoredSizeBar( @@ -1561,6 +1574,7 @@ class overplot_pol(align_maps): while not self.aligned: self.align() self.overplot(levels=levels, P_cut=P_cut, SNRi_cut=SNRi_cut, scale_vec=scale_vec, savename=savename, **kwargs) + plt.show(block=True) def add_vector(self, position="center", pol_deg=1.0, pol_ang=0.0, **kwargs): if isinstance(position, str) and position == "center": @@ -2148,7 +2162,7 @@ class image_lasso_selector(object): self.mask = np.zeros(self.img.shape[:2], dtype=bool) self.mask[self.indices] = True if hasattr(self, "cont"): - for coll in self.cont: + for coll in self.cont.collections: coll.remove() self.cont = self.ax.contour(self.mask.astype(float), levels=[0.5], colors="white", linewidths=1) if not self.embedded: @@ -2261,7 +2275,7 @@ class slit(object): for p in self.pix: self.mask[tuple(p)] = (np.abs(np.dot(rot2D(-self.angle), p - self.rect.get_center()[::-1])) < (self.height / 2.0, self.width / 2.0)).all() if hasattr(self, "cont"): - for coll in self.cont: + for coll in self.cont.collections: try: coll.remove() except AttributeError: @@ -2364,7 +2378,7 @@ class aperture(object): x0, y0 = self.circ.center self.mask = np.sqrt((xx - x0) ** 2 + (yy - y0) ** 2) < self.radius if hasattr(self, "cont"): - for coll in self.cont: + for coll in self.cont.collections: try: coll.remove() except AttributeError: @@ -2499,7 +2513,7 @@ class pol_map(object): self.selected = False self.region = deepcopy(self.select_instance.mask.astype(bool)) self.select_instance.displayed.remove() - for coll in self.select_instance.cont: + for coll in self.select_instance.cont.collections: coll.remove() self.select_instance.lasso.set_active(False) self.set_data_mask(deepcopy(self.region)) @@ -2543,7 +2557,7 @@ class pol_map(object): self.select_instance.update_mask() self.region = deepcopy(self.select_instance.mask.astype(bool)) self.select_instance.displayed.remove() - for coll in self.select_instance.cont: + for coll in self.select_instance.cont.collections: coll.remove() self.select_instance.circ.set_visible(False) self.set_data_mask(deepcopy(self.region)) @@ -2601,7 +2615,7 @@ class pol_map(object): self.select_instance.update_mask() self.region = deepcopy(self.select_instance.mask.astype(bool)) self.select_instance.displayed.remove() - for coll in self.select_instance.cont: + for coll in self.select_instance.cont.collections: coll.remove() self.select_instance.rect.set_visible(False) self.set_data_mask(deepcopy(self.region)) @@ -3324,7 +3338,7 @@ class pol_map(object): ) if hasattr(self, "cont"): - for coll in self.cont: + for coll in self.cont.collections: try: coll.remove() except AttributeError: diff --git a/package/src/Combine.py b/package/src/Combine.py index 87b9632..59bc2ce 100755 --- a/package/src/Combine.py +++ b/package/src/Combine.py @@ -173,7 +173,7 @@ def main(infiles, target=None, output_dir="./data/"): # Reduction parameters kwargs = {} # Polarization map output - kwargs["SNRp_cut"] = 3.0 + kwargs["P_cut"] = 0.99 kwargs["SNRi_cut"] = 1.0 kwargs["flux_lim"] = 1e-19, 3e-17 kwargs["scale_vec"] = 5 @@ -186,9 +186,7 @@ def main(infiles, target=None, output_dir="./data/"): new_infiles = [] for i, group in enumerate(grouped_infiles): - new_infiles.append( - FOC_reduction(target=target + "-" + str(i + 1), infiles=["/".join([data_folder, file]) for file in group], interactive=True)[0] - ) + new_infiles.append(FOC_reduction(target=target + "-" + str(i + 1), infiles=["/".join([data_folder, file]) for file in group], interactive=True)[0]) infiles = new_infiles diff --git a/package/src/overplot_IC5063.py b/package/src/overplot_IC5063.py index f9ac6c3..1d58c88 100755 --- a/package/src/overplot_IC5063.py +++ b/package/src/overplot_IC5063.py @@ -11,46 +11,46 @@ from lib.plots import overplot_pol, overplot_radio from matplotlib.colors import LogNorm Stokes_UV = fits.open("./data/IC5063/5918/IC5063_FOC_b0.10arcsec_c0.20arcsec.fits") -Stokes_18GHz = fits.open("./data/IC5063/radio/IC5063_18GHz.fits") -Stokes_24GHz = fits.open("./data/IC5063/radio/IC5063_24GHz.fits") -Stokes_103GHz = fits.open("./data/IC5063/radio/IC5063_103GHz.fits") -Stokes_229GHz = fits.open("./data/IC5063/radio/IC5063_229GHz.fits") -Stokes_357GHz = fits.open("./data/IC5063/radio/IC5063_357GHz.fits") +# Stokes_18GHz = fits.open("./data/IC5063/radio/IC5063_18GHz.fits") +# Stokes_24GHz = fits.open("./data/IC5063/radio/IC5063_24GHz.fits") +# Stokes_103GHz = fits.open("./data/IC5063/radio/IC5063_103GHz.fits") +# Stokes_229GHz = fits.open("./data/IC5063/radio/IC5063_229GHz.fits") +# Stokes_357GHz = fits.open("./data/IC5063/radio/IC5063_357GHz.fits") # Stokes_S2 = fits.open("./data/IC5063/POLARIZATION_COMPARISON/S2_rot_crop.fits") Stokes_IR = fits.open("./data/IC5063/IR/u2e65g01t_c0f_rot.fits") # levelsMorganti = np.array([1.,2.,3.,8.,16.,32.,64.,128.]) -levelsMorganti = np.logspace(-0.1249, 1.97, 7) / 100.0 +# levelsMorganti = np.logspace(-0.1249, 1.97, 7) / 100.0 -levels18GHz = levelsMorganti * Stokes_18GHz[0].data.max() -A = overplot_radio(Stokes_UV, Stokes_18GHz) -A.plot(levels=levels18GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/18GHz_overplot.pdf", scale_vec=None) +# levels18GHz = levelsMorganti * Stokes_18GHz[0].data.max() +# A = overplot_radio(Stokes_UV, Stokes_18GHz) +# A.plot(levels=levels18GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/18GHz_overplot.pdf", scale_vec=None) -levels24GHz = levelsMorganti * Stokes_24GHz[0].data.max() -B = overplot_radio(Stokes_UV, Stokes_24GHz) -B.plot(levels=levels24GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/24GHz_overplot.pdf", scale_vec=None) +# levels24GHz = levelsMorganti * Stokes_24GHz[0].data.max() +# B = overplot_radio(Stokes_UV, Stokes_24GHz) +# B.plot(levels=levels24GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/24GHz_overplot.pdf", scale_vec=None) -levels103GHz = levelsMorganti * Stokes_103GHz[0].data.max() -C = overplot_radio(Stokes_UV, Stokes_103GHz) -C.plot(levels=levels103GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/103GHz_overplot.pdf", scale_vec=None) +# levels103GHz = levelsMorganti * Stokes_103GHz[0].data.max() +# C = overplot_radio(Stokes_UV, Stokes_103GHz) +# C.plot(levels=levels103GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/103GHz_overplot.pdf", scale_vec=None) -levels229GHz = levelsMorganti * Stokes_229GHz[0].data.max() -D = overplot_radio(Stokes_UV, Stokes_229GHz) -D.plot(levels=levels229GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/229GHz_overplot.pdf", scale_vec=None) +# levels229GHz = levelsMorganti * Stokes_229GHz[0].data.max() +# D = overplot_radio(Stokes_UV, Stokes_229GHz) +# D.plot(levels=levels229GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/229GHz_overplot.pdf", scale_vec=None) -levels357GHz = levelsMorganti * Stokes_357GHz[0].data.max() -E = overplot_radio(Stokes_UV, Stokes_357GHz) -E.plot(levels=levels357GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/357GHz_overplot.pdf", scale_vec=None) +# levels357GHz = levelsMorganti * Stokes_357GHz[0].data.max() +# E = overplot_radio(Stokes_UV, Stokes_357GHz) +# E.plot(levels=levels357GHz, P_cut=2.0, SNRi_cut=10.0, savename="./plots/IC5063/357GHz_overplot.pdf", scale_vec=None) # F = overplot_pol(Stokes_UV, Stokes_S2) # F.plot(P_cut=3.0, SNRi_cut=80.0, savename='./plots/IC5063/S2_overplot.pdf', norm=LogNorm(vmin=5e-20,vmax=5e-18)) G = overplot_pol(Stokes_UV, Stokes_IR, cmap="inferno") G.plot( - P_cut=2.0, - SNRi_cut=10.0, + P_cut=0.99, + SNRi_cut=1.0, savename="./plots/IC5063/IR_overplot.pdf", scale_vec=None, norm=LogNorm(Stokes_IR[0].data.max() * Stokes_IR[0].header["photflam"] / 1e3, Stokes_IR[0].data.max() * Stokes_IR[0].header["photflam"]), - cmap="inferno_r", + cmap="inferno", ) diff --git a/package/src/overplot_MRK463E.py b/package/src/overplot_MRK463E.py index a05bfb0..d89845e 100755 --- a/package/src/overplot_MRK463E.py +++ b/package/src/overplot_MRK463E.py @@ -11,16 +11,33 @@ from lib.plots import overplot_chandra, overplot_pol from matplotlib.colors import LogNorm Stokes_UV = fits.open("./data/MRK463E/5960/MRK463E_FOC_b0.05arcsec_c0.10arcsec.fits") -Stokes_IR = fits.open("./data/MRK463E/WFPC2/IR_rot_crop.fits") -Stokes_Xr = fits.open("./data/MRK463E/Chandra/X_ray_crop.fits") +# Stokes_IR = fits.open("./data/MRK463E/WFPC2/IR_rot_crop.fits") +# Stokes_Xr = fits.open("./data/MRK463E/Chandra/X_ray_crop.fits") +Radio = fits.open("./data/MRK463E/EMERLIN/Voorwerpjes_1356+1822_1356+1822_uniform-image.fits") -levels = np.geomspace(1.0, 99.0, 7) +# levels = np.geomspace(1.0, 99.0, 7) -A = overplot_chandra(Stokes_UV, Stokes_Xr, norm=LogNorm()) -A.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, zoom=1, savename="./plots/MRK463E/Chandra_overplot.pdf") -A.write_to(path1="./data/MRK463E/FOC_data_Chandra.fits", path2="./data/MRK463E/Chandra_data.fits", suffix="aligned") +# A = overplot_chandra(Stokes_UV, Stokes_Xr, norm=LogNorm()) +# A.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, zoom=1, savename="./plots/MRK463E/Chandra_overplot.pdf") +# A.write_to(path1="./data/MRK463E/FOC_data_Chandra.fits", path2="./data/MRK463E/Chandra_data.fits", suffix="aligned") + +# levels = np.array([0.8, 2, 5, 10, 20, 50]) / 100.0 * Stokes_UV[0].header["photflam"] +# B = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm()) +# B.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, norm=LogNorm(8.5e-18, 2.5e-15), savename="./plots/MRK463E/IR_overplot.pdf") +# B.write_to(path1="./data/MRK463E/FOC_data_WFPC.fits", path2="./data/MRK463E/WFPC_data.fits", suffix="aligned") levels = np.array([0.8, 2, 5, 10, 20, 50]) / 100.0 * Stokes_UV[0].header["photflam"] -B = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm()) -B.plot(levels=levels, P_cut=0.99, SNRi_cut=1.0, scale_vec=5, norm=LogNorm(8.5e-18, 2.5e-15), savename="./plots/MRK463E/IR_overplot.pdf") -B.write_to(path1="./data/MRK463E/FOC_data_WFPC.fits", path2="./data/MRK463E/WFPC_data.fits", suffix="aligned") +C = overplot_pol(Stokes_UV, Radio, norm=LogNorm()) +C.other_im.set(norm=LogNorm(1e-4, 2e-2)) +C.plot( + levels=levels, + P_cut=0.999, + SNRi_cut=5.0, + scale_vec=3, + norm=LogNorm(1e-4, 2e-2), + cmap="inferno_r", + width=0.5, + linewidth=0.5, + savename="./plots/MRK463E/EMERLIN_overplot.pdf", +) +C.write_to(path1="./data/MRK463E/FOC_data_EMERLIN.fits", path2="./data/MRK463E/EMERLIN_data.fits", suffix="aligned") diff --git a/package/src/overplot_NGC1068.py b/package/src/overplot_NGC1068.py index fa4f31f..52ca188 100755 --- a/package/src/overplot_NGC1068.py +++ b/package/src/overplot_NGC1068.py @@ -30,10 +30,10 @@ A.add_vector( A.other_wcs.celestial.wcs.crpix - (1.0, 1.0), pol_deg=0.124, pol_ang=100.7, - width=2, - linewidth=2, + width=2.0, + linewidth=1.0, scale=1.0 / (A.px_scale * 6.0), - edgecolor="b", + edgecolor="w", color="b", label=r"IXPE torus: P = $12.4 (\pm 3.6)$%, PA = $100.7 (\pm 8.3)$°", )