Extinção Interstelar

é a dispersão, desvio ou avermelhamento (ou todos esses efeitos) causados na luz emitida por uma fonte devido ao meio interstelar.

A extinção (ou reddening) é simbolizado pela letra $A_\lambda$ com um subscrito que identifica a banda da correção.

Existem alguns mapas de extinção, sendo Zaritsky, (2004) um dos mais citados. Porém, vou utilizar um mapa criado a partir dos dados de RRLyr da LMC contidos no OGLE-III (Pejcha, 2009) e disponíveis em http://www.astro.princeton.edu/~pejcha/lmc_extmap/.

O arquivo lmcextmap.dat possui os dados do mapa de extinção na seguinte ordem:

MinX,   MaxX,   MinY,   MaxY,   E(V-I),     E_rms,    Delta E(V-I),     numero de estrelas

em que as coordenadas $x$ e $y$ são obtidas pela conversão, $$ x = (\alpha -80.35)\cdot \cos(\delta) $$ $$ y = \delta +69.68 $$ onde $\alpha$ é ascenção reta e $\delta$ é a declinação, ambas dadas em graus.

Para obter a extinção, devo converter os dados do OGLE para as coordenadas $x$ e $y$ e obter E(V-I) da tabela lmcextmap.dat. A extinção na banda I é obtida pelo calculo, $$A_I = R_{VI} E(V-I)$$ em que $R_{VI} = 1.10$ como utilizado no paper.

Porem, no catálogo OGLE, $\alpha$ é dado em horas. Primeiramente vou converter $\alpha$ em graus.


In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
#Ler a tabela de dados com RA e Dec de cada RRLyr AB
# a estrela OGLE-LMC-RRLYR-15485 (posição 11013 no arquivo) não possui DatabaseNumber, estava dando problema para ler os dados. Coloquei manualmente o valor 0000
from astropy.io import ascii
tbl = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/identification_onlyRRAB.dat', guess=False, Reader=ascii.Tab, header_start = 0, data_start = 1)
tbl


Out[1]:
NameFieldNTypeRADeccol6col7
OGLE-LMC-RRLYR-00001LMC158.8216RRab04:27:45.47-70:43:12.00----
OGLE-LMC-RRLYR-00002LMC158.5211RRab04:27:58.46-70:10:59.40----
OGLE-LMC-RRLYR-00003LMC158.6247RRab04:28:08.50-70:21:22.80----
OGLE-LMC-RRLYR-00004LMC158.6271RRab04:28:19.31-70:22:12.00----
OGLE-LMC-RRLYR-00005LMC158.5262RRab04:28:21.07-70:08:54.40----
OGLE-LMC-RRLYR-00006LMC158.8286RRab04:28:28.19-70:41:44.40----
OGLE-LMC-RRLYR-00007LMC158.7320RRab04:28:38.81-70:28:19.80----
OGLE-LMC-RRLYR-00008LMC158.6346RRab04:28:57.78-70:17:16.60----
OGLE-LMC-RRLYR-00010LMC158.82019RRab04:29:19.47-70:41:24.10----
OGLE-LMC-RRLYR-00011LMC158.81805RRab04:29:22.01-70:39:49.80----
OGLE-LMC-RRLYR-00012LMC157.7245RRab04:29:24.99-69:50:38.50----
........................
OGLE-LMC-RRLYR-24895LMC214.34219RRab06:17:07.43-70:20:05.70----
OGLE-LMC-RRLYR-24896LMC214.34223RRab06:17:08.86-70:25:08.10----
OGLE-LMC-RRLYR-24897LMC215.41147RRab06:17:09.27-70:49:20.20----
OGLE-LMC-RRLYR-24898LMC214.44923RRab06:17:18.69-70:12:32.20----
OGLE-LMC-RRLYR-24899LMC215.3818RRab06:17:22.52-70:53:34.40----
OGLE-LMC-RRLYR-24900LMC215.11168RRab06:17:22.54-71:14:11.80----
OGLE-LMC-RRLYR-24901LMC215.34373RRab06:17:36.32-70:58:51.20----
OGLE-LMC-RRLYR-24902LMC215.34516RRab06:18:23.05-70:57:28.20----
OGLE-LMC-RRLYR-24903LMC215.24621RRab06:18:28.97-71:06:12.50----
OGLE-LMC-RRLYR-24905LMC215.13991RRab06:19:05.65-71:11:09.60----
OGLE-LMC-RRLYR-24906LMC215.13996RRab06:19:06.15-71:11:28.30----

In [2]:
#Converter RA de horas para graus
import astropy.coordinates as coord
import astropy.units as u
ra = coord.Angle(tbl['RA'], unit = u.hour)
ra.degree


Out[2]:
array([ 66.93945833,  66.99358333,  67.03541667, ...,  94.62070833,
        94.77354167,  94.775625  ])

In [3]:
#Converter Dec em graus
dec = coord.Angle(tbl['Dec'], unit = u.degree)
dec.degree


Out[3]:
array([-70.72      , -70.18316667, -70.35633333, ..., -71.10347222,
       -71.186     , -71.19119444])

Agora, tendo os valores de RA e Dec em graus, vou converter para o sistema X e Y utilizado no paper


In [4]:
#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)
#np.savetxt('ra_dec.txt',ra_dec, fmt='%s')


[[ 0.5 -1. ]
 [-6.4 -0.5]
 [-4.3 -0.7]
 ..., 
 [-5.8 -1.4]
 [-6.9 -1.5]
 [-7.  -1.5]]

Ler os dados do extinction map em forma de tabela


In [5]:
#ext_map = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/extinction_map/lmcextmap.dat', 
#                     guess=False, Reader=ascii.Basic)
ext_map = np.array(np.loadtxt('/home/glauffer/Dropbox/FURG/final_project/data/extinction_map/lmcextmap.dat', skiprows=29))
ext_map


Out[5]:
array([[-4.7, -4.5, -3.4, ...,  0. ,  0. ,  0. ],
       [-4.7, -4.5, -3.2, ...,  0. ,  0. ,  0. ],
       [-4.7, -4.5, -3. , ...,  0. ,  0. ,  0. ],
       ..., 
       [ 4.9,  5.1,  2.8, ...,  0. ,  0. ,  0. ],
       [ 4.9,  5.1,  3. , ...,  0. ,  0. ,  0. ],
       [ 4.9,  5.1,  3.2, ...,  0. ,  0. ,  0. ]])

In [6]:
#Criar array 2d com os dados min e max para X e Y
#x_map = np.column_stack((ext_map['minX'], ext_map['maxX']))
#y_map = np.column_stack((ext_map['minY'], ext_map['maxY']))
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[6]:
array([ 0.144,  0.   ,  0.   , ...,  0.   ,  0.   ,  0.   ])

In [7]:
#np.savetxt('teste_ext.txt', ext, fmt='%s')
#Ler os dados de magnitude das RRLYR AB
#Nome MagI MagV OGLE_period RA DEC CE_period...
magI = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/ID_mag_Period_RA_DE.dat')
magI


Out[7]:
col1col2col3col4col5col6col7
OGLE-LMC-RRLYR-0000118.79819.3910.634752204:27:45.47-70:43:12.000.6347
OGLE-LMC-RRLYR-0000218.83219.4940.625525404:27:58.46-70:10:59.400.6255
OGLE-LMC-RRLYR-0000318.65119.30.656497504:28:08.50-70:21:22.800.6565
OGLE-LMC-RRLYR-0000418.76419.3350.512374504:28:19.31-70:22:12.000.5124
OGLE-LMC-RRLYR-0000518.93919.640.643352404:28:21.07-70:08:54.400.6434
OGLE-LMC-RRLYR-0000618.76819.4210.600017904:28:28.19-70:41:44.400.6
OGLE-LMC-RRLYR-0000718.75919.4370.5528504:28:38.81-70:28:19.800.5528
OGLE-LMC-RRLYR-0000818.68919.3880.787739504:28:57.78-70:17:16.600.7877
OGLE-LMC-RRLYR-0001018.83919.4950.594074704:29:19.47-70:41:24.100.5941
OGLE-LMC-RRLYR-0001118.84219.4270.553567904:29:22.01-70:39:49.800.5536
OGLE-LMC-RRLYR-0001218.71219.2930.609436604:29:24.99-69:50:38.500.6094
.....................
OGLE-LMC-RRLYR-2489518.75419.4920.700194106:17:07.43-70:20:05.700.7002
OGLE-LMC-RRLYR-2489618.71319.3640.591790206:17:08.86-70:25:08.100.5918
OGLE-LMC-RRLYR-2489718.99619.6170.53910906:17:09.27-70:49:20.200.5391
OGLE-LMC-RRLYR-2489818.82119.4290.507191306:17:18.69-70:12:32.200.5072
OGLE-LMC-RRLYR-2489918.71119.2350.557475406:17:22.52-70:53:34.400.5575
OGLE-LMC-RRLYR-2490018.81919.360.540443306:17:22.54-71:14:11.800.5404
OGLE-LMC-RRLYR-2490118.56219.1390.605990406:17:36.32-70:58:51.200.606
OGLE-LMC-RRLYR-2490218.73219.2760.55704506:18:23.05-70:57:28.200.557
OGLE-LMC-RRLYR-2490318.81919.4470.629634506:18:28.97-71:06:12.500.6296
OGLE-LMC-RRLYR-2490518.64619.2180.513162706:19:05.65-71:11:09.600.5132
OGLE-LMC-RRLYR-2490618.69519.3130.611652606:19:06.15-71:11:28.300.6116

In [8]:
len(magI['col2'])


Out[8]:
17693

In [9]:
#corrigir a extincao 
magI_correct = magI['col2'] - 1.10*ext

In [10]:
len(magI_correct)


Out[10]:
17693

In [11]:
#relacao pl -> magI - ext = a*log(P) + B + u
#em que u eh o modulo da distancia, ou seja, quanto se distancia do eixo x

#plot da dispersao
plt.plot(np.log(magI['col7']), magI_correct, 'bo', ms=0.5)
plt.gca().invert_yaxis()



In [12]:
#Least squares regression for PL relation
A = np.column_stack((np.log(magI['col7']), np.ones_like(magI['col7'])))
x, res, rank, s = np.linalg.lstsq(A,magI_correct)
print('m - A= ', x[0], '*log P + ', x[1])


m - A=  -0.720619071784 *log P +  18.3639869006

In [13]:
plt.plot(np.log(magI['col7']), magI_correct, 'bo', ms=0.4)
plt.plot(np.log(magI['col7']), x[0]*np.log(magI['col7']) + x[1], 'r-', label='m - A= %f*log P +%f'%(x[0],x[1]))
plt.gca().invert_yaxis()
plt.legend()
plt.xlim(-2.0,0)


Out[13]:
(-2.0, 0)

In [14]:
#u = magI_correct - x[0]*np.log(P) - x[1]
u = np.zeros(len(magI_correct))
u = magI_correct - x[0]*np.log(magI['col7'])-x[1]

In [15]:
u


Out[15]:
<Column name='col2' unit=None format=None description=None>
array([-0.05198237,  0.12989579, -0.01624689, ...,  0.12160385,
       -0.19870442, -0.02329858])

In [16]:
#passando para distancia -> r = 10**(0.2*(u+5))
r = np.zeros(len(u))
r = 10**(0.2*u+0.2*5)

In [17]:
r


Out[17]:
<Column name='col2' unit=None format=None description=None>
array([  9.76345493,  10.61644605,   9.92545952, ...,  10.57598365,
         9.1255514 ,   9.8932796 ])

In [18]:
#sendo a distancia media 50 kpc -> r_lmc = r+50
r_lmc = np.zeros(len(r))
r_lmc = r + 50e3

In [19]:
#r_lmc
#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.hour), np.mean(dec.degree))
#plt.plot(ra.radian, dec.radian, 'bo', ms=1.)
#plt.plot(np.mean(ra.radian), np.mean(dec.radian), 'go')


5.35708330947 -69.6454403766

In [20]:
#theta = dec.radian
theta = ra.radian
theta


Out[20]:
array([ 1.16831395,  1.16925861,  1.16998874, ...,  1.6514429 ,
        1.65411035,  1.65414671])

In [21]:
ax = plt.subplot(111, polar=True)
#l=17000
#plt.gca(projection='polar')
#plt.plot(2*np.pi*theta[0:l], r_lmc[0:l], 'ro')#dec.degree, 'ro')
#plt.plot(np.sin(theta), r_lmc, 'ro', ms=0.5)
plt.plot(theta, r_lmc, 'bo', ms=1.0)
ax.set_rlim(0, 55e3)
#ax.set_ylim(30,70)
#ax.set_yticks([35, 50, 60, 70]) #['35 kpc', '50 kpc', '60 kpc', '70 kpc']
#plt.thetagrids([20,40,60], labels=None, fmt='%d', frac = 1.1)
#plt.scatter(theta, r, 'ro')
#plt.show()
#plt.gca(projection='polar')


Out[21]:
(0, 55000.0)

In [22]:
#Calculando as coordenadas x, y, z
# x = -r_lmc * np.sin(ra - ra_0) * np.cos(dec)
# y = r_lmc * np.sin(dec)*np.cos(dec_0) − r_lmc * np.sin(dec_0) * np.cos(ra - ra_0) * np.cos(dec)
# z = r_lmc_0 - r_lmc * np.sin(dec) * np.sin(dec_0) - r_lmc * np.cos(dec_0)*np.cos(ra) - \ra_0 * np.cos(dec)

x, y, z = np.zeros(len(r_lmc)), np.zeros(len(r_lmc)), np.zeros(len(r_lmc))
r_lmc_0 = 50e3
ra_0 = np.mean(ra.radian)
dec_0 = np.mean(dec.radian)

x = -r_lmc * np.sin(ra.radian - ra_0) * np.cos(dec.radian)
y = r_lmc * np.sin(dec.radian)*np.cos(dec_0) - r_lmc * np.sin(dec_0) * np.cos(ra.radian - ra_0) * np.cos(dec.radian)
z = r_lmc_0 - r_lmc * np.sin(dec.radian) * np.sin(dec_0) - r_lmc * np.cos(dec_0)*np.cos(ra.radian) - ra_0 * np.cos(dec.radian)

In [24]:
from mpl_toolkits.mplot3d import Axes3D
plt.clf()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.plot(x, y, z, 'bo' ,ms=0.2)
#fig.show()


Out[24]:
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x7f0ac02ef470>]
<matplotlib.figure.Figure at 0x7f0ab8bd0240>

In [23]: