Astropy is a bundle of several kinds of functionality for Astronomy. It consists of the core packages, providing base functionality for everyone, and affiliated packages which ontain more specialized functionality. The core packages have a minimum of dependencies, while affiliated packages are allowed to depend on external libraries.
A selection of core packages:
astropy.constants
astropy.units
astropy.nddata
- Numpy arrays on steroidsastropy.io
- file input/output, i.e. FITS and ASCIIastropy.convolution
- smoothing, filtering etc.astropy.cosmology
astropy.modeling
- modelling and fittingExamples of affiliated packages:
...etc. etc.
Obviously, there is way too much functionality in Astropy to do more than just scratch the paint.
We have already looked at astropy.io.fits
. We will now take a little look at:
astropy.units
astropy.cosmology
astropy.units
¶import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from astropy import units as u
from astropy import constants as con
# Only when running in n IPython notebook:
%matplotlib inline
# Just for good looks, no need to worry about this now.
# You are welcome to do this but it is not necessary.
c1 = cm.Reds(0.6)
c2 = cm.Blues(0.6)
c3 = cm.PuOr(0.3)
plt.rc('xtick', labelsize=13)
plt.rc('ytick', labelsize=13)
plt.rc('font', size = 12)
plt.rc('lines', linewidth=3)
plt.rc('figure', figsize=(10.0, 7.0))
plt.rcParams['axes.formatter.limits'] = (-3,3) #
dist = 3.4159 * u.km
dist
time = 2.7183 * u.h
time
dist / time
(_).to(u.km / u.s)
We can of course define custom units for our convenience:
kmph = u.Unit(u.km / u.hour)
kms = u.Unit(u.km / u.s)
mph = u.Unit(u.imperial.mile / u.hour)
(200 * kmph).to(mph)
It handles NumPy arrays as well as single scalars:
speeds = np.arange(100, 200) * kmph
print(speeds.to(mph))
[ 62.13711922 62.75849042 63.37986161 64.0012328 64.62260399 65.24397518 65.86534638 66.48671757 67.10808876 67.72945995 68.35083115 68.97220234 69.59357353 70.21494472 70.83631592 71.45768711 72.0790583 72.70042949 73.32180068 73.94317188 74.56454307 75.18591426 75.80728545 76.42865665 77.05002784 77.67139903 78.29277022 78.91414141 79.53551261 80.1568838 80.77825499 81.39962618 82.02099738 82.64236857 83.26373976 83.88511095 84.50648214 85.12785334 85.74922453 86.37059572 86.99196691 87.61333811 88.2347093 88.85608049 89.47745168 90.09882287 90.72019407 91.34156526 91.96293645 92.58430764 93.20567884 93.82705003 94.44842122 95.06979241 95.6911636 96.3125348 96.93390599 97.55527718 98.17664837 98.79801957 99.41939076 100.04076195 100.66213314 101.28350433 101.90487553 102.52624672 103.14761791 103.7689891 104.3903603 105.01173149 105.63310268 106.25447387 106.87584506 107.49721626 108.11858745 108.73995864 109.36132983 109.98270103 110.60407222 111.22544341 111.8468146 112.46818579 113.08955699 113.71092818 114.33229937 114.95367056 115.57504176 116.19641295 116.81778414 117.43915533 118.06052653 118.68189772 119.30326891 119.9246401 120.54601129 121.16738249 121.78875368 122.41012487 123.03149606 123.65286726] mi / h
print speeds.to(u.lightyear/u.millisecond)
[ 2.93611343e-18 2.96547456e-18 2.99483570e-18 3.02419683e-18 3.05355796e-18 3.08291910e-18 3.11228023e-18 3.14164137e-18 3.17100250e-18 3.20036364e-18 3.22972477e-18 3.25908590e-18 3.28844704e-18 3.31780817e-18 3.34716931e-18 3.37653044e-18 3.40589158e-18 3.43525271e-18 3.46461384e-18 3.49397498e-18 3.52333611e-18 3.55269725e-18 3.58205838e-18 3.61141952e-18 3.64078065e-18 3.67014178e-18 3.69950292e-18 3.72886405e-18 3.75822519e-18 3.78758632e-18 3.81694746e-18 3.84630859e-18 3.87566972e-18 3.90503086e-18 3.93439199e-18 3.96375313e-18 3.99311426e-18 4.02247540e-18 4.05183653e-18 4.08119766e-18 4.11055880e-18 4.13991993e-18 4.16928107e-18 4.19864220e-18 4.22800334e-18 4.25736447e-18 4.28672560e-18 4.31608674e-18 4.34544787e-18 4.37480901e-18 4.40417014e-18 4.43353128e-18 4.46289241e-18 4.49225354e-18 4.52161468e-18 4.55097581e-18 4.58033695e-18 4.60969808e-18 4.63905922e-18 4.66842035e-18 4.69778148e-18 4.72714262e-18 4.75650375e-18 4.78586489e-18 4.81522602e-18 4.84458716e-18 4.87394829e-18 4.90330942e-18 4.93267056e-18 4.96203169e-18 4.99139283e-18 5.02075396e-18 5.05011510e-18 5.07947623e-18 5.10883736e-18 5.13819850e-18 5.16755963e-18 5.19692077e-18 5.22628190e-18 5.25564304e-18 5.28500417e-18 5.31436530e-18 5.34372644e-18 5.37308757e-18 5.40244871e-18 5.43180984e-18 5.46117098e-18 5.49053211e-18 5.51989324e-18 5.54925438e-18 5.57861551e-18 5.60797665e-18 5.63733778e-18 5.66669892e-18 5.69606005e-18 5.72542118e-18 5.75478232e-18 5.78414345e-18 5.81350459e-18 5.84286572e-18] lyr / ms
type(speeds)
astropy.units.quantity.Quantity
Astropy's Quantity
objects plays nicely with the rest of the SciPy ecosystem by disguising themselves as numpy.ndarray
's.
No need to convert before sending them to e.g. Matplotlib. They also play nicely with most simple numpy functions:
plt.plot(speeds, np.sqrt(speeds.to(mph)))
[<matplotlib.lines.Line2D at 0x7f8e6a278f90>]
...But not all, which makes sense because log(km/h) is meaningless:
np.sqrt(speeds)
<Quantity [ 10. , 10.04987562, 10.09950494, 10.14889157, 10.19803903, 10.24695077, 10.29563014, 10.34408043, 10.39230485, 10.44030651, 10.48808848, 10.53565375, 10.58300524, 10.63014581, 10.67707825, 10.72380529, 10.77032961, 10.81665383, 10.86278049, 10.90871211, 10.95445115, 11. , 11.04536102, 11.09053651, 11.13552873, 11.18033989, 11.22497216, 11.26942767, 11.3137085 , 11.35781669, 11.40175425, 11.44552314, 11.48912529, 11.53256259, 11.5758369 , 11.61895004, 11.66190379, 11.70469991, 11.74734012, 11.78982612, 11.83215957, 11.87434209, 11.91637529, 11.95826074, 12. , 12.04159458, 12.08304597, 12.12435565, 12.16552506, 12.20655562, 12.24744871, 12.28820573, 12.32882801, 12.36931688, 12.40967365, 12.4498996 , 12.489996 , 12.52996409, 12.56980509, 12.60952021, 12.64911064, 12.68857754, 12.72792206, 12.76714533, 12.80624847, 12.84523258, 12.88409873, 12.92284798, 12.9614814 , 13. , 13.03840481, 13.07669683, 13.11487705, 13.15294644, 13.19090596, 13.22875656, 13.26649916, 13.3041347 , 13.34166406, 13.37908816, 13.41640786, 13.45362405, 13.49073756, 13.52774926, 13.56465997, 13.60147051, 13.6381817 , 13.67479433, 13.7113092 , 13.74772708, 13.78404875, 13.82027496, 13.85640646, 13.89244399, 13.92838828, 13.96424004, 14. , 14.03566885, 14.07124728, 14.10673598] km(1/2) / h(1/2)>
np.log10(speeds)
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-9-abd8bbabc22f> in <module>() ----> 1 np.log10(speeds) /home/trive/anaconda/lib/python2.7/site-packages/astropy/units/quantity.pyc in __array_prepare__(self, obj, context) 281 # the unit the output from the ufunc will have. 282 if function in UFUNC_HELPERS: --> 283 scales, result_unit = UFUNC_HELPERS[function](function, *units) 284 else: 285 raise TypeError("Unknown ufunc {0}. Please raise issue on " /home/trive/anaconda/lib/python2.7/site-packages/astropy/units/quantity_helper.pyc in helper_dimensionless_to_dimensionless(f, unit) 87 raise TypeError("Can only apply '{0}' function to " 88 "dimensionless quantities" ---> 89 .format(f.__name__)) 90 return [scale], dimensionless_unscaled 91 TypeError: Can only apply 'log10' function to dimensionless quantities
The solution to this is a little annoying but there's no way around it:
logspeeds = np.log10(speeds.to(mph).value) # No units anymore!
print(logspeeds)
plt.plot(speeds, logspeeds)
plt.show()
[ 1.79335111 1.79767249 1.80195129 1.80618834 1.81038445 1.81454041 1.81865698 1.82273489 1.82677487 1.83077761 1.8347438 1.83867409 1.84256914 1.84642956 1.85025597 1.85404896 1.8578091 1.86153698 1.86523312 1.86889808 1.87253236 1.87613649 1.87971095 1.88325623 1.8867728 1.89026113 1.89372166 1.89715484 1.90056108 1.90394083 1.90729447 1.91062241 1.91392505 1.91720276 1.92045591 1.92368488 1.92689002 1.93007168 1.9332302 1.93636592 1.93947915 1.94257023 1.94563946 1.94868715 1.95171361 1.95471912 1.95770397 1.96066845 1.96361283 1.96653738 1.96944237 1.97232806 1.9751947 1.97804255 1.98087184 1.98368281 1.98647571 1.98925077 1.9920082 1.99474824 1.9974711 2.00017699 2.00286613 2.00553872 2.00819496 2.01083506 2.0134592 2.01606759 2.0186604 2.02123782 2.02380004 2.02634723 2.02887956 2.03139722 2.03390036 2.03638916 2.03886378 2.04132438 2.04377112 2.04620415 2.04862362 2.05102969 2.0534225 2.0558022 2.05816894 2.06052284 2.06286406 2.06519272 2.06750896 2.06981292 2.07210472 2.07438448 2.07665234 2.07890842 2.08115284 2.08338573 2.08560719 2.08781734 2.09001631 2.09220419]
astropy.cosmology
¶from astropy.cosmology import Planck13, LambdaCDM
pl13 = Planck13
redshifts = np.linspace(0.001, 100, 100)
pl13.lookback_time(7.8)
pl13.angular_diameter_distance(redshifts)
<Quantity [ 4.41823263, 1701.89483502, 1768.15666645, 1622.24519705, 1461.51845089, 1319.11084429, 1198.20127448, 1096.10255984, 1009.43135739, 935.23098014, 871.1271297 , 815.25727217, 766.16458236, 722.70378618, 683.96706805, 649.22804468, 617.89990425, 589.50413756, 563.64706083, 540.00205354, 518.29600201, 498.2988584 , 479.81552582, 462.67949537, 446.74781434, 431.89707476, 418.02019063, 405.02378999, 392.82608989, 381.35515352, 370.54745206, 360.34667099, 350.70271389, 341.5708667 , 332.91109316, 324.68743805, 316.86751948, 309.42209513, 302.32469024, 295.55127716, 289.07999861, 282.8909275 , 276.96585816, 271.28812395, 265.84243768, 260.61475154, 255.59213367, 250.76265928, 246.11531429, 241.6399098 , 237.32700606, 233.16784473, 229.15428827, 225.27876584, 221.53422466, 217.91408632, 214.4122075 , 211.02284444, 207.74062096, 204.56049941, 201.47775434, 198.48794871, 195.58691206, 192.77072083, 190.03568023, 187.37830781, 184.79531837, 182.28361017, 179.84025223, 177.46247274, 175.14764838, 172.89329443, 170.6970558 , 168.55669856, 166.47010232, 164.43525302, 162.45023638, 160.51323176, 158.62250648, 156.77641056, 154.97337186, 153.2118915 , 151.49053965, 149.80795155, 148.16282391, 146.55391143, 144.9800236 , 143.44002173, 141.93281615, 140.4573636 , 139.01266476, 137.59776197, 136.21173707, 134.85370939, 133.52283384, 132.21829911, 130.93932603, 129.68516599, 128.45509941, 127.24843438] Mpc>
plt.plot(redshifts, pl13.angular_diameter_distance(redshifts))
[<matplotlib.lines.Line2D at 0x7f8e668c3150>]
pl13.comoving_distance(redshifts).to(u.lightyear)
<Quantity [ 1.44247579e+07, 1.11632409e+10, 1.74230218e+10, 2.13296904e+10, 2.40313519e+10, 2.60355611e+10, 2.75965894e+10, 2.88561678e+10, 2.98999955e+10, 3.07832293e+10, 3.15431515e+10, 3.22059726e+10, 3.27907228e+10, 3.33115945e+10, 3.37794155e+10, 3.42026086e+10, 3.45878370e+10, 3.49404502e+10, 3.52647993e+10, 3.55644639e+10, 3.58424199e+10, 3.61011641e+10, 3.63428092e+10, 3.65691562e+10, 3.67817516e+10, 3.69819313e+10, 3.71708567e+10, 3.73495423e+10, 3.75188792e+10, 3.76796537e+10, 3.78325625e+10, 3.79782253e+10, 3.81171958e+10, 3.82499701e+10, 3.83769943e+10, 3.84986709e+10, 3.86153640e+10, 3.87274038e+10, 3.88350910e+10, 3.89386994e+10, 3.90384795e+10, 3.91346608e+10, 3.92274539e+10, 3.93170526e+10, 3.94036356e+10, 3.94873677e+10, 3.95684017e+10, 3.96468788e+10, 3.97229304e+10, 3.97966785e+10, 3.98682366e+10, 3.99377106e+10, 4.00051995e+10, 4.00707957e+10, 4.01345857e+10, 4.01966505e+10, 4.02570663e+10, 4.03159045e+10, 4.03732323e+10, 4.04291130e+10, 4.04836063e+10, 4.05367683e+10, 4.05886522e+10, 4.06393081e+10, 4.06887835e+10, 4.07371235e+10, 4.07843704e+10, 4.08305648e+10, 4.08757449e+10, 4.09199472e+10, 4.09632062e+10, 4.10055549e+10, 4.10470246e+10, 4.10876451e+10, 4.11274448e+10, 4.11664510e+10, 4.12046894e+10, 4.12421849e+10, 4.12789610e+10, 4.13150403e+10, 4.13504443e+10, 4.13851939e+10, 4.14193088e+10, 4.14528080e+10, 4.14857096e+10, 4.15180312e+10, 4.15497894e+10, 4.15810003e+10, 4.16116794e+10, 4.16418416e+10, 4.16715010e+10, 4.17006714e+10, 4.17293660e+10, 4.17575976e+10, 4.17853783e+10, 4.18127199e+10, 4.18396339e+10, 4.18661311e+10, 4.18922222e+10, 4.19179173e+10] lyr>
plt.plot(redshifts, pl13.comoving_distance(redshifts).to(u.lightyear))
[<matplotlib.lines.Line2D at 0x7f8e667dded0>]
pl13.angular_diameter_distance(8.4)
plt.plot(redshifts, pl13.luminosity_distance(redshifts))
[<matplotlib.lines.Line2D at 0x7f8e6678fb50>]
redshift = 0.3
pl13.comoving_distance(redshift)
Question: What is the current recession velocity at this distance?
def vel_now(z):
return pl13.H0 * pl13.comoving_distance(z)
vel_now(redshift)
But not to worry, the more interesting question is: What was the recession velocity at the time of emission:
def vel_em(z):
return pl13.H(z) * pl13.scale_factor(z) * pl13.comoving_distance(z)
vel_em(redshift)
pl13.H(redshift)
pl13.scale_factor(redshift)
0.7692307692307692
pl13.comoving_distance(redshift)
plt.plot((redshifts), vel_em(redshifts), label='Velovity now')
plt.plot((redshifts), vel_now(redshifts), label='Velocity at emission')
plt.axhline(con.c.to(kms).value, color='k', ls='--')
plt.legend(loc='best')
#plt.axis((-.1, 3, 0, 4e5))
#plt.plot(redshifts, vel_em(redshifts))
<matplotlib.legend.Legend at 0x7f8e66455810>
MatDom = LambdaCDM(0.7, 0.3, 0)
MatDom
LambdaCDM(H0=0.7 km / (Mpc s), Om0=0.3, Ode0=0, Tcmb0=2.725 K, Neff=3.04, m_nu=[ 0. 0. 0.] eV)