The gas streaming around a young star, or in a galactic disk is dominated by gravity. So we can simply compute the orbits of a point mass around a star, or in the more complex potential of a galactic disk. Although there is a wonderful package in python for this, galpy, we will first try and use this from scratch in simple pure python.
In [1]:
%matplotlib inline
In [3]:
# python 2-3 compatibility
from __future__ import print_function
In [4]:
import numpy as np
import math
First we need to define a function that tells us the speed of the gas at a given distance from the center of the star or galaxy. We consider only three simple cases here, always based on $$ { v^2 \over r } = {{ G. M(<r) } \over r^2} $$ or $$ v = \sqrt{ {G. M(<r) } \over r} $$
In [5]:
def velocity(radius, model='galaxy'):
"""describe the streaming velocity as function of radius in or around an object
such as a star or a galaxy. We define the velocity to be/data2/teuben/HASC-software/
1 at a radius of 1.
"""
if model == 'star':
# A star has a keplerian rotation curve. The planets around our sun obey this law.
if radius == 0.0:
return 0.0
else:
return 1.0/math.sqrt(radius)
elif model == 'galaxy':
# Most disk galaxies have a flat rotation curve with a linear slope in the center.
if radius > 1.0:
# flat rotation curve outside radius 1.0
return 1.0
else:
# solid body inside radius 1.0, linearly rising rotation curve
return radius
elif model == 'plummer':
# A plummer sphere was an early 1800s description of clusters, and is not
# a bad description for the inner portions of a galaxy. You can also view it
# as a hybrid and softened version of the 'star' and 'galaxy' described above.
# Note: not 1 at 1 yet
return radius / (1+radius*radius)**0.75
else:
return 0.0
In [6]:
model = 'star'
#model = 'galaxy'
model = 'plummer'
rad = np.arange(0.0,4.0,0.05)
vel = rad * 0.0
for i in range(len(rad)):
vel[i] = velocity(rad[i],model)
In [7]:
import matplotlib.pyplot as plt
In [8]:
plt.plot(rad,vel)
plt.xlabel("Radius")
plt.ylabel("Velocity")
# plt.show()
Out[8]:
This curve of velocity as function of radius is called a Rotation Curve, and extracting such a curve from an observation is crucial to understanding the mass distribution within a galaxy, or the mass of the young star at the center of the disk. We are assuming the gas is on circular orbits, which turns out is not always correct. For this experiment we will keep that assumption.
In [9]:
# set the inclination of the disk with the line of sigh (0 means face-on, 90 means edge-on)
inc = 60
cosi = math.cos(inc*math.pi/180.0)
sini = math.sin(inc*math.pi/180.0)
# radius of the disk, and steps in radius
r0 = 4.0
dr = 0.5
In [10]:
# get random points along a regular series of circles, which describes the model
# we assume the gas is streaming along those circles
radii = np.arange(0.0,r0,dr)
xobs = np.arange(0)
yobs = np.arange(0)
vobs = np.arange(0)
for r in radii:
vrot = velocity(r,model)
phi = np.random.random(45)*2*math.pi
# orbit properties we need
x = r*np.cos(phi)
y = r*np.sin(phi)
v = vrot*np.cos(phi)
# project
xobs = np.append(xobs,x)
yobs = np.append(yobs,y*cosi)
vobs = np.append(vobs,v*sini)
Now we have a set of (xobs,yobs,vobs) points that need to be gridded on the sky.
In [11]:
plt.scatter(xobs,yobs)
Out[11]:
In [12]:
plt.scatter(xobs,vobs)
Out[12]:
In [13]:
import scipy.interpolate
In [14]:
# f = scipy.interpolate.interp2d(xobs,yobs, vobs)
Perhaps you can already see some flaws in this approach. But it's a more general approach to simulate an observation. We will now do the reverse operation, since the rotation disk is simple to describe.
This is where we take a point in the sky, and deproject back where in the galaxy this point came from and compute the velocity and projected velocity. The big advantage is the simplicity of computing the observable at each picked point in the sky. The big drawback is that the deprojection may not be trivial in cases where the model is not simple, e.g. non-circular motion and/or non-planar disks.
In [15]:
dr = 0.05
x = np.arange(-r0,r0,dr)
y = np.arange(-r0,r0,dr)
xx,yy = np.meshgrid(x,y)
rr = np.sqrt(xx*xx+(yy/cosi)**2)
if r0/dr < 20:
plt.scatter(xx,yy)
Although we have defined a function velocity to compute the rotation velocity at any radius, this function cannot easily compute from a numpy array, as we just created on a grid on the sky. Thus we need a convenience function to do just that.
In [16]:
def velocity2d(rad2d, model):
""" convenience function to take a 2d array of radii
and return the same-shaped velocities
"""
(ny,nx) = rad2d.shape
vel2d = rad2d.copy()
for y in range(ny):
for x in range(nx):
vel2d[y,x] = velocity(rad2d[y,x],model)
return vel2d
In [19]:
model = 'star'
#model = 'plummer'
#model = 'galaxy'
vv = velocity2d(rr,model)
vvmasked = np.ma.masked_where(rr>r0,vv)
vobs = vvmasked * xx / rr * sini
print("max v:",vobs.max())
In [20]:
vmax = 1
if vmax > 0:
plt.imshow(vobs,origin=['Lower'],vmin=-vmax, vmax=vmax)
else:
plt.imshow(vobs,origin=['Lower'])
plt.colorbar()
Out[20]:
This plot will not look very good if you don't make the step size in radius (dr) small enough.
In [ ]: