4. Radio Source Maps#
4.1. Notebook setup#
import os
import sys
from cmcrameri import cm as cmc
import emcee
from IPython import display
import matplotlib as mpl
from matplotlib import animation
import matplotlib.pyplot as plt
import numpy as np
import smplotlib
epochs = ["A", "B", "C", "D", "E", "F"]
data_path = "../data/"
fig_path = "../figures/"
calibrator = "J0340"
target = "HR1099"
cm = [
"#377eb8",
"#e41a1c",
"#4daf4a",
"#dede00",
"#ff7f00",
"#999999",
"#984ea3",
"#f781bf",
"#a65628",
]
marker_cycle = ["o", "v", "X", "<", "D", "^"]
sys.path.append(os.path.join(os.getcwd(), ".."))
from library import HR1099_astrometry, utils
deg = np.pi / 180
mas = deg / (60 * 60 * 1000)
4.2. Import data#
epochs, mean_jd = np.genfromtxt(
data_path + target + "_I_positions.txt",
skip_header=2,
dtype="U1,f8,f8,f8,f8,f8,f8,f8",
usecols=(0, 5),
unpack=True,
)
mean_mjd = np.round(mean_jd - 2400000.5, 1)
models = {}
mean_x = np.array([])
mean_y = np.array([])
for epoch in epochs:
beam = utils.get_beam(os.path.join(data_path, f"{epoch}_{target}_sc.log"))
models[epoch] = utils.vlb_model(
os.path.join(data_path, f"{epoch}_{target}_sc.mod"), beam
)
models[epoch].excise_components(200)
x, y, jd, x_err, y_err = np.genfromtxt(
data_path + epoch + "_" + target + "_I_sub-epoch_positions.txt",
skip_header=2,
dtype="U5,f8,f8,f8,f8,f8,f8,f8",
usecols=(2, 3, 5, 6, 7),
unpack=True,
)
mean_x = np.append(mean_x, np.mean(x))
mean_y = np.append(mean_y, np.mean(y))
sampler = emcee.backends.HDFBackend(data_path + target + "_orbital_chain.h5")
flat_samples = sampler.get_chain(flat=True)
flat_samples = np.divide(flat_samples, np.array([1, 1, deg, deg]))
med_val = np.percentile(flat_samples, [16, 50, 84], axis=0)
lower_val = med_val[1] - med_val[0]
upper_val = med_val[2] - med_val[1]
med_val = med_val[1]
hr1099 = HR1099_astrometry.HR1099_astrometry(med_val[2] * deg, med_val[3] * deg)
T0, P, d, a, m1, m2, R1, R2 = hr1099.hr1099_info()
4.3. Plot#
epoch_orbit_phases = []
for i, epoch in enumerate(epochs):
epoch_orbit_phases.append(hr1099.orbit_phase(mean_jd[i]))
sorted_index = np.argsort(epoch_orbit_phases)
fig, ax = plt.subplots(3, 2, figsize=(7, 9.5), sharex=True, sharey=True)
fig.subplots_adjust(wspace=0, hspace=0.1)
xx, yy = np.meshgrid(
np.linspace(-2 * mas, 2 * mas, 256, endpoint=True),
np.linspace(-2 * mas, 2 * mas, 256, endpoint=True),
)
model_sums = [np.sum(m.get_image_plane_model(xx, yy)) for m in models.values()]
max_lev_scalar = np.max(model_sums) / model_sums
for i, epoch in enumerate(epochs):
utils.plot_binary(
ax[i // 2, i % 2],
hr1099,
mean_jd[i],
corotate=False,
centroid=[mean_x[i] + med_val[0], mean_y[i] + med_val[1], 0, 0],
label=mean_mjd[i],
cmap="viridis",
model=models[epoch],
mapsize=4,
cells=256,
levs=[2, 4, 8, 16, 32, 64],
max_lev_scalar=max_lev_scalar[i],
beam=models[epoch].restoring_beam,
show_coord_axes=True,
hour_range=2,
fontsize=8,
d=d,
bar_pos="lower left",
)
ax[i // 2, i % 2].annotate(
"%s (%.1f)" % (epoch, mean_mjd[i]),
xy=(0.6, 0.93),
xycoords="axes fraction",
fontsize=10,
)
ax[i // 2, i % 2].annotate(
r"$\phi$=%.2f" % epoch_orbit_phases[i],
xy=(0.72, 0.86),
xycoords="axes fraction",
fontsize=10,
)
ax[i // 2, i % 2].set_aspect("equal")
ax[i // 2, i % 2].grid(True, linestyle="dotted", alpha=0.5)
fig.set_facecolor("white")
fig.set_dpi(300)
plt.savefig(fig_path + "radio_maps.pdf", bbox_inches="tight")
plt.show()
4.4. Generate movie#
fps = 30 # frames per second
phaseps = 0.03 # phase per second
phase_width = 0.03 # phase width of each epoch
P = P / (24 * 3600) # period in days
epoch_orbit_phases = []
for i, epoch in enumerate(epochs):
epoch_orbit_phases.append(hr1099.orbit_phase(mean_jd[i]))
sorted_index = np.argsort(epoch_orbit_phases)
phase_sorted_indices = np.argsort(epoch_orbit_phases)
show_map = []
jd_range = []
for i in range(len(phase_sorted_indices)):
start = epoch_orbit_phases[phase_sorted_indices[i]] - phase_width
end = epoch_orbit_phases[phase_sorted_indices[i]] + phase_width
if i == 0:
rng = np.linspace(T0, T0 + P * start, int(fps * (start - 0) / phaseps))
s_rng = np.zeros_like(rng)
for j in range(len(rng)):
jd_range.append(rng[j])
show_map.append(s_rng[j])
else:
prev_end = epoch_orbit_phases[phase_sorted_indices[i - 1]] + phase_width
if not prev_end > start:
rng = np.linspace(
T0 + P * prev_end,
T0 + P * start,
int(fps * (start - prev_end) / phaseps),
)
s_rng = np.zeros_like(rng)
for j in range(len(rng)):
jd_range.append(rng[j])
show_map.append(s_rng[j])
else:
start = prev_end
rng = np.linspace(
T0 + P * start, T0 + P * end, int(2.5 * fps * (end - start) / phaseps)
)
s_rng = np.ones_like(rng)
for j in range(len(rng)):
jd_range.append(rng[j])
show_map.append(s_rng[j])
if i == len(phase_sorted_indices) - 1:
rng = np.linspace(T0 + P * end, T0 + P, int(fps * (1 - end) / phaseps))
s_rng = np.zeros_like(rng)
for j in range(len(rng)):
jd_range.append(rng[j])
show_map.append(s_rng[j])
def show_frame(frame):
ax.clear()
phase = hr1099.orbit_phase(jd_range[frame])
if show_map[frame]:
idx = np.argmin(np.abs(epoch_orbit_phases - phase))
utils.plot_binary(
ax,
hr1099,
jd_range[frame],
centroid=[mean_x[idx] + med_val[0], mean_y[idx] + med_val[1], 0, 0],
label=mean_mjd[idx],
model=models[epochs[idx]],
mapsize=4,
cells=256,
levs=[4, 8, 16, 32, 64],
beam=models[epochs[idx]].restoring_beam,
cmap="viridis",
d=d,
fontsize=14,
bar_pos="lower left",
)
ax.annotate(
"MJD=%.1f" % mean_mjd[idx],
xy=(0.05, 0.9),
xycoords="axes fraction",
fontsize=12,
)
else:
utils.plot_binary(ax, hr1099, jd_range[frame], mapsize=4, d=d, fontsize=14)
ax.annotate(
r"$\phi = %.3f$" % phase, xy=(0.05, 0.95), xycoords="axes fraction", fontsize=14
)
ax.grid(True, linestyle="dotted", alpha=0.5)
fig, ax = plt.subplots(figsize=(7, 7))
fig.set_facecolor("white")
anim = animation.FuncAnimation(
fig, show_frame, frames=len(jd_range), interval=1e3 / fps
)
html = anim.to_html5_video()
html = display.HTML(html)
display.display(html)
writervideo = animation.FFMpegWriter(fps=fps)
anim.save(fig_path + "radio_maps.mp4", writer=writervideo)
print("Saved animation to " + fig_path + "radio_maps.mp4")
plt.close()
MovieWriter stderr:
[vost#0:0/libx264 @ 0x12e004920] Error submitting a packet to the muxer: Immediate exit requested
[out#0/mp4 @ 0x6000020dc0c0] Error muxing a packet
Received > 3 system signals, hard exiting
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
File ~/miniforge3/envs/hr1099-timelapse-vlbi/lib/python3.10/site-packages/matplotlib/animation.py:243, in AbstractMovieWriter.saving(self, fig, outfile, dpi, *args, **kwargs)
242 try:
--> 243 yield self
244 finally:
File ~/miniforge3/envs/hr1099-timelapse-vlbi/lib/python3.10/site-packages/matplotlib/animation.py:1105, in Animation.save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs, progress_callback)
1103 for anim, d in zip(all_anim, data):
1104 # TODO: See if turning off blit is really necessary
-> 1105 anim._draw_next_frame(d, blit=False)
1106 if progress_callback is not None:
File ~/miniforge3/envs/hr1099-timelapse-vlbi/lib/python3.10/site-packages/matplotlib/animation.py:1140, in Animation._draw_next_frame(self, framedata, blit)
1139 self._pre_draw(framedata, blit)
-> 1140 self._draw_frame(framedata)
1141 self._post_draw(framedata, blit)
File ~/miniforge3/envs/hr1099-timelapse-vlbi/lib/python3.10/site-packages/matplotlib/animation.py:1768, in FuncAnimation._draw_frame(self, framedata)
1766 # Call the func with framedata and args. If blitting is desired,
1767 # func needs to return a sequence of any artists that were modified.
-> 1768 self._drawn_artists = self._func(framedata, *self._args)
1770 if self._blit:
Cell In[5], line 62, in show_frame(frame)
61 idx = np.argmin(np.abs(epoch_orbit_phases - phase))
---> 62 utils.plot_binary(
63 ax,
64 hr1099,
65 jd_range[frame],
66 centroid=[mean_x[idx] + med_val[0], mean_y[idx] + med_val[1], 0, 0],
67 label=mean_mjd[idx],
68 model=models[epochs[idx]],
69 mapsize=4,
70 cells=256,
71 levs=[4, 8, 16, 32, 64],
72 beam=models[epochs[idx]].restoring_beam,
73 cmap="viridis",
74 d=d,
75 fontsize=14,
76 bar_pos="lower left",
77 )
78 ax.annotate(
79 "MJD=%.1f" % mean_mjd[idx],
80 xy=(0.05, 0.9),
81 xycoords="axes fraction",
82 fontsize=12,
83 )
File ~/Library/CloudStorage/Dropbox/research/HR1099/HR1099-timelapse-vlbi/notebooks/../library/utils.py:147, in plot_binary(ax, binary, jd, corotate, r1, r2, plot_stars, centroid, color, label, model, mapsize, cells, levs, max_lev_scalar, cmap, bg, write_levs, beam, show_coord_axes, hour_range, fontsize, marker, d, bar_pos)
143 xx, yy = np.meshgrid(
144 np.linspace(-mapsize * mas, mapsize * mas, cells, endpoint=True),
145 np.linspace(-mapsize * mas, mapsize * mas, cells, endpoint=True),
146 )
--> 147 im_model = model.get_image_plane_model(
148 xx, yy, x_shift=centroid[0] * mas, y_shift=centroid[1] * mas
149 )
150 if corotate:
KeyboardInterrupt:
During handling of the above exception, another exception occurred:
CalledProcessError Traceback (most recent call last)
Cell In[5], line 103
101 display.display(html)
102 writervideo = animation.FFMpegWriter(fps=fps)
--> 103 anim.save(fig_path + "radio_maps.mp4", writer=writervideo)
104 print("Saved animation to " + fig_path + "radio_maps.mp4")
105 plt.close()
File ~/miniforge3/envs/hr1099-timelapse-vlbi/lib/python3.10/site-packages/matplotlib/animation.py:1089, in Animation.save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs, progress_callback)
1085 savefig_kwargs['transparent'] = False # just to be safe!
1086 # canvas._is_saving = True makes the draw_event animation-starting
1087 # callback a no-op; canvas.manager = None prevents resizing the GUI
1088 # widget (both are likewise done in savefig()).
-> 1089 with writer.saving(self._fig, filename, dpi), \
1090 cbook._setattr_cm(self._fig.canvas, _is_saving=True, manager=None):
1091 for anim in all_anim:
1092 anim._init_draw() # Clear the initial frame
File ~/miniforge3/envs/hr1099-timelapse-vlbi/lib/python3.10/contextlib.py:153, in _GeneratorContextManager.__exit__(self, typ, value, traceback)
151 value = typ()
152 try:
--> 153 self.gen.throw(typ, value, traceback)
154 except StopIteration as exc:
155 # Suppress StopIteration *unless* it's the same exception that
156 # was passed to throw(). This prevents a StopIteration
157 # raised inside the "with" statement from being suppressed.
158 return exc is not value
File ~/miniforge3/envs/hr1099-timelapse-vlbi/lib/python3.10/site-packages/matplotlib/animation.py:245, in AbstractMovieWriter.saving(self, fig, outfile, dpi, *args, **kwargs)
243 yield self
244 finally:
--> 245 self.finish()
File ~/miniforge3/envs/hr1099-timelapse-vlbi/lib/python3.10/site-packages/matplotlib/animation.py:360, in MovieWriter.finish(self)
356 _log.log(
357 logging.WARNING if self._proc.returncode else logging.DEBUG,
358 "MovieWriter stderr:\n%s", err)
359 if self._proc.returncode:
--> 360 raise subprocess.CalledProcessError(
361 self._proc.returncode, self._proc.args, out, err)
CalledProcessError: Command '['ffmpeg', '-f', 'rawvideo', '-vcodec', 'rawvideo', '-s', '1400x1400', '-pix_fmt', 'rgba', '-framerate', '30', '-loglevel', 'error', '-i', 'pipe:', '-vcodec', 'h264', '-pix_fmt', 'yuv420p', '-y', '../figures/radio_maps.mp4']' returned non-zero exit status 123.
