Skip to content

Commit

Permalink
Merge pull request #80 from mhvk/axis-extent
Browse files Browse the repository at this point in the history 
Move axis_extent from examples/docs to actual utility.
  • Loading branch information
@mhvk
mhvk authored Sep 2, 2024
2 parents 208f0ed + 2a59b91 commit fefe764
Show file tree
Hide file tree
Showing 8 changed files with 60 additions and 60 deletions.
26 changes: 13 additions & 13 deletions docs/background/velocities.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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}
Expand Down Expand Up @@ -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()


Expand Down Expand Up @@ -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[
Expand All @@ -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')
Expand Down Expand Up @@ -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)
Expand All @@ -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}`.
Expand All @@ -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
Expand Down Expand Up @@ -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 },
Expand Down Expand Up @@ -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:
Expand Down
8 changes: 1 addition & 7 deletions docs/tutorials/different_transforms.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
16 changes: 8 additions & 8 deletions docs/tutorials/error_analysis.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,27 +187,27 @@ 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))


def display_ufloats(upars, par_strs, fmts):
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))


Expand Down Expand Up @@ -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

Expand Down
15 changes: 4 additions & 11 deletions docs/tutorials/screen1d.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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')})")

Expand Down Expand Up @@ -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.)
Expand Down
18 changes: 5 additions & 13 deletions docs/tutorials/two_screens.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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')})")

Expand Down Expand Up @@ -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.)
Expand All @@ -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.)
Expand Down
11 changes: 3 additions & 8 deletions examples/two_screens.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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()

Expand Down Expand Up @@ -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'))
Expand All @@ -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)
Expand Down
1 change: 1 addition & 0 deletions screens/visualization/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
25 changes: 25 additions & 0 deletions screens/visualization/utils.py
View file
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit fefe764

Please sign in to comment.