Python 1 solutions with `astropy`

Astropy is an open source collection of tools for astronomers written in Python. One of astropy's most useful features for using Python as a calculator is its units module. This module lets you assign units to variables and understands how they propagate through arithmetic.

5) Schwarzschild radii

a) Print the masses:



In [11]:

    
import numpy as np

# This is how we'll access astropy's units module: 
import astropy.units as u

# astropy stores fundamental physical constants like 
# the mass of the sun M_sun
from astropy.constants import M_sun

masses = u.Quantity([4.31e6*M_sun, 1*M_sun, 1.4*M_sun])
masses









    Out[11]:





$[8.573021 \times 10^{36},~1.9891 \times 10^{30},~2.78474 \times 10^{30}] \; \mathrm{kg}$

b) Retrieve constants from astropy's constants module



In [23]:

    
# Newton's gravitational constant G, 
# speed of light c, radius of sun R_sun
from astropy.constants import c, G, R_sun

# Convert the masses to units of grams using the `to` method
masses.to(u.gram)









    Out[23]:





$[8.573021 \times 10^{39},~1.9891 \times 10^{33},~2.78474 \times 10^{33}] \; \mathrm{g}$

c) Compute $r_s$, convert to cm



In [20]:

    
r_s = 2*G*masses/c**2
print "r_s in centimeters: {0}".format(r_s.to(u.cm))









    



r_s in centimeters: [  1.27320480e+12   2.95407147e+05   4.13570005e+05] cm

d) $r_s$ in units of solar radii:



In [19]:

    
print "r_s in solar radii: {0}".format(r_s/R_sun)









    



r_s in solar radii: [  1.83061130e+01   4.24735800e-06   5.94630120e-06]

e) $r_s$ in units of AU



In [25]:

    
print "r_s in AU: {0}".format(r_s.to(u.AU))









    



r_s in AU: [  8.51084842e-02   1.97467481e-08   2.76454473e-08] AU

f) Compare $r_s$ of Sag A* to orbits of solar system planets:



In [69]:

    
%matplotlib inline
import matplotlib.pyplot as plt

planet_names = ['Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune']
semimajor_axes = u.Quantity([0.38, 0.72, 1, 1.5, 5.2, 9.5, 19.1, 30.0], unit=u.AU)

print "Any planets interior to r_s?: {0}".format(any(r_s[0] > semimajor_axes))

fig, ax = plt.subplots()
ax.bar(range(len(semimajor_axes)), semimajor_axes.to(u.AU).value, 
       log=True, color='k')
ax.set_xticklabels(planet_names, ha='left')
ax.axhline(r_s[0].to(u.AU).value, color='r', label='$r_s$')
ax.legend()
ax.set_ylabel('Semimajor Axis [AU]');









    



Any planets interior to r_s?: False

g) How much bigger is neutron star radius than $r_s$?



In [73]:

    
radius_neutron_star = 10*u.km
print("how much bigger is neutron star's radius?: {0:.2f}x bigger"
      .format(float(radius_neutron_star/r_s[2])))









    



how much bigger is neutron star's radius?: 2.42x bigger

6) Planets

a) Semimajor axes in cm:



In [47]:

    
planet_names = ['Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto']
semimajor_axes = u.Quantity([0.38, 0.72, 1, 1.5, 5.2, 9.5, 19.1, 30.0, 40], unit=u.AU)
semimajor_axes.to(u.cm)









    Out[47]:





$[5.6847191 \times 10^{12},~1.0771047 \times 10^{13},~1.4959787 \times 10^{13},~2.2439681 \times 10^{13},~7.7790893 \times 10^{13},~1.4211798 \times 10^{14},~2.8573193 \times 10^{14},~4.4879361 \times 10^{14},~5.9839148 \times 10^{14}] \; \mathrm{cm}$

b-c) Orbital period, in years

Use numpy's np.pi variable to get a high-precision value for $\pi$



In [64]:

    
T = (4 * np.pi**2 * semimajor_axes**3 / (G*M_sun))**0.5
T.to(u.year)









    Out[64]:





$[0.23421959,~0.61086685,~0.99987985,~1.8368966,~11.8564,~29.277448,~83.463745,~164.29702,~252.95182] \; \mathrm{yr}$

d) Light travel time in hours, given

$$\Delta t = \frac{\Delta x}{c}$$



In [65]:

    
(semimajor_axes/c).to(u.hour)









    Out[65]:





$[0.052672727,~0.099800957,~0.13861244,~0.20791866,~0.72078469,~1.3168182,~2.6474976,~4.1583732,~5.5444976] \; \mathrm{h}$



In [ ]:

Python 1 solutions with astropy

5) Schwarzschild radii

6) Planets

Python 1 solutions with `astropy`