In [ ]:
import datetime
# Third-party
from astropy.io import ascii
import astropy.coordinates as coord
import astropy.units as u
import astropy.time as at
import matplotlib as mpl
import matplotlib.pyplot as pl
import numpy as np
%matplotlib inline
In [ ]:
import astroplan
from astroplan import Observer, FixedTarget
from astropy.time import Time
In [ ]:
mdm = Observer.at_site("MDM", timezone="America/Phoenix")
t1 = Time(datetime.datetime(2015, 10, 14, 18, 0))
t2 = t1 + 12*u.hour
time_range = Time([t1, t2])
In [ ]:
def coords_in_rect(c, corner_c):
if not c.frame.is_equivalent_frame(corner_c[0].frame):
raise ValueError("Frame mismatch.")
min_lon = corner_c[0].spherical.lon
min_lat = corner_c[0].spherical.lat
max_lon = corner_c[1].spherical.lon
max_lat = corner_c[1].spherical.lat
return ((c.spherical.lon > min_lon) & (c.spherical.lon < max_lon) &
(c.spherical.lat > min_lat) & (c.spherical.lat < max_lat))
In [ ]:
css = ascii.read("/Users/adrian/projects/triand-rrlyrae/data/catalina.csv")
print(css.colnames)
In [ ]:
linear = ascii.read("/Users/adrian/projects/triand-rrlyrae/data/linear.csv")
print(linear.colnames)
In [ ]:
css_c = coord.SkyCoord(ra=css['ra']*u.deg, dec=css['dec']*u.deg, distance=css['helio_dist']*u.kpc)
lin_c = coord.SkyCoord(ra=linear['ra']*u.deg, dec=linear['dec']*u.deg, distance=linear['helio_dist']*u.kpc)
In [ ]:
window_corners = [coord.SkyCoord(l=100*u.deg, b=-20*u.deg,frame='galactic'),
coord.SkyCoord(l=220*u.deg, b=-10*u.deg,frame='galactic')]
css_ix = coords_in_rect(css_c.galactic, window_corners)
lin_ix = coords_in_rect(lin_c.galactic, window_corners)
print("{} CSS targets, {} LINEAR targets in this window.".format(css_ix.sum(), lin_ix.sum()))
In [ ]:
gass_tbl = css[css_ix]
gass = coord.SkyCoord(l=gass_tbl['l']*u.deg, b=gass_tbl['b']*u.deg,
distance=gass_tbl['helio_dist']*u.kpc, frame='galactic')
In [ ]:
ix = (gass.distance > 2.5*u.kpc) & (gass.distance < 8.*u.kpc)
gass = gass[ix]
gass_tbl = gass_tbl[ix]
In [ ]:
print("{} GASS targets".format(len(gass)))
In [ ]:
pl.figure(figsize=(10,6))
pl.scatter(gass.l.degree, gass.b.degree, c=gass_tbl['VmagAvg'], marker='o')
pl.xlim(220,100)
pl.ylim(-20,-10)
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(12,5))
n,bins,pa = axes[0].hist(gass.distance, bins=np.linspace(2,10,15))
axes[0].set_xlabel("Radial dist. [kpc]")
axes[0].set_xlim(bins.min(), bins.max())
n,bins,pa = axes[1].hist(gass.galactocentric.represent_as(coord.CylindricalRepresentation).rho, bins=bins+8)
axes[1].set_xlabel("Cylindrical dist. [kpc]")
axes[1].set_xlim(bins.min(), bins.max())
axes[0].set_ylabel("Number RRL")
In [ ]:
fig,ax = pl.subplots(1,1,figsize=(6,6),subplot_kw =dict(polar=True))
ax.add_artist(mpl.patches.Circle((-8.,0), radius=2.5, transform=ax.transData._b, facecolor='r', alpha=0.2))
ax.add_artist(mpl.patches.Circle((-8.,0), radius=0.5, transform=ax.transData._b, facecolor='y', alpha=1.))
gass_cyl = gass.galactocentric.represent_as(coord.CylindricalRepresentation)
css_cyl = css_c.galactocentric.represent_as(coord.CylindricalRepresentation)
ax.plot(css_cyl.phi.to(u.radian)[np.abs(css_cyl.z) < 5*u.kpc],
css_cyl.rho.to(u.kpc)[np.abs(css_cyl.z) < 5*u.kpc],
color='k', linestyle='none', marker='.', alpha=0.1)
# ax.plot(gass_cyl.phi.to(u.radian), gass_cyl.rho.to(u.kpc),
# color='k', linestyle='none', marker='o')
ax.scatter(gass_cyl.phi.to(u.radian), gass_cyl.rho.to(u.kpc),
c=gass_tbl['VmagAvg'], marker='o')
ax.set_rmax(20.0)
ax.grid(True)
ticks = [5,10,15]
ax.set_rticks(ticks)
ax.set_yticklabels(['{0:d} kpc'.format(x) for x in ticks])
ax.set_xlabel("Galactic Longitude", labelpad=15)
ax.tick_params(axis='y', labelsize=20)
# fig.savefig("/Users/adrian/papers/proposals/MDM-2015/GASS.pdf")
# ------
fig,ax = pl.subplots(1,1,figsize=(7,6))
gass_cyl = gass.galactocentric.represent_as(coord.CylindricalRepresentation)
css_cyl = css_c.galactocentric.represent_as(coord.CylindricalRepresentation)
ax.plot(css_cyl.rho.to(u.kpc)[np.abs(css_cyl.z) < 5*u.kpc],
css_cyl.z.to(u.kpc)[np.abs(css_cyl.z) < 5*u.kpc],
color='k', linestyle='none', marker='.', alpha=0.25)
cc = ax.scatter(gass_cyl.rho.to(u.kpc),
gass_cyl.z.to(u.kpc),
c=gass_tbl['VmagAvg'], marker='o')
ax.set_xlim(0,25)
ax.set_ylim(-12.5,12.5)
fig.colorbar(cc)
# ticks = [10,20]
# ax.set_rticks(ticks)
# ax.set_yticklabels(['{0:d} kpc'.format(x) for x in ticks])
# ax.set_xlabel("Galactic Longitude", labelpad=15)
# ax.tick_params(axis='y', labelsize=20)
# # fig.savefig("/Users/adrian/papers/proposals/MDM-2015/GASS.pdf")
In [ ]:
pl.hist(gass_tbl['VmagAvg'])
In [ ]:
gass_tbl['ID2015'] = ["GASS2015RR{0:d}".format(x+1) for x in np.arange(len(gass_tbl)).astype(int)]
In [ ]:
# ascii.write(gass_tbl, "/Users/adrian/projects/triand-rrlyrae/data/targets/gass.txt")
In [ ]:
# ascii.write(gass_tbl[['ID2015','ra','dec']], "/Users/adrian/projects/triand-rrlyrae/data/targets/gass_targets_2015_short.txt")#, format="ascii")
In [ ]:
gass_tbl_sex = gass_tbl.copy()
ra = coord.Longitude(gass_tbl_sex['ra']*u.deg)
gass_tbl_sex['ra_sex'] = ra.to_string(unit=u.hour, precision=5, sep=' ')
dec = coord.Latitude(gass_tbl_sex['dec']*u.deg)
gass_tbl_sex['dec_sex'] = dec.to_string(unit=u.degree, precision=5, sep=' ')
ascii.write(gass_tbl_sex[['ID2015','ra_sex','dec_sex']], "/Users/adrian/projects/triand-rrlyrae/data/targets/gass_targets_2015_short_sexa.txt")
In [ ]:
window_corners = [coord.SkyCoord(l=100*u.deg, b=-35*u.deg,frame='galactic'),
coord.SkyCoord(l=160*u.deg, b=-15*u.deg,frame='galactic')]
css_ix2 = coords_in_rect(css_c.galactic, window_corners)
print("{} CSS targets in this window.".format(css_ix2.sum()))
In [ ]:
triand_tbl = css[css_ix2]
triand = coord.SkyCoord(l=triand_tbl['l']*u.deg, b=triand_tbl['b']*u.deg,
distance=triand_tbl['helio_dist']*u.kpc, frame='galactic')
In [ ]:
ix = (triand.distance > 15*u.kpc) & (triand.distance < 21*u.kpc)
triand = triand[ix]
triand_tbl = triand_tbl[ix]
In [ ]:
print("{} TriAnd targets".format(len(triand)))
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(12,5))
axes[0].hist(triand_tbl['helio_dist'], bins=np.linspace(15,25,12))
axes[0].set_xlabel("Radial dist. [kpc]")
axes[0].set_xlim(11,25)
# axes[1].hist(gc_cyl_triand.rho, bins=8)
# axes[1].set_xlabel("Cylindrical dist. [kpc]")
# axes[1].set_xlim(17,24)
axes[0].set_ylabel("Number RRL")
In [ ]:
# fig,ax = pl.subplots(1,1,figsize=(8,8),subplot_kw =dict(polar=True))
# ax.add_artist(mpl.patches.Circle((-8.,0), radius=0.5, transform=ax.transData._b, facecolor='y', alpha=1.))
# ax.plot(all_gc_cyl.phi.to(u.radian)[sky_window], all_gc_cyl.rho.to(u.kpc)[sky_window],
# color='k', linestyle='none', marker='o', alpha=0.25)
# ax.plot(gc_cyl_triand.phi.to(u.radian), gc_cyl_triand.rho.to(u.kpc),
# color='k', linestyle='none', marker='o')
# ax.set_rmax(30.0)
# ax.grid(True)
In [ ]:
# triand2013_tbl = ascii.read("/Users/adrian/projects/triand-rrlyrae/data/TriAnd_RRL_26mar15.csv")
triand2013_tbl = ascii.read("/Users/adrian/projects/triand-rrlyrae/data/publication_data.csv")
triand2013 = coord.SkyCoord(ra=triand2013_tbl['ra']*u.deg,
dec=triand2013_tbl['dec']*u.deg)
In [ ]:
pl.hist(triand_tbl['VmagAvg'])
Match to old data
In [ ]:
idx, sep2d, _ = triand.match_to_catalog_sky(triand2013.galactic)
In [ ]:
match_ix = (sep2d < 5*u.arcsec)
match_ix.sum()
In [ ]:
triand2013_tbl
In [ ]:
names = list()
old_names = list()
for i,x in enumerate(sep2d):
names.append("TriAnd2015RR{0:d}".format(i+1))
if x < 5*u.arcsec:
old_names.append(triand2013_tbl[idx[i]]['name'])
else:
old_names.append("--")
triand_tbl['ID2015'] = names
triand_tbl['ID2013'] = old_names
triand_tbl
In [ ]:
# ascii.write(triand_tbl, "/Users/adrian/projects/triand-rrlyrae/data/targets/triand1_targets_2015.txt")#, format="ascii")
In [ ]:
# ascii.write(triand_tbl[triand_tbl['ID2013']=='--'][['ID2015','ra','dec']], "/Users/adrian/projects/triand-rrlyrae/data/targets/triand1_targets_2015_short.txt")#, format="ascii")
In [ ]:
triand_tbl_sex = triand_tbl.copy()
ra = coord.Longitude(triand_tbl_sex['ra']*u.deg)
triand_tbl_sex['ra_sex'] = ra.to_string(unit=u.hour, precision=5, sep=' ')
dec = coord.Latitude(triand_tbl_sex['dec']*u.deg)
triand_tbl_sex['dec_sex'] = dec.to_string(unit=u.degree, precision=5, sep=' ')
ascii.write(triand_tbl_sex[triand_tbl['ID2013']=='--'][['ID2015','ra_sex','dec_sex']], "/Users/adrian/projects/triand-rrlyrae/data/targets/triand1_targets_2015_short_sexa.txt")
In [ ]:
from gary.observation import distance, apparent_magnitude
from gary.observation.rrlyrae import M_V
In [ ]:
distance(12.5 - M_V(-1.5)).to(u.kpc)
In [ ]:
apparent_magnitude(M_V(-1.5), 10.*u.kpc)
In [ ]:
from scipy.interpolate import InterpolatedUnivariateSpline
In [ ]:
pl.semilogy([14.5,16,17.5], [150,600,2400])
pl.xlim(14, 18)
In [ ]:
def Vmag_to_exptime(V):
s = InterpolatedUnivariateSpline([14.5,16,17.5], np.log10([150,600,2400]), k=1)
y = 10**s(V)
return y
In [ ]:
triand_exptimes = [Vmag_to_exptime(V) for V in triand_tbl['VmagAvg']]
GASS_exptimes = [Vmag_to_exptime(V) for V in gass_tbl['VmagAvg']]
In [ ]:
len(triand_exptimes), len(GASS_exptimes)
In [ ]:
# print("TriAnd", (sum(triand_exptimes)*u.second).to(u.hour))
print("GASS", (sum(GASS_exptimes)*u.second).to(u.hour))
In [ ]: