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, where we actually want to discover the mass distribution from the gas (or star) streaming.
If we assume the mass distribution is spherical, we know circular orbits are the simple solution to gas flow or periodic orbits in such a potential. This will allow us to predict the velocity field we are observing in galaxies such as NGC 6503.
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
In [ ]:
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 the balance of gravitation and centrifugal force in a spherical mass distribution: $$ { v^2 \over r } = {{ G M(<r) } \over r^2} $$ or $$ v = \sqrt{ {G M(<r) } \over r} $$
Of course this implies (and that's what we eventually want to do) that for a giving rotation curve, $v$, we can find out the mass distribution: $$ G M(<r) = v^2 r $$
In [ ]:
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 usually define the velocity to be 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/np.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 1900s description of clusters, and is also 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 quite 1 at 1 yet
# return radius / (1+radius*radius)**0.75
return radius / (0.5+0.5*radius*radius)**0.75
else:
return 0.0
In [ ]:
#model = 'star'
#model = 'galaxy'
model = 'plummer'
rad = np.arange(0.0,4.0,0.05)
vel = np.zeros(len(rad)) # this also works: vel = rad * 0.0
for i in range(len(rad)):
vel[i] = velocity(rad[i],model)
print("First, peak and Last value:",vel[0],vel.max(),vel[-1])
In [ ]:
plt.plot(rad,vel)
plt.xlabel("Radius")
plt.ylabel("Velocity")
plt.title("Rotation Curve (%s)" % model);
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 galaxies. However, for this experiment we will keep that assumption.
In [ ]:
# set the inclination of the disk with the line of sigh
inc = 60 # (0 means face-on, 90 means edge-on)
# some helper variables
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.1
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. Since we have a simple model here, let's take this approach.
The so-called forward projection we would need to use some extra steps that only add to the complexity.
In [ ]:
dr = 0.5
x = np.arange(-r0,r0,dr)
y = np.arange(-r0,r0,dr)
xx,yy = np.meshgrid(x,y) # helper variables for interpolations
rr = np.sqrt(xx*xx+(yy/cosi)**2)
if r0/dr < 20:
plt.scatter(xx,yy)
else:
print("not plotting too many gridpoints/dimension",r0/dr)
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.
You could also try and modify the velocity function so it takes a numpy array as input, and return a numpy !!!
In [ ]:
def velocity2d(rad2d, model):
""" convenient helper function to take a 2d array of radii
and return the same-shaped velocities
"""
(ny,nx) = rad2d.shape
vel2d = rad2d.copy() # could also do np.zeros(nx*ny).reshape(ny,nx)
for y in range(ny):
for x in range(nx):
vel2d[y,x] = velocity(rad2d[y,x],model)
return vel2d
In [ ]:
vv = velocity2d(rr,model)
vvmasked = np.ma.masked_where(rr>r0,vv)
vobs = vvmasked * xx / rr * sini
print("V_max:",vobs.max())
In [ ]:
vmax = 1
vmax = vobs.max()
if vmax > 0:
plt.imshow(vobs,origin=['Lower'],vmin=-vmax, vmax=vmax)
#plt.matshow(vobs,origin=['Lower'],vmin=-vmax, vmax=vmax)
else:
plt.imshow(vobs,origin=['Lower'])
plt.colorbar()
This plot will not look very good :