In [232]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import astropy.units as u
from astropy import constants
from astropy.io import fits
from astropy.coordinates import Galactic, Galactocentric, ICRS, LSR
from astropy.coordinates import (CartesianRepresentation,
CartesianDifferential)
%matplotlib inline
Test angular momentum from Sales+2017 and Fritz+2018 combo data.
In [233]:
_coord = Galactic(l=265.9*u.degree, b=-49.6*u.degree).transform_to(ICRS)
ra = _coord.ra
dec = _coord.dec
In [ ]:
In [234]:
coords_3d = ICRS(
ra=ra,
dec=dec,
distance=30*u.kpc,
pm_ra_cosdec=2.398*u.mas/u.yr*np.cos(ra.rad),
pm_dec=-1.319*u.mas/u.yr,
radial_velocity=62.8*u.km/u.s
)
In [235]:
coords_3d.transform_to(Galactocentric)
Out[235]:
First calculate the angular momentum direction given in Sales's paper:
In [236]:
LMC_sales = np.array([-0.97, 0.14, -0.18])
In [237]:
LMC_sales_r = np.linalg.norm(LMC_sales)
In [238]:
theta_LMC_sales = np.arctan2(LMC_sales[1], LMC_sales[0])
phi_LMC_sales = np.arccos(LMC_sales[2]/LMC_sales_r)
In [239]:
print(theta_LMC_sales, np.degrees(theta_LMC_sales))
print(phi_LMC_sales, np.degrees(phi_LMC_sales))
Read in a table of Milky Way dwarf masses.
In [250]:
dwarf_tbl = pd.read_table('/Users/spardy/Research/data/Auriga_data/McConnachie2012.tsv', delimiter=';')
Strip out any extra spaces to allow merging.
In [251]:
dwarf_tbl['Name'] = dwarf_tbl['Name'].apply(lambda e: e.strip())
Read in a table of positions and velocities.
In [252]:
df = pd.read_csv('/Users/spardy/Research/data/local_group_pos_vel_GaiaDR2_tableC4.txt')
Change Sagittarious's name to match with table above.
In [253]:
df.loc[7, 'Name'] = 'Sagittarius dSph'
Examine this table.
Merge the two tables together using the name column.
In [255]:
combined_df = pd.merge(dwarf_tbl, df, on='Name')
combined_df['Mass'] = combined_df['Mass'].astype(float)*1e+6
Examine the new table.
In [256]:
combined_df
Out[256]:
And just get the order of the dwarfs.
In [257]:
combined_df.Name
Out[257]:
Make vectors of position and velocity.
Offset for the position of the sun
In [258]:
combined_df['Y'] += 8.3
combined_df['Z'] += 27e-3
In [259]:
R = combined_df[['X', 'Y', 'Z']]
In [260]:
np.sqrt(np.sum(R**2, axis=1))
Out[260]:
In [261]:
np.linalg.norm(R - R.iloc[1], axis=1)
Out[261]:
And the solar motion.
In [262]:
combined_df['Vx'] = combined_df['U'] + 11.1
combined_df['Vy'] = combined_df['V'] + 232.24
combined_df['Vz'] = combined_df['W'] + 7.25
In [263]:
V = combined_df[['Vx', 'Vy', 'Vz']]
In [264]:
V
Out[264]:
Take the cross product of R and V
In [306]:
#L = np.cross(R.values, V.values*combined_df['Mass'].values[:, None])
L = np.cross(R.values, V.values)
And compute the norm of the angular momentum for each galaxy.
In [307]:
L_r = np.linalg.norm(L, axis=1) # 2-norm is default
Compute the unit angular momentum for each galaxy.
In [308]:
unit_L = L/L_r[:, None]
Add this unit vector into our dataframe
In [309]:
combined_df['Lx'] = unit_L[:, 0]
combined_df['Ly'] = unit_L[:, 1]
combined_df['Lz'] = unit_L[:, 2]
And examine the vectors.
In [310]:
combined_df[['Name', 'Lx', 'Ly', 'Lz']]
Out[310]:
In [110]:
galactic_coords = ICRS(
ra=combined_df['_RAJ2000'].values.astype('float')*u.degree,
dec=combined_df['_DEJ2000'].values.astype('float')*u.degree
).transform_to(Galactic)
In [111]:
combined_df['l'] = galactic_coords.l
combined_df['b'] = galactic_coords.b
Compute the first angle.
In [29]:
theta = np.arctan2(L[:, 1], L[:, 0])
And the second angle
In [30]:
phi = np.arccos(L[:, 2]/L_r)
Examine the angles in radians.
In [31]:
theta, phi
Out[31]:
And degrees.
In [32]:
np.degrees(theta), np.degrees(phi)
Out[32]:
We appear to be very close to the angle provided by Sales.
Offset the angles to make the LMC appear at the center of the Aitoff projection.
In [76]:
x = theta-np.pi
x[x < -np.pi] += 2*np.pi
y = phi-np.pi/2.
Look at smaller galaxies from Simon et al. 2018, using mass values from Sales
In [240]:
sales_masses = pd.read_csv('/Users/spardy/Research/data/Auriga_data/sales_18_tbl1_masses.csv')
simon_phase_space = pd.read_csv('/Users/spardy/Research/data/Auriga_data/simon_2018_tbl1and3.csv')
In [241]:
simon_coords = ICRS(
ra=simon_phase_space['RA'].values*u.deg,
dec=simon_phase_space['DEC'].values*u.deg,
distance=simon_phase_space['DIST'].values*u.kpc,
).transform_to(Galactocentric)
In [242]:
simon_phase_space['X'] = simon_coords.x.value
simon_phase_space['Y'] = simon_coords.y.value
simon_phase_space['Z'] = simon_coords.z.value
In [283]:
simon_phase_space['Vx'] = simon_phase_space['U'] + 11.1
simon_phase_space['Vy'] = simon_phase_space['V'] + 232.24
simon_phase_space['Vz'] = simon_phase_space['W'] + 7.25
Look at the smaller dwarfs from Fritz and Sales papers
In [244]:
fritz = pd.read_csv('/Users/spardy/Research/data/Auriga_data/fritz_2018_tbl2_selected.csv')
sales = pd.read_csv('/Users/spardy/Research/data/Auriga_data/sales_18_tbl1_selected.csv')
fritz_sales = pd.merge(left=fritz, right=sales, left_on='name', right_on='name')
Select only those that are not in the Simon table
In [281]:
fritz_sales = fritz_sales[~fritz_sales.name.isin(simon_phase_space.short_name)]
In [282]:
fritz_sales_ICRS = Galactic(
l=fritz_sales['l'].values*u.degree,
b=fritz_sales['b'].values*u.degree
).transform_to(ICRS)
fritz_sales_3d = ICRS(
ra=fritz_sales_ICRS.ra,
dec=fritz_sales_ICRS.dec,
distance=fritz_sales.d_sun.values*u.kpc,
pm_ra_cosdec=fritz_sales.mu_alpha.values*u.mas/u.yr*np.cos(fritz_sales_ICRS.ra),
pm_dec=fritz_sales.mu_delta.values*u.mas/u.yr,
radial_velocity=fritz_sales.v_sun.values*u.km/u.s
)
fritz_sales_names = fritz_sales['name']
fritz_sales_cartesian = fritz_sales_3d.transform_to(Galactocentric)
In [246]:
fritz_sales_V = np.stack((
fritz_sales_cartesian.v_x.value,
fritz_sales_cartesian.v_y.value,
fritz_sales_cartesian.v_z.value
), axis=1)
fritz_sales_R = np.stack((
fritz_sales_cartesian.x.value,
fritz_sales_cartesian.y.value,
fritz_sales_cartesian.z.value
), axis=1)
In [293]:
simon_selected = simon_phase_space[['short_name', 'X', 'Y', 'Z', 'Vx', 'Vy', 'Vz']]
In [299]:
uf_df = simon_selected.append(
pd.DataFrame({'short_name': fritz_sales_names,
'X': fritz_sales_cartesian.x.value,
'Y': fritz_sales_cartesian.y.value,
'Z': fritz_sales_cartesian.z.value,
'Vx': fritz_sales_cartesian.v_x.value,
'Vy': fritz_sales_cartesian.v_y.value,
'Vz': fritz_sales_cartesian.v_z.value})
).reset_index(drop=True)
In [300]:
uf_df
Out[300]:
In [311]:
uf_L = np.cross(uf_df[['X', 'Y', 'Z']], uf_df[['Vx', 'Vy', 'Vz']])
uf_L_r = np.linalg.norm(uf_L, axis=1)
In [312]:
uf_theta = np.arctan2(uf_L[:, 1], uf_L[:, 0])
uf_phi = np.arccos(uf_L[:, 2]/uf_L_r)
In [313]:
uf_x = uf_theta - np.pi
uf_x[uf_x < -np.pi] += 2*np.pi
uf_y = uf_phi - np.pi/2.
Plot this in aitoff coordinates.
In [314]:
associated = [7, 2, 6]
In [315]:
unassociated = [5, 8, 0, 9, 3, 4]
In [394]:
additional_x_offsets = {
'Seg 1': np.radians(-0.7),
'UM 2': np.radians(-0.9),
'Com': np.radians(-0.5),
'Boo 1': np.radians(0)
}
additional_y_offsets = {
'Seg 1': np.radians(2.2),
'UM 2': np.radians(2.2),
'Com': np.radians(3),
'Boo 1': np.radians(1)
}
In [396]:
fig, axis = plt.subplots(1, subplot_kw={'projection': 'aitoff'}, dpi=200)
axis.plot(x[unassociated], y[unassociated], '.', color='gray', ms=10)
axis.plot(x[associated], y[associated], '.', color='red', ms=10)
axis.plot(x[1], y[1], 'r^', ms=6)# LMC
axis.plot(uf_x, uf_y, '+', color='gray', ms=10)
for i, (name, x_i, phi_i) in enumerate(zip(combined_df.Name, x, y)):
# change some labels so that they don't overlap
if phi_i < np.radians(10):
va='top'
else:
va='bottom'
if name in ['Draco', 'Carina', 'Sculptor']:
ha = 'right'
else:
ha = 'left'
offset = -0.085 if va == 'top' else 0.085
text_color = 'black' if i in unassociated else 'red'
axis.text(x_i, phi_i+offset, name, fontsize='xx-small', ha=ha, va=va, color=text_color)
for i, (name, x_i, phi_i) in enumerate(zip(uf_df['short_name'], uf_x, uf_y)):
ha = 'right' if name not in ['Dra 2', 'Seg 2', 'Com',
'Boo 2', 'Tri 2', 'Hyd 1'] else 'left'
va = 'bottom' if name not in ['Seg 2', 'Boo 2', 'Car 3'] else 'top'
offset_y = 0.025 if ha == 'bottom' else -0.025
offset_x = -0.05 if ha == 'right' else 0.05
if name in additional_x_offsets:
offset_x += additional_x_offsets[name]
offset_y += additional_y_offsets[name]
axis.text(x_i+offset_x, phi_i+offset_y, f'{i+1}', fontsize='xx-small', ha=ha, va=va)
for i, name in enumerate(uf_df['short_name']):
axis.annotate(xy=(0.1*(i//2), -0.1-0.08*(i%2)), s=f'{i+1}: {name}',
xycoords='axes fraction', fontsize='xx-small')
axis.set_xticklabels(['%d$^{\circ}$' % num for num in np.arange(30, 360, 30)], fontsize='x-small')
axis.set_yticklabels(['%d$^{\circ}$' % num for num in np.arange(15, 180, 15)], fontsize='x-small')
axis.set_title('Directon of angular momentum vector\nGalactic coordinates', fontsize='xx-small')
axis.grid(True)
x0 = x[1]
y0 = y[1]
#(x - x0)^2 + (y - y0)^2 == r^2
for r in np.arange(10, 40, 10):
r = np.radians(r)
x_circ = np.linspace(x0 - r, x0 + r, 1000)
y_circ = np.sqrt(r**2 - (x_circ - x0)**2) + y0
y_circ = np.append(y_circ, -np.sqrt(r**2 - (x_circ[::-1] - x0)**2) + y0)
x_circ = np.append(x_circ, x_circ[::-1])
axis.plot(x_circ, y_circ, color='blue', linewidth=0.5)
fig.savefig('/Users/spardy/Research/projects/Auriga/for_paper/angular_momentum.pdf', bbox_inches='tight')
Compute the cosine distance between each galaxy's angular momentum vector.
In [37]:
cos_distance = np.dot(unit_L, unit_L[1])
Sort these.
In [38]:
s = np.argsort(cos_distance)
And plot them with the names appearing on the x-axis.
In [39]:
fig, axis = plt.subplots(1)
axis.bar(np.arange(len(s)), cos_distance[s])
axis.set_xticks(np.arange(len(s)))
_ = axis.set_xticklabels(combined_df.Name.iloc[s], rotation=90)
axis.set_title('Angular momentum vector similarity with LMC')
axis.set_ylabel('Cosine Similarity')
fig.tight_layout()
fig.savefig('/Users/spardy/Research/projects/Auriga/for_paper/angular_momentum_cosine.pdf')
In [ ]: