See also: https://ui.adsabs.harvard.edu/#abs/2015MNRAS.454.1453R/abstract
Downloaded PM GC catalogs from: http://www.stsci.edu/~marel/hstpromo.html#Projects
For dSph PMs: Piatek et al. 2003, 2005, 2006, 2007; Lepine et al. 2011; Pryor et al. 2014
In [ ]:
import os
# Third-party
from astropy.io import ascii
import astropy.coordinates as coord
import astropy.units as u
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
%matplotlib inline
import gary.coordinates as gc
import gary.dynamics as gd
import gary.potential as gp
from ophiuchus import galactocentric_frame, vcirc, vlsr
import ophiuchus.potential as op
from astropy.utils.data import download_file
In [ ]:
data_path = "/Users/adrian/projects/globmfer/data"
In [ ]:
import astropy.table as at
In [ ]:
pm_gc_main = np.genfromtxt(os.path.join(data_path,"gl_2012_J2000.cat1.txt"), dtype=None,
skip_header=2,
usecols=[0,2,3,6,7,8,9,10,11,12,13],
names=['ngc_num','ra','dec','dist','dist_err','mu_ra','mu_ra_err',
'mu_dec','mu_dec_err', 'vr', 'vr_err'])
pm_gc_main = at.Table(pm_gc_main)
In [ ]:
len(pm_gc_main)
In [ ]:
pm_gc_bulge = np.genfromtxt(os.path.join(data_path,"bulge_J2000.cat1.txt"), dtype=None,
skip_header=2,
usecols=[0,2,3,6,7,8,9,10,11,12,13],
names=['ngc_num','ra','dec','dist','dist_err','mu_ra','mu_ra_err',
'mu_dec','mu_dec_err', 'vr', 'vr_err'])
pm_gc_bulge = at.Table(pm_gc_bulge)
In [ ]:
go = ascii.read(os.path.join(data_path,"go97_table1.txt"))
In [ ]:
all_gc = at.vstack((pm_gc_main, pm_gc_bulge))
all_gc['name'] = np.array(["NGC {}".format(x) for x in all_gc['ngc_num']])
all_gc = at.join(all_gc, go, keys='name')
In [ ]:
c = coord.ICRS(ra=all_gc['ra']*u.degree,
dec=all_gc['dec']*u.degree,
distance=all_gc['dist']*u.kpc)
In [ ]:
xyz = c.transform_to(galactocentric_frame).cartesian.xyz
vxyz = gc.vhel_to_gal(c, pm=(all_gc['mu_ra']*u.mas/u.yr,
all_gc['mu_dec']*u.mas/u.yr),
rv=all_gc['vr']*u.km/u.s,
galactocentric_frame=galactocentric_frame,
vcirc=vcirc, vlsr=vlsr)
In [ ]:
pot = op.load_potential('static_mw')
# pot = op.load_potential('barred_mw_8')
In [ ]:
w0 = gd.CartesianPhaseSpacePosition(pos=xyz, vel=vxyz)
orbit = pot.integrate_orbit(w0, dt=-0.5, nsteps=12000)
In [ ]:
pers = u.Quantity([orbit[:,i].pericenter() for i in range(orbit.norbits)])
apos = u.Quantity([orbit[:,i].apocenter() for i in range(orbit.norbits)])
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(10,5))
bins = np.logspace(-1,2,16)
axes[0].hist(pers.value, bins=bins)
axes[0].set_xscale('log')
bins = np.logspace(0,3.,16)
axes[1].hist(apos.value, bins=bins)
axes[1].set_xscale('log');
axes[0].set_xlabel("Pericenter [kpc]")
axes[1].set_xlabel("Apocenter [kpc]")
axes[0].set_ylabel("N")
Phases of the clusters
In [ ]:
periods = u.Quantity([np.abs(orbit[:,i].estimate_period()) for i in range(orbit.norbits)])
In [ ]:
from scipy.signal import argrelmin
In [ ]:
phases = np.zeros(periods.size)
for i in range(orbit.norbits):
sph,_ = orbit[:,i].represent_as(coord.SphericalRepresentation)
r = sph.distance
idx, = argrelmin(r)
phase = (orbit.t - orbit.t[idx[0]]) / periods[i]
phases[i] = phase[0]
In [ ]:
pl.hist(phases, bins=np.linspace(0,1,16));
pl.xlabel("Orbital phase today")
In [ ]:
pers_xyz = np.zeros((3,len(pers)))
pers_xyz[0] = pers.value
pers_xyz = pers_xyz*u.kpc
mx = pot.mass_enclosed(pers_xyz)
pers_xyz = np.zeros((3,len(pers)))
pers_xyz[2] = pers.value
pers_xyz = pers_xyz*u.kpc
mz = pot.mass_enclosed(pers_xyz)
In [ ]:
rtide_x = pers.to(u.pc) * (all_gc['M'] / (3*mx))**(1/3.)
rtide_z = pers.to(u.pc) * (all_gc['M'] / (3*mz))**(1/3.)
In [ ]:
rtide = np.mean(np.vstack((rtide_x, rtide_z)).value*rtide_z.unit, axis=0)
err_rtide = np.std(np.vstack((rtide_x, rtide_z)).value*rtide_z.unit, axis=0)
In [ ]:
core_radius = all_gc['Rc']*u.pc
ratio = rtide / core_radius
In [ ]:
all_gc['peri_rtide_to_rcore'] = ratio
sort_idx = np.argsort(ratio)[:][:10]
# ratio[sort_idx]
In [ ]:
import astropy.coordinates as coord
import astropy.units as u
In [ ]:
coord.Galactic(l=255*u.degree, b=48*u.degree).transform_to(coord.ICRS)
In [ ]:
all_gc['name','peri_rtide_to_rcore','ra','dec'][sort_idx]
In [ ]:
all_gc[sort_idx[0:10]]
In [ ]:
ngc5897_orbit = orbit[:,sort_idx[0]]
In [ ]:
fig = ngc5897_orbit.plot()
In [ ]:
ngc5897_c_back,ngc5897_v_back = ngc5897_orbit[:100].to_frame(coord.ICRS, vcirc=vcirc, vlsr=vlsr, # orbit -50 Myr
galactocentric_frame=galactocentric_frame)
ngc5897_orbit_forw = pot.integrate_orbit(ngc5897_orbit[0], dt=0.5, nsteps=100)
ngc5897_c_forw,ngc5897_v_forw = ngc5897_orbit_forw.to_frame(coord.ICRS, vcirc=vcirc, vlsr=vlsr, # orbit +50 Myr
galactocentric_frame=galactocentric_frame)
There are 293039 SDSS stars with good photometry in the grey selection box
In [ ]:
sdss_select = mpl.patches.Rectangle((210,-25), width=20, height=10, alpha=0.1)
fig,ax = pl.subplots(1,1,figsize=(8,6))
ax.plot(ngc5897_c_back.ra.degree, ngc5897_c_back.dec.degree, ls='none')
ax.plot(ngc5897_c_forw.ra.degree, ngc5897_c_forw.dec.degree, ls='none')
ax.set_xlim(325,200)
ax.set_ylim(-45,0)
ax.xaxis.set_ticks(np.arange(200,330,10))
ax.add_patch(sdss_select)
pl.minorticks_on()
pl.grid()
Or, in RA - distance
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(8,6))
ax.plot(ngc5897_c_back.ra.degree, ngc5897_c_back.distance, ls='none')
ax.plot(ngc5897_c_forw.ra.degree, ngc5897_c_forw.distance, ls='none')
ax.set_ylabel("Helio. distance [kpc]")
ax.set_xlim(325,200)
# ax.set_ylim(-45,0)
ax.xaxis.set_ticks(np.arange(200,330,10))
pl.minorticks_on()
pl.grid()
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(8,6))
ax.plot(ngc5897_c_back.ra.degree, ngc5897_v_back[2].to(u.km/u.s).value, ls='none')
ax.plot(ngc5897_c_forw.ra.degree, ngc5897_v_forw[2].to(u.km/u.s).value, ls='none')
ax.set_ylabel(r"$v_{\rm los}$ [km/s]")
ax.set_xlim(325,200)
# ax.set_ylim(-45,0)
ax.xaxis.set_ticks(np.arange(200,330,10))
pl.minorticks_on()
pl.grid()
Here I plot the positions of all stars selected from SDSS in the grey box:
In [ ]:
ngc5897_sdss = ascii.read("/Users/adrian/Downloads/NGC5897_adrn.csv")
In [ ]:
sdss_select = mpl.patches.Rectangle((210,-25), width=20, height=20, alpha=0.1)
fig,ax = pl.subplots(1,1,figsize=(8,6))
ax.plot(ngc5897_c_back.ra.degree, ngc5897_c_back.dec.degree, ls='none')
ax.plot(ngc5897_c_forw.ra.degree, ngc5897_c_forw.dec.degree, ls='none')
ax.plot(ngc5897_sdss['ra'], ngc5897_sdss['dec'], marker=',', ls='none', alpha=0.25)
ax.set_xlim(325,200)
ax.set_ylim(-45,0)
ax.xaxis.set_ticks(np.arange(200,330,10))
ax.add_patch(sdss_select)
pl.minorticks_on()
pl.grid()
In [ ]:
control_idx = (ngc5897_sdss['ra'] < 220.) & (ngc5897_sdss['dec'] < -10.) & (ngc5897_sdss['dec'] > -15.)
targets_idx = (ngc5897_sdss['dec'] < -15.) & (ngc5897_sdss['dec'] > -20.)
control = ngc5897_sdss[control_idx]
targets = ngc5897_sdss[targets_idx]
In [ ]:
# isochrone
ngc5897_iso = ascii.read("/Users/adrian/Downloads/ngc5897_iso.txt", header_start=8)
ngc5897_iso.colnames
Distance modulus = 15.55 mag (from http://arxiv.org/pdf/1403.1262v1.pdf)
In [ ]:
dm = 15.55
In [ ]:
_u = control['dered_u']
_g = control['dered_g']
_r = control['dered_r']
_i = control['dered_i']
alpha = 0.02
fig,axes = pl.subplots(1,3,figsize=(10,5),sharey=True)
axes[1].set_title("Control field", fontsize=20)
axes[0].plot(_u-_g, _g, ls='none', marker='.', alpha=alpha)
axes[1].plot(_g-_r, _g, ls='none', marker='.', alpha=alpha)
axes[2].plot(_g-_i, _g, ls='none', marker='.', alpha=alpha)
axes[0].plot(ngc5897_iso['sdss_u'] - ngc5897_iso['sdss_g'], ngc5897_iso['sdss_g'] + dm, marker=None)
axes[1].plot(ngc5897_iso['sdss_g'] - ngc5897_iso['sdss_r'], ngc5897_iso['sdss_g'] + dm, marker=None)
axes[2].plot(ngc5897_iso['sdss_g'] - ngc5897_iso['sdss_i'], ngc5897_iso['sdss_g'] + dm, marker=None)
axes[0].set_ylim(22,14)
axes[0].set_xlim(0.5,2.5)
axes[1].set_xlim(-0.25,1.1)
axes[2].set_xlim(-0.5,2)
axes[0].set_ylabel("$g_0$")
axes[0].set_xlabel("$(u-g)_0$")
axes[1].set_xlabel("$(g-r)_0$")
axes[2].set_xlabel("$(g-i)_0$")
fig.tight_layout()
In [ ]:
_u = targets['dered_u']
_g = targets['dered_g']
_r = targets['dered_r']
_i = targets['dered_i']
alpha = 0.02
fig,axes = pl.subplots(1,3,figsize=(10,5),sharey=True)
axes[1].set_title("NGC 5897 field", fontsize=20)
axes[0].plot(_u-_g, _g, ls='none', marker='.', alpha=alpha)
axes[1].plot(_g-_r, _g, ls='none', marker='.', alpha=alpha)
axes[2].plot(_g-_i, _g, ls='none', marker='.', alpha=alpha)
axes[0].plot(ngc5897_iso['sdss_u'] - ngc5897_iso['sdss_g'], ngc5897_iso['sdss_g'] + dm, marker=None)
axes[1].plot(ngc5897_iso['sdss_g'] - ngc5897_iso['sdss_r'], ngc5897_iso['sdss_g'] + dm, marker=None)
axes[2].plot(ngc5897_iso['sdss_g'] - ngc5897_iso['sdss_i'], ngc5897_iso['sdss_g'] + dm, marker=None)
axes[0].set_ylim(22,14)
axes[0].set_xlim(0.5,2.5)
axes[1].set_xlim(-0.25,1.1)
axes[2].set_xlim(-0.5,2)
axes[0].set_ylabel("$g_0$")
axes[0].set_xlabel("$(u-g)_0$")
axes[1].set_xlabel("$(g-r)_0$")
axes[2].set_xlabel("$(g-i)_0$")
fig.tight_layout()
In [ ]:
from scipy.misc import logsumexp
def matched_filter(data, isochrone, dm, smooth=0.02, threshold=2):
data = data[(data['dered_g'] > 14) & (data['dered_g'] < 21.5)]
isochrone = isochrone[((isochrone['sdss_g'] + dm) > 14) & ((isochrone['sdss_g'] + dm) < 21.5)]
i_g = isochrone['sdss_g'] + dm
i_ug = isochrone['sdss_u'] - isochrone['sdss_g']
i_gr = isochrone['sdss_g'] - isochrone['sdss_r']
i_gi = isochrone['sdss_g'] - isochrone['sdss_i']
d_g = data['dered_g']
d_ug = data['dered_u'] - data['dered_g']
d_gr = data['dered_g'] - data['dered_r']
d_gi = data['dered_g'] - data['dered_i']
d_g_var = data['err_g']**2
d_ug_var = data['err_u']**2 + data['err_g']**2
d_gr_var = data['err_g']**2 + data['err_r']**2
d_gi_var = data['err_g']**2 + data['err_i']**2
s_var = smooth**2
const = 0.5*np.log(2*np.pi*s_var)
dist = -0.5 * (d_g[None] - i_g[:,None])**2 / (s_var + d_g_var[None]) - const
dist += -0.5 * (d_ug[None] - i_ug[:,None])**2 / (s_var + d_ug_var[None]) - const
dist += -0.5 * (d_gr[None] - i_gr[:,None])**2 / (s_var + d_gr_var[None]) - const
dist += -0.5 * (d_gi[None] - i_gi[:,None])**2 / (s_var + d_gi_var[None]) - const
log_prob = logsumexp(dist, axis=0) - np.log(len(isochrone))
pl.hist(log_prob,bins=np.linspace(-1,12,16))
ix = log_prob > threshold
return data[ix]
In [ ]:
filter_dm = 15.2
filter_smooth = 0.01 # mag
filter_threshold = 9.5
In [ ]:
filtered_control = matched_filter(control, ngc5897_iso,
dm=filter_dm, smooth=filter_smooth,
threshold=filter_threshold)
In [ ]:
filtered_targets = matched_filter(targets, ngc5897_iso,
dm=filter_dm, smooth=filter_smooth,
threshold=filter_threshold)
In [ ]:
print(len(filtered_control), len(filtered_targets))
In [ ]:
In [ ]:
alpha = 0.02
fig,axes = pl.subplots(1,3,figsize=(10,5),sharey=True)
axes[1].set_title("NGC 5897 field", fontsize=20)
_u = targets['dered_u']
_g = targets['dered_g']
_r = targets['dered_r']
_i = targets['dered_i']
axes[0].plot(_u-_g, _g, ls='none', marker='.', alpha=alpha)
axes[1].plot(_g-_r, _g, ls='none', marker='.', alpha=alpha)
axes[2].plot(_g-_i, _g, ls='none', marker='.', alpha=alpha)
_u = filtered_targets['dered_u']
_g = filtered_targets['dered_g']
_r = filtered_targets['dered_r']
_i = filtered_targets['dered_i']
axes[0].plot(_u-_g, _g, ls='none', marker='.')
axes[1].plot(_g-_r, _g, ls='none', marker='.')
axes[2].plot(_g-_i, _g, ls='none', marker='.')
axes[0].plot(ngc5897_iso['sdss_u'] - ngc5897_iso['sdss_g'], ngc5897_iso['sdss_g'] + dm, marker=None)
axes[1].plot(ngc5897_iso['sdss_g'] - ngc5897_iso['sdss_r'], ngc5897_iso['sdss_g'] + dm, marker=None)
axes[2].plot(ngc5897_iso['sdss_g'] - ngc5897_iso['sdss_i'], ngc5897_iso['sdss_g'] + dm, marker=None)
axes[0].set_ylim(22,14)
axes[0].set_xlim(0.5,2.5)
axes[1].set_xlim(-0.25,1.1)
axes[2].set_xlim(-0.5,2)
axes[0].set_ylabel("$g_0$")
axes[0].set_xlabel("$(u-g)_0$")
axes[1].set_xlabel("$(g-r)_0$")
axes[2].set_xlabel("$(g-i)_0$")
fig.tight_layout()
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(8,6))
# ax.plot(ngc5897_c_back.ra.degree, ngc5897_c_back.dec.degree, ls='none')
# ax.plot(ngc5897_c_forw.ra.degree, ngc5897_c_forw.dec.degree, ls='none')
ax.plot(filtered_control['ra'], filtered_control['dec'], marker='.', ls='none', alpha=0.5)
ax.plot(filtered_targets['ra'], filtered_targets['dec'], marker='.', ls='none', alpha=0.5)
ax.set_xlim(220,210)
ax.set_ylim(-20,-10)
pl.minorticks_on()
pl.grid()
Hm, well, the above isn't very conclusive...maybe I'll try looking in teh catalina data too?
In [ ]:
css = ascii.read("/Users/adrian/Downloads/catalina_south.csv")
print(css.colnames)
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(8,6))
galcen = coord.Galactic(l=0*u.deg, b=0*u.deg).transform_to(coord.ICRS)
ax.scatter(galcen.ra.degree, galcen.dec.degree)
ax.plot(css['RAJ2000'], css['DEJ2000'], marker='.', ls='none')
ax.plot(ngc5897_c_back.ra.degree, ngc5897_c_back.dec.degree, ls='none')
ax.plot(ngc5897_c_forw.ra.degree, ngc5897_c_forw.dec.degree, ls='none')
ax.set_xlim(325,200)
ax.set_ylim(-45,0)
ax.xaxis.set_ticks(np.arange(200,330,10))
pl.minorticks_on()
pl.grid()
In [ ]:
box_ix = ((css['RAJ2000'] > 210) & (css['RAJ2000'] < 230) &
(css['DEJ2000'] > -22) & (css['DEJ2000'] < -15) &
(css['dH'] > 9.) & (css['dH'] < 13.))
box_ix.sum()
In [ ]:
pl.plot(css['RAJ2000'][box_ix], css['DEJ2000'][box_ix], ls='none', marker='o')
In [ ]: