Description
In this assignment you will determine the mass distribution of each galaxy at SnapNumber
0 and use this to determine each galaxy’s rotation curve.
• You will need the files “MW 000.txt”, “M31 000.txt” and “M33 000.txt”.
• Note, you will need to import Read from ReadFile and the class CenterOfMass from
the program CenterOfMass you made in Homework 4. So make sure you have .py files
of both codes in your working directory.
• Import relevant plotting modules – see solution to the Labs.
1 Class: Mass Profile
Create a class called MassProfile. You will follow the structure of the previous assignment
to initialize your class. But you will not specify the particle type at this stage. Also, instead
of the requiring the full filename as input, only supply the snapshot number and the galaxy
name. i.e. Inputs:
• galaxy: A String with Galaxy Name, e.g. “MW”, “M31” or “M33”
• snap: Snapshot number, e.g. 0, 1, etc
You will use these inputs to reconstruct the filename as follows:
# add a string of the filenumber to the value “000”
ilbl = ‘000’ + str(Snap)
# remove all but the last 3 digits
ilbl = ilbl[-3:]
self.filename=“%s ”%(galaxy) + ilbl + ‘.txt’
So if you want the file “MW 010.txt”, you would input “MW” and 10.
Next, read in the data for the x,y,z positions and mass. Store x,y,z, with appropriate units.
Don’t assign the mass units just yet.
Store the name of the galaxy as a global property self.gname.
1
2 Mass Enclosed
Create a function called MassEnclosed that will compute the mass enclosed within a given
radius of the COM position for a specified galaxy and a specified component of that galaxy.
The function:
• takes as input: the particle type and an array with radii.
• returns an array of masses (in units of M )
Hints:
• You will need to first determine the COM position by creating a CenterOfMass object
and calling COM P.
• You will need to loop over the Radius array to define particles that are enclosed within
the radius given at each array element. You can initialize this array using np.zeros.
3 Total Mass Enclosed
Create a function called MassEnclosedTotal that:
• Takes as input an array of radii
• Returns an array of masses (in units of M ) representing the total mass (bulge+disk+halo)
at each radius of the input array.
• NOTE: M33 does not have a bulge, so you need to add a caveat in the code to account
for this possibility, e.g. using self.gname, which you stored in step 1.
4 Hernquist Mass Profile
Create a function called HernquistMass that will compute the mass enclosed within a given
radius using the theoretical profile:
ρ(r) = M a
2πr
1
(r + a)
3 M(r) = Mhalor
2
(a + r)
2
(1)
This function takes as input: the radius, the scale factor a, the halo mass Mhalo. It
returns the halo mass in units of M .
5 Circular Velocity
Create a function called CircularVelocity with:
• inputs: the particle type and an array with radii.
2
• returns an array of circular speeds in units of km/s, rounded to two decimal places.
The circular speed is computed using the Mass enclosed at each radius, assuming spherical
symmetry.
Note that astropy stores constants. You can use it for e.g. the Gravitational constant.
from astropy.constants import G. But you will need to adjust the units. For our purposes:
G = G.to(u.kpc*u.km**2/u.s**2/u.Msun). With these units you will get V in km/s.
6 Total CircularVelocity
Create a function called CircularVelocityTotal that:
• Takes as input an array of radii
• Returns an array of circular velocity (in units of km/s) representing the total Vcirc
created by all the galaxy components (bulge+disk+halo) at each radius of the input
array.
NOTE: the total circular velocity is NOT just the circular velocity of each individual
galaxy component summed together.
7 Hernquist Circular Speed
Create a function called HernquistVCirc that computes the circular speed using the Hernquist
mass profile. You can either call HernquistMass or write out the actual equation.
This function takes as input: the radius, the scale factor a, the halo mass Mhalo. It
returns the circular speed in units of km/s, rounded to two decimal places.
8 Plot the Mass Profile for each Galaxy
1. Create a plot of the mass profile (mass enclosed as a function of radius) of each component of the MW to a radius of 30 kpc.
2. You will need to define an array of radii. Don’t start the first element at 0 (say 0.1
instead).
3. Use different line types and/or colors to indicate different components. If you get stuck
with the plotting look at the solutions for InClassLab1 on the class GitHub repo.
4. Y-axis should be in log (plt.semilogy).
5. Overplot the sum of all components (MassEnclosedTotal)
3
6. Determine the best fitting Hernquist profile:
Use the total dark matter mass of each galaxy as Mhalo.
You will need to find the scale radius (a) where the Hernquist profile best matches the
plotted mass profile as a function of radius. Do this by plotting the Hernquist profile
(starting with a guess for a) on top of the simulated dark matter halo mass profile.
Then keep changing a until the mass distributions reasonably agree. Alternatively, you
could print out values for the simulated enclosed mass at 30 kpc and change a until
you get a reasonable agreement with the Hernquist profile.
7. Use a different line-type or color for the best fit Hernquist profile and include text that
states the best fit scale length.
8. Repeat for M33 and M31.
9 Plot the Rotation Curve for each Galaxy
1. Create a plot of the Rotation Curve of each galaxy – i.e. the circular speed of each
galaxy component to a radius of 30 kpc.
2. You will need to define an array of radii. Don’t start the first element at 0 (say 0.1
instead).
3. Use different line types and/or colors to indicate different components.
4. Overplot the total circular speed for all galaxy component (CircularVelocityTotal).
5. Overplot the best fit theoretical Hernquist circular speed, using the scale radius you
determined from the mass enclosed plot.
6. Do the same for M31 and M33.
4