Tibeuleu/FOC_Reduction
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

small fixes and improvments

Browse Source
This commit is contained in:
Thibault Barnouin
2024-10-17 17:05:35 +02:00
parent bd7cad46a1
commit f6d62bff73
8 changed files with 119 additions and 174 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
README.md
View File

@@ -1,2 +1,6 @@
# FOC_Reduction # FOC_Reduction
FOC reduction pipeline 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 ?)

BIN
doc/pipeline.pdf Normal file
View File

Binary file not shown.

134
package/FOC_reduction.py
View File

@@ -41,8 +41,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
# Background estimation # Background estimation
error_sub_type = "freedman-diaconis" # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) error_sub_type = "freedman-diaconis" # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51))
subtract_error = 0.5 subtract_error = 1.0
display_bkg = False display_bkg = True
# Data binning # Data binning
pxsize = 0.05 pxsize = 0.05
@@ -51,7 +51,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
# Alignement # Alignement
align_center = "center" # If None will not align the images align_center = "center" # If None will not align the images
display_align = False display_align = True
display_data = False display_data = False
# Transmittance correction # Transmittance correction
@@ -176,6 +176,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
proj_plots.plot_obs( proj_plots.plot_obs(
data_array, data_array,
headers, headers,
shifts=shifts,
savename="_".join([figname, str(align_center)]), savename="_".join([figname, str(align_center)]),
plots_folder=plots_folder, 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"]), 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]), savename="_".join([figname]),
plots_folder=plots_folder, plots_folder=plots_folder,
) )
proj_plots.polarization_map( for figtype, figsuffix in zip(
deepcopy(Stokes_hdul), ["Intensity", "Pol_flux", "Pol_deg", "Pol_ang", "I_err", "P_err", "SNRi", "SNRp", "confp"],
data_mask, ["I", "P_flux", "P", "PA", "I_err", "P_err", "SNRi", "SNRp", "confP"],
P_cut=P_cut, ):
SNRi_cut=SNRi_cut, try:
flux_lim=flux_lim, proj_plots.polarization_map(
step_vec=step_vec, deepcopy(Stokes_hdul),
scale_vec=scale_vec, data_mask,
savename="_".join([figname, "I"]), P_cut=P_cut,
plots_folder=plots_folder, SNRi_cut=SNRi_cut,
display="Intensity", flux_lim=flux_lim,
) step_vec=step_vec,
proj_plots.polarization_map( scale_vec=scale_vec,
deepcopy(Stokes_hdul), savename="_".join([figname, figsuffix]),
data_mask, plots_folder=plots_folder,
P_cut=P_cut, display=figtype,
SNRi_cut=SNRi_cut, )
flux_lim=flux_lim, except ValueError:
step_vec=step_vec, pass
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",
)
elif not interactive: elif not interactive:
proj_plots.polarization_map( 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" deepcopy(Stokes_hdul), data_mask, P_cut=P_cut, SNRi_cut=SNRi_cut, savename=figname, plots_folder=plots_folder, display="integrate"

60
package/lib/plots.py
View File

@@ -66,7 +66,7 @@ except ImportError:
from utils import PCconf, princ_angle, rot2D, sci_not 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 Plots raw observation imagery with some information on the instrument and
filters. 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 single observation with multiple polarizers of an instrument
headers : header list headers : header list
List of headers corresponding to the images in data_array 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 rectangle : numpy.ndarray, optional
Array of parameters for matplotlib.patches.Rectangle objects that will Array of parameters for matplotlib.patches.Rectangle objects that will
be displayed on each output image. If None, no rectangle displayed. be displayed on each output image. If None, no rectangle displayed.
Defaults to None. 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 savename : str, optional
Name of the figure the map should be saved to. If None, the map won't Name of the figure the map should be saved to. If None, the map won't
be saved (only displayed). 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"]]) 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)) 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) 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) r_pol = dict(pol0=0, pol60=1, pol120=2)
c_pol = dict(pol0=0, pol60=0, pol120=0) c_pol = dict(pol0=0, pol60=0, pol120=0)
for i, (data, head) in enumerate(zip(data_array, headers)): 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 c_pol[filt.lower()] += 1
if shape[1] != 1: if shape[1] != 1:
ax_curr = ax[r_ax][c_ax] ax_curr = ax[r_ax][c_ax]
used[r_ax][c_ax] = True
else: else:
ax_curr = ax[r_ax] ax_curr = ax[r_ax]
ax_curr[r_ax] = True
# plots # plots
if "vmin" in kwargs.keys() or "vmax" in kwargs.keys(): if "vmin" in kwargs.keys() or "vmax" in kwargs.keys():
vmin, vmax = kwargs["vmin"], kwargs["vmax"] vmin, vmax = kwargs["vmin"], kwargs["vmax"]
del kwargs["vmin"], kwargs["vmax"] del kwargs["vmin"], kwargs["vmax"]
else: else:
vmin, vmax = convert * data[data > 0.0].min() / 10.0, convert * data[data > 0.0].max() 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: try:
_ = kwargs[key] _ = kwargs[key]
except KeyError: 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) # im = ax[r_ax][c_ax].imshow(convert*data, origin='lower', **kwargs)
data[data * convert < vmin * 10.0] = vmin * 10.0 / convert data[data * convert < vmin * 10.0] = vmin * 10.0 / convert
im = ax_curr.imshow(convert * data, origin="lower", **kwargs) 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: if rectangle is not None:
x, y, width, height, angle, color = rectangle[i] x, y, width, height, angle, color = rectangle[i]
ax_curr.add_patch(Rectangle((x, y), width, height, angle=angle, edgecolor=color, fill=False)) ax_curr.add_patch(Rectangle((x, y), width, height, angle=angle, edgecolor=color, fill=False))
# position of centroid # 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([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[1] / 2, data.shape[1] / 2], "--", lw=1, color="grey", alpha=0.5) 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( 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(filt, color="white", fontsize=15, 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(
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.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}$]") 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") fig = plt.figure(figsize=(7 * ratiox, 7 * ratioy), layout="constrained")
if ax is None: if ax is None:
ax = fig.add_subplot(111, projection=wcs) 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) # fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02)
# ax.coords.grid(True, color='white', ls='dotted', alpha=0.5) # ax.coords.grid(True, color='white', ls='dotted', alpha=0.5)
@@ -677,7 +690,7 @@ class align_maps(object):
except KeyError: except KeyError:
for key_i, val_i in value: for key_i, val_i in value:
kwargs[key_i] = val_i 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 [ if kwargs["cmap"] in [
"inferno", "inferno",
@@ -766,7 +779,7 @@ class align_maps(object):
except KeyError: except KeyError:
for key_i, val_i in value: for key_i, val_i in value:
other_kwargs[key_i] = val_i 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_size2 = self.other_wcs.wcs.get_cdelt()[0] * 3600.0
px_sc2 = AnchoredSizeBar( px_sc2 = AnchoredSizeBar(
@@ -1561,6 +1574,7 @@ class overplot_pol(align_maps):
while not self.aligned: while not self.aligned:
self.align() self.align()
self.overplot(levels=levels, P_cut=P_cut, SNRi_cut=SNRi_cut, scale_vec=scale_vec, savename=savename, **kwargs) 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): def add_vector(self, position="center", pol_deg=1.0, pol_ang=0.0, **kwargs):
if isinstance(position, str) and position == "center": 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 = np.zeros(self.img.shape[:2], dtype=bool)
self.mask[self.indices] = True self.mask[self.indices] = True
if hasattr(self, "cont"): if hasattr(self, "cont"):
for coll in self.cont: for coll in self.cont.collections:
coll.remove() coll.remove()
self.cont = self.ax.contour(self.mask.astype(float), levels=[0.5], colors="white", linewidths=1) self.cont = self.ax.contour(self.mask.astype(float), levels=[0.5], colors="white", linewidths=1)
if not self.embedded: if not self.embedded:
@@ -2261,7 +2275,7 @@ class slit(object):
for p in self.pix: 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() 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"): if hasattr(self, "cont"):
for coll in self.cont: for coll in self.cont.collections:
try: try:
coll.remove() coll.remove()
except AttributeError: except AttributeError:
@@ -2364,7 +2378,7 @@ class aperture(object):
x0, y0 = self.circ.center x0, y0 = self.circ.center
self.mask = np.sqrt((xx - x0) ** 2 + (yy - y0) ** 2) < self.radius self.mask = np.sqrt((xx - x0) ** 2 + (yy - y0) ** 2) < self.radius
if hasattr(self, "cont"): if hasattr(self, "cont"):
for coll in self.cont: for coll in self.cont.collections:
try: try:
coll.remove() coll.remove()
except AttributeError: except AttributeError:
@@ -2499,7 +2513,7 @@ class pol_map(object):
self.selected = False self.selected = False
self.region = deepcopy(self.select_instance.mask.astype(bool)) self.region = deepcopy(self.select_instance.mask.astype(bool))
self.select_instance.displayed.remove() self.select_instance.displayed.remove()
for coll in self.select_instance.cont: for coll in self.select_instance.cont.collections:
coll.remove() coll.remove()
self.select_instance.lasso.set_active(False) self.select_instance.lasso.set_active(False)
self.set_data_mask(deepcopy(self.region)) self.set_data_mask(deepcopy(self.region))
@@ -2543,7 +2557,7 @@ class pol_map(object):
self.select_instance.update_mask() self.select_instance.update_mask()
self.region = deepcopy(self.select_instance.mask.astype(bool)) self.region = deepcopy(self.select_instance.mask.astype(bool))
self.select_instance.displayed.remove() self.select_instance.displayed.remove()
for coll in self.select_instance.cont: for coll in self.select_instance.cont.collections:
coll.remove() coll.remove()
self.select_instance.circ.set_visible(False) self.select_instance.circ.set_visible(False)
self.set_data_mask(deepcopy(self.region)) self.set_data_mask(deepcopy(self.region))
@@ -2601,7 +2615,7 @@ class pol_map(object):
self.select_instance.update_mask() self.select_instance.update_mask()
self.region = deepcopy(self.select_instance.mask.astype(bool)) self.region = deepcopy(self.select_instance.mask.astype(bool))
self.select_instance.displayed.remove() self.select_instance.displayed.remove()
for coll in self.select_instance.cont: for coll in self.select_instance.cont.collections:
coll.remove() coll.remove()
self.select_instance.rect.set_visible(False) self.select_instance.rect.set_visible(False)
self.set_data_mask(deepcopy(self.region)) self.set_data_mask(deepcopy(self.region))
@@ -3324,7 +3338,7 @@ class pol_map(object):
) )
if hasattr(self, "cont"): if hasattr(self, "cont"):
for coll in self.cont: for coll in self.cont.collections:
try: try:
coll.remove() coll.remove()
except AttributeError: except AttributeError:

6
package/src/Combine.py
View File

@@ -173,7 +173,7 @@ def main(infiles, target=None, output_dir="./data/"):
# Reduction parameters # Reduction parameters
kwargs = {} kwargs = {}
# Polarization map output # Polarization map output
kwargs["SNRp_cut"] = 3.0 kwargs["P_cut"] = 0.99
kwargs["SNRi_cut"] = 1.0 kwargs["SNRi_cut"] = 1.0
kwargs["flux_lim"] = 1e-19, 3e-17 kwargs["flux_lim"] = 1e-19, 3e-17
kwargs["scale_vec"] = 5 kwargs["scale_vec"] = 5
@@ -186,9 +186,7 @@ def main(infiles, target=None, output_dir="./data/"):
new_infiles = [] new_infiles = []
for i, group in enumerate(grouped_infiles): for i, group in enumerate(grouped_infiles):
new_infiles.append( new_infiles.append(FOC_reduction(target=target + "-" + str(i + 1), infiles=["/".join([data_folder, file]) for file in group], interactive=True)[0])
FOC_reduction(target=target + "-" + str(i + 1), infiles=["/".join([data_folder, file]) for file in group], interactive=True)[0]
)
infiles = new_infiles infiles = new_infiles

48
package/src/overplot_IC5063.py
View File

@@ -11,46 +11,46 @@ from lib.plots import overplot_pol, overplot_radio
from matplotlib.colors import LogNorm from matplotlib.colors import LogNorm
Stokes_UV = fits.open("./data/IC5063/5918/IC5063_FOC_b0.10arcsec_c0.20arcsec.fits") 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_18GHz = fits.open("./data/IC5063/radio/IC5063_18GHz.fits")
Stokes_24GHz = fits.open("./data/IC5063/radio/IC5063_24GHz.fits") # Stokes_24GHz = fits.open("./data/IC5063/radio/IC5063_24GHz.fits")
Stokes_103GHz = fits.open("./data/IC5063/radio/IC5063_103GHz.fits") # Stokes_103GHz = fits.open("./data/IC5063/radio/IC5063_103GHz.fits")
Stokes_229GHz = fits.open("./data/IC5063/radio/IC5063_229GHz.fits") # Stokes_229GHz = fits.open("./data/IC5063/radio/IC5063_229GHz.fits")
Stokes_357GHz = fits.open("./data/IC5063/radio/IC5063_357GHz.fits") # Stokes_357GHz = fits.open("./data/IC5063/radio/IC5063_357GHz.fits")
# Stokes_S2 = fits.open("./data/IC5063/POLARIZATION_COMPARISON/S2_rot_crop.fits") # Stokes_S2 = fits.open("./data/IC5063/POLARIZATION_COMPARISON/S2_rot_crop.fits")
Stokes_IR = fits.open("./data/IC5063/IR/u2e65g01t_c0f_rot.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.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() # levels18GHz = levelsMorganti * Stokes_18GHz[0].data.max()
A = overplot_radio(Stokes_UV, Stokes_18GHz) # 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) # 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() # levels24GHz = levelsMorganti * Stokes_24GHz[0].data.max()
B = overplot_radio(Stokes_UV, Stokes_24GHz) # 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) # 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() # levels103GHz = levelsMorganti * Stokes_103GHz[0].data.max()
C = overplot_radio(Stokes_UV, Stokes_103GHz) # 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) # 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() # levels229GHz = levelsMorganti * Stokes_229GHz[0].data.max()
D = overplot_radio(Stokes_UV, Stokes_229GHz) # 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) # 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() # levels357GHz = levelsMorganti * Stokes_357GHz[0].data.max()
E = overplot_radio(Stokes_UV, Stokes_357GHz) # 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) # 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 = 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)) # 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 = overplot_pol(Stokes_UV, Stokes_IR, cmap="inferno")
G.plot( G.plot(
P_cut=2.0, P_cut=0.99,
SNRi_cut=10.0, SNRi_cut=1.0,
savename="./plots/IC5063/IR_overplot.pdf", savename="./plots/IC5063/IR_overplot.pdf",
scale_vec=None, 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"]), 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",
) )

35
package/src/overplot_MRK463E.py
View File

@@ -11,16 +11,33 @@ from lib.plots import overplot_chandra, overplot_pol
from matplotlib.colors import LogNorm from matplotlib.colors import LogNorm
Stokes_UV = fits.open("./data/MRK463E/5960/MRK463E_FOC_b0.05arcsec_c0.10arcsec.fits") 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_IR = fits.open("./data/MRK463E/WFPC2/IR_rot_crop.fits")
Stokes_Xr = fits.open("./data/MRK463E/Chandra/X_ray_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 = 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.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.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"] 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()) C = overplot_pol(Stokes_UV, Radio, 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") C.other_im.set(norm=LogNorm(1e-4, 2e-2))
B.write_to(path1="./data/MRK463E/FOC_data_WFPC.fits", path2="./data/MRK463E/WFPC_data.fits", suffix="aligned") 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")

6
package/src/overplot_NGC1068.py
View File

@@ -30,10 +30,10 @@ A.add_vector(
A.other_wcs.celestial.wcs.crpix - (1.0, 1.0), A.other_wcs.celestial.wcs.crpix - (1.0, 1.0),
pol_deg=0.124, pol_deg=0.124,
pol_ang=100.7, pol_ang=100.7,
width=2, width=2.0,
linewidth=2, linewidth=1.0,
scale=1.0 / (A.px_scale * 6.0), scale=1.0 / (A.px_scale * 6.0),
edgecolor="b", edgecolor="w",
color="b", color="b",
label=r"IXPE torus: P = $12.4 (\pm 3.6)$%, PA = $100.7 (\pm 8.3)$°", label=r"IXPE torus: P = $12.4 (\pm 3.6)$%, PA = $100.7 (\pm 8.3)$°",
) )