-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathhighsnobjects.py
126 lines (83 loc) · 2.66 KB
/
highsnobjects.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
import math
import astropy
import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as fits
from mayavi.mlab import *
#read in data from Gaia catalog
file=fits.open('tgas-source.fits')
data_list=file[1]
#read in data from 2mass
file2=fits.open('tgas-matched-2mass.fits')
data_list_2=file2[1]
#Pull data we want from .fits file
#Parallax data
parallax=data_list.data['parallax'] #in mas
parallax_error=data_list.data['parallax_error'] #error
#proper motion
propermotion_ra=data_list.data['pmra']
propermotion_dec=data_list.data['pmdec']
#RA and DEC
ra=data_list.data['ra']
dec=data_list.data['dec']
#magnitudes
#g_band=data_list.data['phot_g_mean_mag']
j_band=data_list_2.data['j_mag']
k_band=data_list_2.data['k_mag']
#h_band=data_list_2.data['h_mag']
#calculate SN ratio
ratio=parallax/parallax_error
#select high SN data that we want
highSNindices = ratio > 128.
#locations of data we want
#np.where(highSNindices)
#distances we want from high SN data
distance=1./parallax[highSNindices] #in Kpc
#Calculate transverse velocity in RA and DEC then in real space
transv_ra=4.74*propermotion_ra[highSNindices]*distance #km/s
transv_dec=4.74*propermotion_dec[highSNindices]*distance #km/s
transverse_vsqared=transv_ra**2+transv_dec**2
transverse_v=transverse_vsqared**(.5)
#calcilate color G-J
j_k=j_band[highSNindices]-k_band[highSNindices]
equalzero= (j_k!=0.) & (k_band[highSNindices]!=0.)
j_select=j_band[highSNindices][equalzero]
k_select=k_band[highSNindices][equalzero]
j_kprime=j_select-k_select
#calculate absolute magnitude
distance_pc=distance[equalzero]*10**3.
absolute_mag=j_select-(5.*(np.log10(distance_pc)-1))
ra_selected=ra[highSNindices][equalzero]
dec_selected=dec[highSNindices][equalzero]
#conversion
x_cord=distance_pc*np.sin(dec_selected)*np.cos(ra_selected)
y_cord=distance_pc*np.sin(dec_selected)*np.sin(ra_selected)
z_cord=distance_pc*np.cos(dec_selected)
v_z=np.zeros_like(z_cord)
quiver3d(x_cord, y_cord, z_cord, transv_ra[equalzero], transv_dec[equalzero],v_z)
show()
#plotting
#H-R diagram
#make figure
rect1=0.1,0.1,0.75,0.75
fig1=plt.figure(1)
ax1=fig1.add_axes(rect1)
ax1.plot(j_kprime,absolute_mag, "o", lw=0, markersize=1)
ax1.invert_yaxis()
ax1.set_ylabel("Absolute Magnitude, J")
ax1.set_xlabel("J-K")
plt.show()
#import pdb; pdb.set_trace()
#distance distribution
#plt.hist(distance,bins=100,log=True)
#plt.show()
#Relative velocity distribution
#plt.hist(transverse_v,bins=100,log=True)
#plt.show()
#plt.savefig('highsnVelocities.png')
#3D Map of selected stars
#fig=plt.figure()
#ax=fig.add_subplot(111,projection='3d')
#ax.scatter(ra[highSNindices],dec[highSNindices],distance)
#plt.show()
#plt.savefig('highsn3Dmap.png')