In [60]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import astropy.coordinates as coord
import astropy.units as u
from mpl_toolkits.mplot3d import Axes3D
from astropy.io import ascii
tbl = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/cepFU_ID_RA_DEC_I.txt', guess=False, Reader=ascii.Tab)
tbl
Out[60]:
In [61]:
#convertendo RA e DEC
ra = coord.Angle(tbl['RA'], unit = u.hour)
dec = coord.Angle(tbl['DEC'], unit = u.degree)
#Converter RA e Dec em X e Y de acordo com o paper do Pejcha (2009)
import numpy as np
x = (ra.degree - 80.35) * np.cos(dec.degree)
y = dec.degree + 69.68
ra_dec = np.around(np.column_stack((x,y)), decimals=1) #arrendondamento para o primeiro decimal apos a virgula
print(ra_dec)
In [62]:
ext_map = np.array(np.loadtxt('/home/glauffer/Dropbox/FURG/final_project/data/extinction_map/lmcextmap.dat', skiprows=29))
x_map = np.column_stack((ext_map.T[0], ext_map.T[1]))
y_map = np.column_stack((ext_map.T[2], ext_map.T[3]))
ext = np.ones(len(ra_dec))
for i in range(len(ra_dec)):
if ra_dec[i,0]>(-4.7): # or ():
ext[i] = 0
if ra_dec[i,0]>4.9:
ext[i] = 0
if ra_dec[i,1]>(-3.4): #or
ext[i] = 0
if ra_dec[i,1]>3.2:
ext[i] = 0
else:
try:
ext[i] = ext_map[np.where((x_map.T[0] == ra_dec[i,0]) & (y_map.T[0] == ra_dec[i,1]))][0][4]
except IndexError:
ext[i] = 0
ext
Out[62]:
In [63]:
magI_correct = tbl['I'] - 1.10*ext
#magI_correct = tbl['I'] - 0.45
#tbl['I']
magI_correct
Out[63]:
In [64]:
ce_p = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/cep_FU_CEperiod.txt', guess=False, Reader=ascii.Tab)
len(ce_p['CE_p'])
Out[64]:
In [65]:
#Removendo as entrelas que nao possuem banda I
exception = np.where(tbl['I'] <= 10)
new_magI_correct = np.delete(magI_correct, exception)
len(new_magI_correct)
Out[65]:
In [66]:
A = np.column_stack((np.log(ce_p['CE_p']), np.ones_like(ce_p['CE_p'])))
x, res, rank, s = np.linalg.lstsq(A,new_magI_correct)
print('m - A= ', x[0], '*log P + ', x[1])
In [67]:
plt.plot(np.log(ce_p['CE_p']), new_magI_correct, 'go', ms=0.5)
plt.plot(np.log(ce_p['CE_p']), x[0]*np.log(ce_p['CE_p']) + x[1], 'r-', label='m - A= %f*log P +%f'%(x[0],x[1]))
plt.gca().invert_yaxis()
plt.legend()
Out[67]:
In [68]:
#calculando o modulo de distancia
u = new_magI_correct - x[0]*np.log(ce_p['CE_p'])-x[1]
#print(new_magI_correct[1], new_magI_correct[1]-u[1])
#passando para distancia -> r = 10**(0.2*(u+5))
r = np.zeros(len(u))
r =10**(0.2*u+0.2*5)
#new_ext = np.delete(ext, exception)
#r = 10**(0.2*(u+5+1.10*new_ext))
#print(r[1])
#sendo a distancia media 50 kpc -> r_lmc = r+50
r_lmc = np.zeros(len(r))
r_lmc = (r + 50e3)/1000
print(r_lmc)
In [69]:
#plot das coordenadas RA e DEC
plt.plot(ra.degree, dec.degree, 'bo', ms=1.)
plt.plot(np.mean(ra.degree), np.mean(dec.degree), 'go')
print(np.mean(ra.degree), np.mean(dec.degree))
In [75]:
#Calculando as coordenadas x, y, z
new_ra = np.delete(ra.radian, exception)
new_dec =np.delete(dec.radian, exception)
r_lmc_0 = 50
ra_0 = np.mean(ra.radian)
dec_0 = np.mean(dec.radian)
x = -r_lmc * np.sin(new_ra - ra_0) * np.cos(new_dec)
y = r_lmc * np.sin(new_dec)*np.cos(dec_0) - r_lmc * np.sin(dec_0) * np.cos(new_ra - ra_0) * np.cos(new_dec)
z = r_lmc_0 - r_lmc * np.sin(new_dec) * np.sin(dec_0) - r_lmc * np.cos(dec_0)*np.cos(new_ra) - ra_0 * np.cos(new_dec)
In [76]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.plot(x, y, z, 'go' ,ms=1.5)
Out[76]:
In [72]:
exception
Out[72]:
In [73]:
fu_exception = np.array([33, 36, 68, 73, 138, 186, 243, 260, 290, 291, 302, 334, 352, 370, 423, 534, 535, 546, 562,
567, 619, 625, 717, 875, 897, 900, 920, 1245, 1393, 1564,1744])
fu_exception
Out[73]:
In [77]:
np.savetxt('fu50e3_plot.txt', np.dstack((x, y, z))[0], fmt='%s')
In [74]: