In [1]:
%matplotlib inline
import numpy as np
import matplotlib as mpl
from matplotlib import pyplot as plt
from astropy import units as u
from astropy import constants as c
from astropy.coordinates import SkyCoord
from astropy.analytic_functions import blackbody_lambda, blackbody_nu
In [2]:
mpl.rc('figure', figsize=(12, 8))
mpl.rc('lines', linewidth=4)
mpl.rc('font', size=20)
mpl.rc('axes.formatter', limits=(-4, 4))
In [3]:
def texify(unit):
unit_strings = []
for unit, power in zip(intensity_unit.bases, intensity_unit.powers):
if power == 1:
unit_strings.append(unit.to_string())
else:
unit_strings.append(unit.to_string() + "$^{{{:d}}}$".format(power))
return " ".join(unit_strings)
Astropy quantitites are a great way to handle all sorts of messy unit conversions. Careful unit conversions save lives! https://en.wikipedia.org/wiki/Gimli_Glider
The simplest way to create a new quantity object is multiply or divide a number by a Unit instance.
In [4]:
print(type(u.Msun))
u.Msun
Out[4]:
A brief word of warning, the pretty printing shown above will take longer than just printing out the values of an array. This is unnoticable in this case, but is evident if you try to print a large array of quantities.
In [5]:
mass = 1 * u.Msun
print(type(mass))
mass
Out[5]:
In [6]:
# quantities subclass numpy ndarray, so you can make handle arrays of quantities like you would
# any other array object
mass.__class__.__bases__
Out[6]:
In [7]:
# we can convert units to other equivalent units
mass.to(u.kg)
Out[7]:
In [8]:
# there are shortcuts for converting to the relevent system, regardless of what type of quantity it is
mass.cgs
Out[8]:
In [9]:
mass.si
Out[9]:
In [10]:
# we can inspect their unit and their numeric value with that unit
print(mass.value, mass.unit)
In [11]:
# calculations with quantities can produce quantities with new units
average_density = mass / (4 / 3 * np.pi * u.Rearth ** 3)
average_density.cgs
Out[11]:
In [12]:
# Newton's constant
c.G
Out[12]:
In [13]:
# Planck's constant
c.h
Out[13]:
In [14]:
# speed of light
c.c
Out[14]:
In [15]:
# made into a quantity array by multiplying numeric array by unit
R = np.linspace(1, 5) * u.Rearth
v = np.sqrt(2 * c.G * u.Msun / R)
print(v)
v = v.to(u.km / u.s)
print(v)
plt.plot(R, v)
plt.xlabel(r"Radius [R$_\oplus$]")
plt.ylabel(r"Escape velocity [km s$^{-1}$]")
Out[15]:
In [16]:
obscure_quantity = 42 * c.G * c.m_e ** 2 / c.k_B ** 2 * c.c ** 3 * (5700 * u.K) ** -2 * u.Msun / u.Mpc
In [17]:
obscure_quantity
Out[17]:
In [18]:
# what the heck is a m^6 kg Msun / (Mpc J^2 s^5)??
obscure_quantity.decompose()
Out[18]:
In [19]:
# will fail!
obscure_quantity.to(u.m)
In [20]:
# addition works for like units
(1 * u.m) + (1 * u.cm)
Out[20]:
In [21]:
# and fails for the wrong dimensions
(1 * u.m) + (1 * u.s)
Equivalencies allow you to do unit conversions under certain physical assumptions. For instance, it makes sense to talk about converting wavelength to frequency when you are discussing the properties of light waves in a vacuum. See http://docs.astropy.org/en/stable/units/equivalencies.html#unit-equivalencies.
In [22]:
wavelengths = np.linspace(0.1, 1, 100) * u.micron
# will fail without the correct equivalency passed in!
frequencies = wavelengths.to(u.THz, equivalencies=u.spectral())
plt.plot(wavelengths, frequencies)
plt.xlabel(r"$\lambda$ [$\mu$m]")
plt.ylabel(r"$\nu$ [THz]")
Out[22]:
In [23]:
intensity_unit = blackbody_nu(wavelengths[0], temperature=1e3 * u.K).unit
wavelengths = np.logspace(-1, 1, 100) * u.micron
temperatures = np.linspace(5e3, 1e4, 5) * u.K
for T in temperatures:
plt.plot(wavelengths, blackbody_nu(wavelengths, temperature=T),
label='{:.2e}'.format(T))
plt.xscale('log')
plt.legend(loc='best')
plt.xlabel(r'$\lambda$ [$\mu$m]')
plt.ylabel('$I_\\nu$ [{}]'.format(texify(intensity_unit)))
Out[23]:
In [24]:
intensity_unit = blackbody_lambda(wavelengths[0], temperature=1e3 * u.K).unit
for T in temperatures:
plt.plot(wavelengths, blackbody_lambda(wavelengths, temperature=T),
label='{:.2e}'.format(T))
plt.xscale('log')
plt.legend(loc='best')
plt.xlabel(r'$\lambda$ [$\mu$m]')
plt.ylabel('$I_\\lambda$ [{}]'.format(texify(intensity_unit)))
Out[24]:
In [25]:
T = 1e4 * u.K
solid_angle = ((1 * u.arcsec) ** 2).to(u.sr)
f_nu = blackbody_nu(wavelengths, temperature=T) * solid_angle
f_lambda = blackbody_lambda(wavelengths, temperature=T) * solid_angle
print(f_nu.unit)
print(f_lambda.unit)
# I_nu.to(I_lambda.unit) # would fail
# for conversion of spectral energy density, we need to specify what part of the spectra we're looking at
f_lambda_converted = f_nu.to(f_lambda.unit, equivalencies=u.spectral_density(wavelengths))
print(f_lambda_converted.unit)
# shouldn't raise any exceptions!
assert np.all(np.isclose(f_lambda.value, f_lambda_converted.value))
In [26]:
# e.g., from the example above
np.isclose(f_lambda, f_lambda_converted)
Even when you think you have a dimensionless array, it can still be a dimensionless quantity.
In [27]:
print((1 * u.m) / (2 * u.m), type((1 * u.m) / (2 * u.m)))
In [29]:
@u.quantity_input(angle=u.arcsec, distance=u.Mpc)
def angle_to_size(angle, distance):
return angle.to(u.radian).value * distance
# this should work
angle_to_size(1 * u.arcsec, 25 * u.Mpc).to(u.kpc)
Out[29]:
In [30]:
# quantity_input only checks for convertability, not that it's the same unit
angle_to_size(1 * u.arcmin, 25 * u.Mpc).to(u.kpc)
Out[30]:
In [31]:
# this should raise an error
angle_to_size(1 * u.m, 25 * u.Mpc)
In [32]:
coord = SkyCoord(45, 30, unit=u.deg)
In [33]:
# ICRS is the reference frame
coord
Out[33]:
In [34]:
# we can transform between coordinate frames
coord.fk4
Out[34]:
In [35]:
coord.fk5
Out[35]:
In [36]:
coord.galactic
Out[36]:
In [37]:
# latitude and longitude are accessed with ra and dec (when in icrs or fk frames)
coord.ra
Out[37]:
In [38]:
coord.dec
Out[38]:
The attributes ra and dec are Angles. They are subclasses of Quantity, and so they behave similarly, but have more specific functionality. See http://docs.astropy.org/en/stable/coordinates/angles.html#working-with-angles for more details.
You can get nice string representations of angles for all your inane legacy software requirements.
In [40]:
print(coord.to_string())
print(coord.to_string('dms'))
print(coord.to_string('hmsdms'))
print(coord.to_string('hmsdms', sep=':'))
print(coord.to_string('hmsdms', sep=' '))
There are lots of specific use cases outlined here, but let's go over a simple catalog matching exercise.
In [41]:
# need network connection
center_coord = SkyCoord.from_name('M31')
center_coord
Out[41]:
In [42]:
# some mock coordinates
n = 500
ra_values = np.random.randn(n) + center_coord.ra.deg
dec_values = np.random.randn(n) + center_coord.dec.deg
coords = SkyCoord(ra_values, dec_values, unit=u.deg)
plt.scatter(coords.ra.deg, coords.dec.deg, s=100,
edgecolor='k', label='Parent sample')
plt.xlim(plt.xlim()[::-1]) # ra increases right to left
plt.xlabel("Right ascension [deg]")
plt.ylabel("Declination [deg]")
Out[42]:
In [43]:
# mock measurements
n_sample = 100
astrometric_noise = 1 * u.arcsec
sample_indices = np.random.choice(np.arange(len(coords)), n_sample)
sample_ra = coords[sample_indices].ra.deg
sample_dec = coords[sample_indices].dec.deg
angles = 2 * np.pi * np.random.rand(n_sample)
dr = astrometric_noise.to(u.deg).value * np.random.randn(n_sample)
dx = np.cos(angles) * dr - np.sin(angles) * dr
dy = np.sin(angles) * dr + np.cos(angles) * dr
sample_coords = SkyCoord(sample_ra + dx, sample_dec + dy, unit=u.deg)
plt.scatter(coords.ra.deg, coords.dec.deg, s=100,
edgecolor='k', marker='o', alpha=0.8, label='Parent sample')
plt.scatter(sample_coords.ra.deg, sample_coords.dec.deg, s=100,
edgecolor='k', marker='v', alpha=0.8, label='Child sample')
plt.xlim(plt.xlim()[::-1]) # ra increases right to left
plt.xlabel("Right ascension [deg]")
plt.ylabel("Declination [deg]")
plt.legend(bbox_to_anchor=(1, 1))
Out[43]:
In [47]:
# match_to_catalog_sky will return indices into coords of the closest matching objects,
# the angular separation, and the physical distance (ignored here)
idx, sep, dist = sample_coords.match_to_catalog_sky(coords)
ideal_sep = astrometric_noise.to(u.deg) * np.random.randn(int(1e6))
plt.hist(np.abs(ideal_sep.to(u.arcsec)), histtype='step', lw=4, bins='auto', normed=True, label='Ideal')
plt.hist(sep.arcsec, histtype='step', lw=4, bins='auto', normed=True, label='Data')
plt.xlim(0, 3)
plt.xlabel("Separation [arcsec]")
plt.legend(loc='best')
Out[47]:
In [46]:
plt.scatter(coords[idx].ra.deg, coords[idx].dec.deg, s=100,
edgecolor='k', marker='o', alpha=0.8, label='Matched parent sample')
plt.scatter(sample_coords.ra.deg, sample_coords.dec.deg, s=100,
edgecolor='k', marker='v', alpha=0.8, label='Child sample')
plt.xlim(plt.xlim()[::-1]) # ra increases right to left
plt.xlabel("Right ascension [deg]")
plt.ylabel("Declination [deg]")
plt.legend(bbox_to_anchor=(1, 1))
Out[46]:
In [ ]: