In [1]:
%pylab inline
In [2]:
avs = np.array([0,1,2,3,5,10,20])
x = []
for av in avs:
x.append(np.load("../data/GDR3/gaiadr3_%dext.npy" %(av)))
print(len(x[0]),x[0].dtype)
In [3]:
for name in x[0].dtype.names:
print(name,x[0][name],len(np.unique(x[0][name])))
In [4]:
def return_index_feh(feh):
dfeh = 0.05
offset = 1.5
return(int((feh+offset)/dfeh))
min_value = min(x[0]['meh_ini'])
max_value = max(x[0]['meh_ini'])
min_index_feh = return_index_feh(min_value)
max_index_feh = return_index_feh(max_value)
stretch = max_index_feh - min_index_feh
print('teff values in parsec')
print('teff min max value: ', min_value, max_value )
print('teff min max index: ', min_index_feh, max_index_feh)
print('dimension cut into %d pieces' %(stretch))
def return_index_teff(teff):
dteff = 0.01
return(int(teff/dteff))
min_value = min(x[0]['log_teff'])
max_value = max(x[0]['log_teff'])
min_index_teff = return_index_teff(min_value)
max_index_teff = return_index_teff(max_value)
stretch = max_index_teff - min_index_teff
print('teff values in parsec')
print('teff min max value: ', min_value, max_value )
print('teff min max index: ', min_index_teff, max_index_teff)
print('dimension cut into %d pieces' %(stretch))
def return_index_lum(lum):
dlum = 0.02
offset = 5.
return(int((lum+offset)/dlum))
min_value = min(x[0]['log_lum'])
max_value = max(x[0]['log_lum'])
min_index_lum = return_index_lum(min_value)
max_index_lum = return_index_lum(max_value)
stretch = max_index_lum - min_index_lum
print('lum values in parsec')
print('lum min max value: ', min_value, max_value )
print('lum min max index: ', min_index_lum, max_index_lum )
print('dimension cut into %d pieces' %(stretch))
def return_index(feh,teff,lum):
"""
feh given in dex
teff given in log teff
lum given in log lum
"""
index_feh = return_index_feh(feh)
if index_feh < min_index_feh:
index_feh = min_index_feh
elif index_feh > max_index_feh:
index_feh = max_index_feh
index_teff = return_index_teff(teff)
if index_teff < min_index_teff:
index_teff = min_index_teff
elif index_teff > max_index_teff:
index_teff = max_index_teff
index_teff *= 1000
index_lum = return_index_lum(lum)
if index_lum > max_index_lum:
index_lum = max_index_lum
elif index_lum < min_index_lum:
index_lum = min_index_lum
index_lum *= 1000 * 1000
assert(index_feh >= 0)
assert(index_teff >= 0)
assert(index_lum >= 0)
return (index_feh + index_teff + index_lum)
In [5]:
print('indexing isochrones')
indizes = np.zeros(shape=(len(x[0])),dtype = np.int32)
for i in range(len(x[0])):
indizes[i] = return_index(x[0]["meh_ini"][i],x[0]["log_teff"][i],x[0]["log_lum"][i])
In [6]:
from astropy.io import fits
y = fits.getdata("../../../Galaxia_wrap-master/output/mock_cat_old_gdr3_isochrones/GDR2mock_20.7Gmag.fits")
print(len(y))
print('indexing galaxia')
indizes_g = np.zeros(len(y),dtype = np.int32)
for i in range(len(y)):
indizes_g[i] = return_index(y['feh'][i],np.log10(y['teff'][i]),y['lum'][i])
print('how many unique indizes in galaxia: ',len(np.unique(indizes_g)))
print('independent volume elements in parsec isochrones: ',len(np.unique(indizes)))
print('Are all galaxia indizes covered by the isochrones?')
In [10]:
problems = 0
total_fail=0
distances = []
unique = np.unique(indizes)
for i,item in enumerate(np.unique(indizes_g)):
if i % 1000 == 0:
print(i, len(np.unique(indizes_g)))
try:
assert(item in unique)
except:
problems += 1
temp = y[np.where(indizes_g == item)]
#We are staying with the same metallicity
temp_cut = np.where(indizes%1000==item%1000)
temp_x = x[0][temp_cut]
temp_indizes = indizes[temp_cut]
total_fail += len(temp)
distance = (np.median(np.log10(temp['teff']))-temp_x['log_teff'])**2 + (np.median(temp['lum']) - temp_x['log_lum'])**2
min_distance = min(distance)
print(i,np.sqrt(min_distance))
distances.append(np.sqrt(min_distance))
indizes_g[np.where(indizes_g == item)] = temp_indizes[np.where(distance == min_distance)][0]
print("Problems with %d cells in galaxia with a total of %d stars" %(problems,total_fail))