diff --git a/docs/background/velocities.rst b/docs/background/velocities.rst index d4426b3..b1fe938 100644 --- a/docs/background/velocities.rst +++ b/docs/background/velocities.rst @@ -190,7 +190,7 @@ The pulsar's systemic velocity in the plane of the sky The component of this velocity that is parallel to the line of images formed by the lens is then given by - + .. math:: v_\mathrm{p,sys,\parallel} = d_\mathrm{p} @@ -331,7 +331,7 @@ rotation on the sky, with :math:`i_\mathrm{p} = 90^\circ` for an edge-on orbit va='center', vh='left', arrow='-|>', color=col_angmom) rot_lin(ax, th=i_p + 90.*u.deg, s=0.6, name='pulsar\norbit', vh='center', arrow="-", color=col_orbit) - label_angle(ax, th0=180.*u.deg, th1=180.*u.deg + i_p, + label_angle(ax, th0=180.*u.deg, th1=180.*u.deg + i_p, th_name=r"$i_\mathrm{p}$", rad=.2, color=col_angmom) ax.set_aspect('equal') @@ -411,7 +411,7 @@ celestial north through east. th_name=r"$\Omega_\mathrm{p}$", color=col_nodes) rot_lin(ax, th=omega_p+90.*u.deg, name='line of nodes', vh='center', arrow="-|>", color=col_nodes) - + plt.show() In the equatorial coordinate system, the pulsar's orbital sky-plane velocity is @@ -440,7 +440,7 @@ Projecting the sky-plane velocity onto the line of lensed images The component of the pulsar's orbital sky-plane velocity :math:`\vec{v}_\mathrm{p,orb,sky}` that is parallel to the line of images formed by the lens is then given by - + .. math:: \begin{align} @@ -479,7 +479,7 @@ the line of lensed images measured from the ascending node of the pulsar orbit. arrow="-", color=col_screen) rot_lin(ax, th=omega_p+90.*u.deg, name='line of nodes', vh='center', arrow="-|>", color=col_nodes) - + plt.show() @@ -566,7 +566,7 @@ with \qquad \phi_\oplus = 2 \pi \frac{ t - t_\mathrm{asc,\oplus} } { P_\mathrm{orb,\oplus} }, - + .. math:: \chi_\oplus = \arctantwo \left[ @@ -589,7 +589,7 @@ can be derived from the pulsar system's ecliptic coordinates .. plot:: :context: close-figs - + from astropy.coordinates import SkyCoord psr_coord = SkyCoord('04h37m15.99744s -47d15m09.7170s') @@ -639,7 +639,7 @@ can be derived from the pulsar system's ecliptic coordinates plt.plot(1.-0.6, 1., '+', color=col_earth) ax.text(0.2, 1., - "Earth\nat $t_\mathrm{eqx}$", + r"Earth\nat $t_\mathrm{eqx}$", horizontalalignment='center', verticalalignment='center', color=col_earth) @@ -662,11 +662,11 @@ can be derived from the pulsar system's ecliptic coordinates label_angle(ax, th1=90.*u.deg, th0=-beta_p, rad=.5, th_name=r" $i_{\!\!\oplus}$", vh='center', color=col_angmom) ax.set_aspect('equal') - + plt.show() .. container:: align-center - + **Left:** top-down view, looking in the direction of :math:`-\vec{h}_\oplus`. **Right:** side view, looking in the direction of :math:`\hat{x}`. @@ -679,7 +679,7 @@ the Earth's orbital specific angular momentum vector :math:`\vec{h}_\oplus`. It is given by .. math:: - + i_\oplus = \beta_\mathrm{p} + 90^\circ. The restriction on the pulsar's ecliptic latitude :math:`-90^\circ \le @@ -716,7 +716,7 @@ Finally, under the simplifying assumption that Earth's orbit is circular, the time of Earth's passage through the ascending node is given by .. math:: - + t_\mathrm{asc,\oplus} = t_\mathrm{eqx} + P_\mathrm{orb,\oplus} \frac{ \lambda_\mathrm{p} + 90^\circ }{ 360^\circ }, @@ -769,7 +769,7 @@ Filling in the terms for effective proper motion gives + \underbrace{ \vphantom{ \frac{ v_\mathrm{0,p} }{ d_\mathrm{p} } } \frac{ v_{0,\oplus} }{ d_\mathrm{eff} } b_\oplus \sin( \phi_\oplus - \chi_\oplus ) - }_\textrm{Earth's orbital motion}. + }_\textrm{Earth's orbital motion}. This shows that the scaled effective velocity can be written as the normed sum of two sinusoids and a constant offset: diff --git a/docs/tutorials/different_transforms.rst b/docs/tutorials/different_transforms.rst index fcb6395..4543cd6 100644 --- a/docs/tutorials/different_transforms.rst +++ b/docs/tutorials/different_transforms.rst @@ -35,13 +35,7 @@ Start with some standard imports and a handy function for presenting images. from screens.fields import dynamic_field from screens.dynspec import DynamicSpectrum as DS from screens.conjspec import ConjugateSpectrum as CS - - def axis_extent(*args): - result = [] - for a in args: - x = a.squeeze().value - result.extend([x[0] - (dx:=x[1]-x[0])/2, x[-1]+dx/2]) - return result + from screens.visualization import axis_extent Set up a scattering screen diff --git a/docs/tutorials/error_analysis.rst b/docs/tutorials/error_analysis.rst index d15db56..3655b2f 100644 --- a/docs/tutorials/error_analysis.rst +++ b/docs/tutorials/error_analysis.rst @@ -187,15 +187,15 @@ interval for each of the parameters. q_m, q_p = q_50 - q_16, q_84 - q_50 fmt = get_format(fmts, i) if fmt(q_m) == fmt(q_p): - txt = '{0} &= {1} \pm {2} \; {4} \\\\[0.5em]' + txt = r'{0} &= {1} \pm {2} \; {4} \\[0.5em]' else: - txt = '{0} &= {1}_{{-{2}}}^{{+{3}}} \; {4} \\\\[0.5em]' + txt = r'{0} &= {1}_{{-{2}}}^{{+{3}}} \; {4} \\[0.5em]' txt = txt.format(par_strs[i].symbol, fmt(q_50), fmt(q_m), fmt(q_p), par_strs[i].unit) txt_all += txt - txt_all = '\\begin{align}' + txt_all + '\\end{align}' + txt_all = r'\begin{align}' + txt_all + r'\end{align}' display(Math(txt_all)) @@ -203,11 +203,11 @@ interval for each of the parameters. txt_all = '' for i, upar in enumerate(upars): fmt = get_format(fmts, i) - txt_all += (f'{par_strs[i].symbol} &= ' - f'{fmt(upar.n)} \pm {fmt(upar.s)} \; ' - f'{par_strs[i].unit} \\\\[0.5em]') + txt_all += (rf'{par_strs[i].symbol} &= ' + rf'{fmt(upar.n)} \pm {fmt(upar.s)} \; ' + rf'{par_strs[i].unit} \\[0.5em]') - txt_all = '\\begin{align}' + txt_all + '\\end{align}' + txt_all = r'\begin{align}' + txt_all + r'\end{align}' display(Math(txt_all)) @@ -307,7 +307,7 @@ unit strings into label strings for a list of parameters. .. jupyter-execute:: def gen_label_strs(par_strs): - label_strs = [f'${par_str.symbol} \; ({par_str.unit})$' + label_strs = [rf'${par_str.symbol} \; ({par_str.unit})$' for par_str in par_strs] return label_strs diff --git a/docs/tutorials/screen1d.rst b/docs/tutorials/screen1d.rst index 844efc6..b6960c9 100644 --- a/docs/tutorials/screen1d.rst +++ b/docs/tutorials/screen1d.rst @@ -37,15 +37,8 @@ Imports. from screens.screen import Source, Screen1D, Telescope from screens.fields import phasor - -Define a handy function to create extents to use with imshow. - -.. jupyter-execute:: - - def axis_extent(x): - x = x.ravel().value - dx = x[1]-x[0] - return x[0]-0.5*dx, x[-1]+0.5*dx + # A handy function to create extents to use with imshow. + from screens.visualization import axis_extent Construct the system's components @@ -246,7 +239,7 @@ Plot the dynamic spectrum. plt.imshow(dynspec.T, origin='lower', aspect='auto', interpolation='none', - cmap='Greys', extent=axis_extent(t) + axis_extent(f), vmin=0.) + cmap='Greys', extent=axis_extent(t, f), vmin=0.) plt.xlabel(rf"time $t$ ({t.unit.to_string('latex')})") plt.ylabel(rf"frequency $f$ ({f.unit.to_string('latex')})") @@ -279,7 +272,7 @@ Plot the secondary spectrum. plt.imshow(secspec.T, origin='lower', aspect='auto', interpolation='none', - cmap='Greys', extent=axis_extent(fd) + axis_extent(tau), + cmap='Greys', extent=axis_extent(fd, tau), norm=LogNorm(vmin=1.e-4, vmax=1.)) plt.xlim(-5., 5.) plt.ylim(-15., 15.) diff --git a/docs/tutorials/two_screens.rst b/docs/tutorials/two_screens.rst index 41381d1..94130c6 100644 --- a/docs/tutorials/two_screens.rst +++ b/docs/tutorials/two_screens.rst @@ -49,16 +49,8 @@ Imports. from screens.screen import Source, Screen1D, Telescope from screens.fields import phasor - -Define a handy function to help create an `extent` for -:py:func:`matplotlib.pyplot.imshow`. - -.. jupyter-execute:: - - def axis_extent(x): - x = x.ravel().value - dx = x[1]-x[0] - return x[0]-0.5*dx, x[-1]+0.5*dx + # A handy function to help create an extent for matplotlib.pyplot.imshow + from screens.visualization import axis_extent Define a matrix to transform between coordinate systems. This is needed because Astropy's :py:class:`~astropy.coordinates.SkyOffsetFrame` yields frames where @@ -414,7 +406,7 @@ Plot the dynamic spectrum. plt.imshow(dynspec.T, origin='lower', aspect='auto', interpolation='none', - cmap='Greys', extent=axis_extent(t) + axis_extent(f), vmin=0.) + cmap='Greys', extent=axis_extent(t, f), vmin=0.) plt.xlabel(rf"time $t$ ({t.unit.to_string('latex')})") plt.ylabel(rf"frequency $f$ ({f.unit.to_string('latex')})") @@ -447,7 +439,7 @@ Plot the secondary spectrum. plt.imshow(secspec.T.value, origin='lower', aspect='auto', interpolation='none', - cmap='Greys', extent=axis_extent(fd) + axis_extent(tau), + cmap='Greys', extent=axis_extent(fd, tau), norm=LogNorm(vmin=1.e-8, vmax=1.)) plt.xlim(-50., 50.) @@ -470,7 +462,7 @@ Overplot the arclet apices. plt.imshow(secspec.T.value, origin='lower', aspect='auto', interpolation='none', - cmap='Greys', extent=axis_extent(fd) + axis_extent(tau), + cmap='Greys', extent=axis_extent(fd, tau), norm=LogNorm(vmin=1.e-8, vmax=1.)) plt.xlim(-50., 50.) diff --git a/examples/two_screens.py b/examples/two_screens.py index 2d70aca..4bcbbb7 100644 --- a/examples/two_screens.py +++ b/examples/two_screens.py @@ -20,6 +20,7 @@ from screens.screen import Source, Screen1D, Telescope from screens.fields import phasor +from screens.visualization import axis_extent dp = 0.75*u.kpc @@ -41,12 +42,6 @@ [0.85]*u.AU, magnification=0.05)[..., np.newaxis] -def axis_extent(x): - x = x.ravel().value - dx = x[1]-x[0] - return x[0]-0.5*dx, x[-1]+0.5*dx - - def unit_vector(c): return c.represent_as(UnitSphericalRepresentation).to_cartesian() @@ -194,7 +189,7 @@ def plot_screen(ax, s, d, color='black', **kwargs): ds = np.abs(dw.sum(0))**2 ax_ds = plt.subplot(233) ax_ds.imshow(ds.T, cmap='Greys', - extent=axis_extent(t) + axis_extent(f), + extent=axis_extent(t, f), origin='lower', interpolation='none', aspect='auto') ax_ds.set_xlabel(t.unit.to_string('latex')) ax_ds.set_ylabel(f.unit.to_string('latex')) @@ -206,7 +201,7 @@ def plot_screen(ax, s, d, color='black', **kwargs): fd = np.fft.fftshift(np.fft.fftfreq(t.size, t[1]-t[0])).to(u.mHz) ax_ss = plt.subplot(236) ax_ss.imshow(np.log10(np.abs(cs.T)**2), vmin=-7, vmax=0, cmap='Greys', - extent=axis_extent(fd) + axis_extent(tau), + extent=axis_extent(fd, tau), origin='lower', interpolation='none', aspect='auto') ax_ss.set_xlim(-5, 5) ax_ss.set_ylim(-10, 10) diff --git a/screens/visualization/__init__.py b/screens/visualization/__init__.py index 911ef54..a52b10c 100644 --- a/screens/visualization/__init__.py +++ b/screens/visualization/__init__.py @@ -3,3 +3,4 @@ from .thetatheta import ThetaTheta # noqa from .sketch import make_sketch # noqa from .hue_cycle_cmap import phasecmap # noqa +from .utils import axis_extent # noqa diff --git a/screens/visualization/utils.py b/screens/visualization/utils.py new file mode 100644 index 0000000..a8d3fbf --- /dev/null +++ b/screens/visualization/utils.py @@ -0,0 +1,25 @@ +from astropy.units import Quantity + + +def axis_extent(*args): + """Treat the arguments as axis arrays of an image, and get their extent. + + Just the first and last element of each array, but properly taking into + account that the limits of the image are half a pixel before and after. + + Parameters + ---------- + args : array + Arrays with coordinates of the images. Assumed to be contiguous, + and have only one dimension with non-unity shape. + + Returns + ------- + extent : list of int + Suitable for passing into matplotlib's ``imshow``. + """ + result = [] + for a in args: + x = Quantity(a).squeeze().value + result.extend([x[0] - (x[1]-x[0])/2, x[-1] + (x[-1]-x[-2])/2]) + return result