In [12]:
#!/usr/bin/python
import sys
sys.path.insert(0, '../Newby-tools/utilities')
import astro_coordinates as ac
import numpy as np
import math as ma
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter, MaxNLocator
import scipy as sc
deg = 180.0 / ma.pi
rad = ma.pi / 180.0
lCP = 122.932 * rad #123.932 is the wrong value in B&M
surveyCenterRa = 185.0
surveyCenterDec = 32.5
raGP = 192.8594813 * rad
decGP = 27.1282511 * rad
arr = sc.array([0.0])
dsun = 8.5 # 8.0 is current IAU standard!
def checkGlobularClusters():
clusters = [ [], [], [], [] ]
f = open("/home/weissj3/Desktop/HarrisCatalog.txt", 'r')
#File Format 46-51 and 53-58
for line in f:
if line[46:51] != " " and line[53:58] != " " and line[60:64] != " ":
templ, tempb, tempr = float(line[46:51]), float(line[53:58]), float(line[60:64])
if templ > 40 and templ < 210 and tempb < -30 and tempb > -80:
clusters[0].append(templ)
clusters[1].append(tempb)
clusters[2].append(getg(tempr))
clusters[3].append(line[1:9])
print(line[1:9], templ, tempb, tempr)
return clusters
def getg (r, M=4.2):
#"""converts a distance (kpc) into a magnitude"""
return M + 5.*(sc.log10(r*1000) - 1.)
def readStarFile_lb(x):
f = open(x, 'r');
f.readline()
stars2 = np.array([ list(map(float, ln.replace(',', '').split())) for ln in f ])
return stars2[:,0].tolist(), stars2[:,1].tolist(), stars2[:,2].tolist()
def lb2GC(l, b, wedge):
ra, dec = lbToEq(l, b)
return EqToGC(ra, dec, wedge)
def lbToEq (l_deg, b_deg):
#""" Converts galactic l,b in to Equatorial ra, dec; from Binney and Merrifield, p. 31;
#l, b must be arrays of same shape"""
l, b = (l_deg*rad), (b_deg*rad)
# Conversion Code
t = lCP - l
dec = sc.arcsin(sc.sin(decGP)*sc.sin(b) + sc.cos(decGP)*sc.cos(b)*sc.cos(t) )
r = sc.arctan2( (sc.cos(b)*sc.sin(t)),
( (sc.cos(decGP)*sc.sin(b)) - (sc.sin(decGP)*sc.cos(b)*sc.cos(t))) )
if type(r) != type(arr): r = sc.array([r])
for i in range(len(r)):
r[i] = angle_bounds((r[i] + raGP)*deg)
return r, (dec*deg)
def angle_bounds (angle, min=0.0, max=360.0):
#""" Keeps an angle, in degrees, in a 360 degree region"""
while angle < min: angle = angle + 360.0
while angle > max: angle = angle - 360.0
return angle
def angle_bounds2 (theta, phi):
#""" Sets two spherical angles in bounds, -90<theta<90; 0<phi<360"""
if type(theta) != type(arr): theta = sc.array([theta])
if type(phi) != type(arr): phi = sc.array([phi])
for i in range(len(theta)):
theta[i] = angle_bounds(theta[i], -180.0, 180.0)
for i in range(len(theta)):
if sc.fabs(theta[i]) > 90.0:
theta[i] = 180.0 - theta[i]
phi[i] = phi[i] + 180.0
for i in range(len(theta)):
theta[i] = angle_bounds(theta[i], -180.0, 180.0)
phi[i] = angle_bounds(phi[i], 0.0, 360.0)
for i in range(len(theta)):
if (sc.fabs(theta[i]) == 90.0): phi[i] = 0.0
if len(theta)==1:
theta, phi = theta[0], phi[0]
return theta, phi
def get_eta (wedge):
#""" Get the eta value that corresponds to the given stripe value """ #wedge_eta?
ss = 2.5
if wedge <= 46: eta = wedge * ss - 57.5
else: eta = wedge * ss - 57.5 - 180.0
return eta
def EqToGC (ra_deg, dec_deg, wedge): #produces lists... anglebounds2!!! ROTATE
#""" Converts equatorial ra,dec into Great Circle mu, nu; 'atSurveyGeometry.c' in
#m31.phys.rpi.edu:/p/prd/astrotools/v5_18/Linux-2-4-2-3-2/src"""
node = (surveyCenterRa - 90.0)*rad
eta = get_eta(wedge)
inc = (surveyCenterDec + eta)*rad
ra, dec = (ra_deg*rad), (dec_deg*rad)
# Rotation
x1 = sc.cos(ra-node)*sc.cos(dec)
y1 = sc.sin(ra-node)*sc.cos(dec)
z1 = sc.sin(dec)
x2 = x1
y2 = y1*sc.cos(inc) + z1*sc.sin(inc)
z2 = -y1*sc.sin(inc) + z1*sc.cos(inc)
mu = sc.arctan2(y2,x2) + node
nu = sc.arcsin(z2)
nu, mu = angle_bounds2((nu*deg), (mu*deg))
return mu,nu
def makeArrows(data, coord = 'dec'):
sgrArrow = {'ra' : [], 'dec' : [], 'r' : [], 'norm' : []}
for k in range(len(data['mu'])):
tmpra, tmpdec, tmpr = ac.streamToEqR(0.0,0.0,1.0,data['mu'][k], data['r'][k], data['theta'][k]*ac.deg, data['phi'][k]*ac.deg, data['wedge'][k])
sgrArrow['ra'][k].append(tmpra[0] - data['ra'][k])
sgrArrow['dec'][k].append(tmpdec - data['dec'][k])
sgrArrow['r'][k].append(tmpr - data['r'][k])
sgrArrow['norm'][k].append(ma.sqrt(sgrArrow['ra'][k]**2. + sgrArrow['dec'][k]**2. + sgrArrow['r'][k]**2.))
for k in range(len(SgrRA)):
plt.arrow(data['ra'][k], data[coord][k], sgrArrow['ra'][k]/sgrArrow['norm'][k], sgrArrow[coord][k]/sgrArrow['norm'][k], length_includes_head=True, head_width=.1, head_length=.4)
In [5]:
stars = [ [], [], [] ]
for i in range(80, 87):
#temp1, temp2, temp3 = readStarFile_lb("/home/weissj3/Desktop/SDSSSouth/DR_14_Stripe_%d_Rev_4.stars" % i)
#temp1, temp2, temp3 = readStarFile_lb("/home/weissj3/Desktop/SDSSSouth/stars-%d.txt" % i)
temp1, temp2, temp3 = readStarFile_lb("/home/weissj3/Desktop/SDSSSouth/stars-%d_2.txt" % i)
#temp = np.where((np.array(stars[2]) < 20.0) & (np.array(stars[2]) > 10.0))# and (stars[2] > 24.0))
temp = np.where((np.array(temp2) < -30.0))# and (stars[2] > 24.0))
temp1 = list(np.array(temp1)[temp])
temp2 = list(np.array(temp2)[temp])
temp3 = list(np.array(temp3)[temp])
stars[0] = stars[0] + temp1
stars[1] = stars[1] + temp2
stars[2] = stars[2] + temp3
clusters = checkGlobularClusters()
In [6]:
from astropy.coordinates import SkyCoord
import astropy.units as u
from astropy.io import ascii
Mask = dict()
Mask["Hole"] = ascii.read("../Data/SDSSSkyHoleMask.csv")
Mask["Bright"] = ascii.read("../Data/SDSSSkyBrightMask.csv")
Mask["Bleed"] = ascii.read("../Data/SDSSSkyBleedMask.csv")
Mask["Trail"] = ascii.read("../Data/SDSSSkyTrailMask.csv")
c_low = dict()
for i in Mask:
c = SkyCoord(Mask[i]["ra"], Mask[i]["dec"], frame='icrs', unit='deg').transform_to('galactic')
c_mask = c.b < (-30 * u.deg)
c_low[i] = c[c_mask]
plt.plot(c_low["Hole"].l, c_low["Hole"].b, 'o')
plt.show()
In [4]:
plt.figure(1,figsize=(15,5))
plt.plot(stars[0], stars[1], "bo", ms=.05, alpha=.5)
plt.plot(clusters[0], clusters[1], "go")
plt.plot(c_low["Trail"].l, c_low["Trail"].b, 'go', ms=.4, alpha=.5)
plt.show()
In [5]:
plt.figure(1, figsize=(16, 8))
for i in range(79, 87):
plt.subplot(2,4,i-78)
temp1, temp2, temp3 = readStarFile_lb("/home/weissj3/Desktop/SDSSSouth/DR_14_Stripe_%d_Rev_4.stars" % i)
temp = np.where((np.array(temp2) < -30.0))# and (stars[2] > 24.0))
temp1 = list(np.array(temp1)[temp])
temp2 = list(np.array(temp2)[temp])
temp3 = list(np.array(temp3)[temp])
temp1, temp2 = lb2GC(np.array(temp1), np.array(temp2), i)
for j in range(len(temp1)):
if temp1[j] < 150:
temp1[j] = temp1[j] + 360.0
plt.plot(temp1, temp3, 'o', ms=0.1, alpha=0.5)
print(str(i) + "Mu lims")
print(min(temp1), max(temp1))
print("g Lim")
print(getg(min(temp3)), getg(max(temp3)))
print("Nu Lim")
print(min(temp2), max(temp2))
plt.title("Stripe %d" % i)
if i - 78 > 4:
plt.xlabel("Mu")
if (i - 78) == 1 or (i - 78) == 5:
plt.ylabel("R")
plt.show()
In [7]:
fitParameters1 = [[0.99610854796299, 0.602707695537048, -1.91409127023908, 394.999536157082, 26.3034791728196, 2.70612761879803, 2.85144592613987, 1.3969235362309, -0.23656017282301, 385.32961593779, 23.8400908073428, 0.4521374080478, 1.01738245820873, 3.79651317783094, -1.89765799828133, 365.013719108822, 45.3828724029551, 1.41743008998293, 0.675894712848845, 3.30838010526825],
[0.999944612079743, 0.305901485053103, -1.61861075316423, 397.112028958749, 26.9669879667783, 1.78816018208519, 0.237181579753778, 0.698941446484881, -0.15027266337857, 377.224029474032, 35.0499388278361, 1.09124117905718, 0.325627431055257, 24.9992751352876, 0.412951991703798, 392.385461143937, 24.1080968204511, 0.827109082148601, 0.537533461647184, 4.86748749793852],
[0.98991378963277, 0.308289981214212, -1.54382724939751, 363.312100226944, 39.8615768517978, 0.928860509563496, 0.0736519965323057, 16.7270141001641, -0.5446031371613, 399.869256268919, 24.8882388994142, 0.635058202623914, 1.81890872116164, 6.22180674708641, -4.56421789620042, 391.056876195451, 7.99796529512058, 0.73136614238237, 2.86445727474234, 22.5991697823511],
[0.99994759953964, 0.252833477693611, 0.183960129525568, 340.141570004638, 12.4598340724202, 1.67456004755789, 1.30164830039372, 6.84689873624907, -0.401329375713388, 362.457389587379, 40.9237849950465, 1.29500975831299, 0.41331563896243, 10.6841892946136, 0.632459559003521, 389.959052282875, 24.3200903237175, 1.72972116118681, 2.8132762718471, 4.18107368079411],
[0.999606274155892, 0.283807028588464, 0.884551229674617, 384.758080221005, 22.9512072305092, 1.74559499366151, 2.70674029862282, 5.05004302665926, -0.52962588812335, 399.92689633255, 67.0978032693911, 1.24931282938646, 0.608042133517162, 15.4413182362025, -0.0523015479247553, 340.000714861477, 2.15252890832674, 2.37347256036285, 0.557351411140675, 5.56553078066461],
[0.997768949410222, 0.506205313961172, -0.424565531941287, 340.008064867463, 26.7863877685715, 1.80776977585074, 0.382493321582206, 24.9966498350949, -1.72119765195003, 382.728981125017, 26.6522293241025, 2.58627633073296, 0.62934242535028, 1.05417274142966, -0.168595254298105, 374.119447132412, 20.1852846232421, 1.89900401819141, 3.13997610064033, 3.2405824698739],
[0.998890926635154, 0.396327451595263, 0.472434509625133, 378.937922438113, 21.7590460985581, 1.81776806778579, 2.68056020974052, 4.78316538013624, -0.833958369362128, 344.432868891238, 16.6107690999911, 0.376441773644166, 2.84831751984169, 3.84391483141543, -0.978129743249635, 340.205371874273, 57.8836148531743, 1.04479216936567, 0.573387147455888, 21.2988225130648]]
fitParameters2 = [[0.99610854796299, 0.602707695537048, -1.91409127023908, 394.999536157082, 26.3034791728196, 2.70612761879803, 2.85144592613987, 1.3969235362309, -0.23656017282301, 385.32961593779, 23.8400908073428, 0.4521374080478, 1.01738245820873, 3.79651317783094, -1.89765799828133, 365.013719108822, 45.3828724029551, 1.41743008998293, 0.675894712848845, 3.30838010526825],
[0.999944612079743, 0.305901485053103, -1.61861075316423, 397.112028958749, 26.9669879667783, 1.78816018208519, 0.237181579753778, 0.698941446484881, -0.15027266337857, 377.224029474032, 35.0499388278361, 1.09124117905718, 0.325627431055257, 24.9992751352876, 0.412951991703798, 392.385461143937, 24.1080968204511, 0.827109082148601, 0.537533461647184, 4.86748749793852],
[0.98991378963277, 0.308289981214212, -1.54382724939751, 363.312100226944, 39.8615768517978, 0.928860509563496, 0.0736519965323057, 16.7270141001641, -0.5446031371613, 399.869256268919, 24.8882388994142, 0.635058202623914, 1.81890872116164, 6.22180674708641, -4.56421789620042, 391.056876195451, 7.99796529512058, 0.73136614238237, 2.86445727474234, 22.5991697823511],
[0.998433175667497, 0.44604689398518, 0.143675394294249, 389.52394314553, 23.8289292057711, 1.87353948260878, 2.82652333236972, 3.33297683358306, -2.04618986982103, 356.772369133057, 15.8960758430977, 0.451143105962654, 0.86870402768177, 2.15922260769236, -0.0876352409689074, 365.072139074424, 18.8431691892643, 0.731447574817498, 0.97779694432965, 24.1448451689983],
[0.999436532700952, 0.359163875586857, 0.484484658413185, 382.18559487859, 21.1662645858236, 1.84341873896455, 2.99779926008825, 4.41959389621051, -0.266884172977324, 346.36239981359, 17.1387853364029, 2.46685582424965, 2.78176016830994, 21.1059483160892, -1.19401126823636, 385.440402805411, 27.3078224365302, 2.55620809469092, 0.531497802808212, 1.09846140783909],
[0.996420608418262, 0.566012322728858, 0.22291579757387, 379.216675961992, 21.6612287120313, 1.94190378995657, 2.99409932944404, 3.38732441365857, -0.689382556697996, 354.398289008745, 31.5919805062695, 1.28000761202547, 3.11272719317788, 23.0911827799794, -3.88767169457418, 375.865412697147, 38.3064561785442, 0.622246528900129, 1.94212105866804, 0.364013294273391],
[0.998051890513442, 0.483944937834668, 0.537046366593952, 372.294448057594, 20.9250162746688, 1.87364820818759, 2.94873748959537, 4.6103832068514, -1.83108880535357, 386.461857784932, 21.5631473137614, 2.39424656574779, 1.78269207444944, 1.8899303925527, -1.30457094133997, 347.223893684292, 61.9839144108112, 1.05361415595827, 0.31855690926199, 19.3284464785803]]
In [10]:
streamCenters = {"mu" : [], "r" : [], "lbr" : [], "radec" : [], 'ra' : [], 'dec' : [], 'theta' : [], 'phi' : [], 'wedge' : []}
wedgeCounter = 80
for i in fitParameters2:
streamCenters["mu"].append(i[3])
streamCenters["mu"].append(i[9])
streamCenters["mu"].append(i[15])
streamCenters["r"].append(i[4])
streamCenters["r"].append(i[10])
streamCenters["r"].append(i[16])
streamCenters["lbr"].append(ac.GC2lbr(i[3], 0.0, i[4], wedgeCounter))
streamCenters["lbr"].append(ac.GC2lbr(i[9], 0.0, i[10], wedgeCounter))
streamCenters["lbr"].append(ac.GC2lbr(i[15], 0.0, i[16], wedgeCounter))
streamCenters["radec"].append(ac.GCToEq(i[3], 0.0, wedgeCounter))
streamCenters["radec"].append(ac.GCToEq(i[9], 0.0, wedgeCounter))
streamCenters["radec"].append(ac.GCToEq(i[15], 0.0, wedgeCounter))
wedgeCounter+=1
streamCenters["lbr"] = np.array(streamCenters["lbr"], dtype=[("l",float),("b",float),("r",float)])
streamCenters["radec"] = np.array(streamCenters["radec"], dtype=[("ra",float),("dec",float)])
c = SkyCoord(stars[0], stars[1], frame='galactic', unit='deg').transform_to('icrs')
correctedRA = []
for i in c.ra:
# print(i * u.deg)
if i > 100 * u.deg:
correctedRA.append((i - 360*u.deg).deg)
else:
correctedRA.append(i.deg)
for i in range(len(streamCenters["radec"]["ra"])):
if streamCenters["radec"]["ra"][i] > 100:
streamCenters["radec"]["ra"][i] -= 360
plt.figure(1,figsize=(15,5))
plt.plot(correctedRA, c.dec.deg, "bo", ms=.05, alpha=.5)
#plt.plot(clusters[0], clusters[1], "go")
plt.plot(streamCenters["radec"]['ra'], streamCenters["radec"]['dec'], 'ro', alpha=.5)
plt.xlabel("ra", fontsize=20)
plt.ylabel("dec", fontsize=20)
plt.xlim([45,-25])
name = list(map(str, range(1,len(streamCenters['r'])+1)))
for i in range(len(name)):
plt.text(streamCenters["radec"]['ra'][i]+1, streamCenters["radec"]['dec'][i], name[i], fontsize = 15)
makeArrows(data)
plt.show()
In [9]:
plt.figure(1,figsize=(15,5))
plt.plot(correctedRA, stars[2], "bo", ms=.05, alpha=.5)
#plt.plot(clusters[0], clusters[2], "go")
plt.plot(streamCenters["radec"]['ra'], streamCenters["lbr"]['r'], 'ro', alpha=.5)
plt.xlabel("ra", fontsize=20)
plt.ylabel("r", fontsize=20)
plt.xlim([45,-25])
for i in range(len(name)):
plt.text(streamCenters["radec"]['ra'][i]+1, streamCenters["lbr"]['r'][i], name[i], fontsize = 15)
plt.show()
In [ ]: