In [65]:

    

%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/rrc_ID_RA_DEC_I.txt', guess=False, Reader=ascii.Tab)
tbl


    

    
Out[65]:




ID
RA
DEC
I
OGLE-LMC-RRLYR-00015
04:29:43.11
-70:41:11.30
18.833
OGLE-LMC-RRLYR-00020
04:30:00.36
-69:59:09.90
18.9
OGLE-LMC-RRLYR-00029
04:30:43.58
-69:06:48.40
18.56
OGLE-LMC-RRLYR-00032
04:30:52.62
-69:36:38.40
18.764
OGLE-LMC-RRLYR-00035
04:31:07.16
-69:41:07.30
18.91
OGLE-LMC-RRLYR-00042
04:31:34.12
-70:38:49.00
18.805
OGLE-LMC-RRLYR-00046
04:31:45.25
-70:22:39.20
18.789
OGLE-LMC-RRLYR-00054
04:32:01.74
-69:09:49.50
18.998
OGLE-LMC-RRLYR-00078
04:33:04.56
-69:50:01.20
19.142
OGLE-LMC-RRLYR-00085
04:33:11.67
-69:44:39.60
18.905
OGLE-LMC-RRLYR-00089
04:33:19.24
-70:17:29.70
19.922
...
...
...
...
OGLE-LMC-RRLYR-24838
06:13:47.24
-70:50:21.30
18.562
OGLE-LMC-RRLYR-24839
06:13:48.02
-71:06:46.90
18.944
OGLE-LMC-RRLYR-24841
06:13:50.10
-70:11:32.40
18.773
OGLE-LMC-RRLYR-24842
06:13:56.37
-69:01:09.70
18.814
OGLE-LMC-RRLYR-24846
06:14:06.28
-69:09:10.20
18.617
OGLE-LMC-RRLYR-24852
06:14:18.77
-69:53:33.40
18.71
OGLE-LMC-RRLYR-24866
06:14:35.30
-70:26:31.10
18.828
OGLE-LMC-RRLYR-24873
06:14:51.39
-70:36:05.80
19.044
OGLE-LMC-RRLYR-24877
06:15:03.79
-70:11:12.30
18.778
OGLE-LMC-RRLYR-24892
06:16:42.14
-70:46:55.10
18.668
OGLE-LMC-RRLYR-24904
06:18:56.51
-70:50:05.20
18.601

In [66]:

    

#convertendo RA e DEC
ra = coord.Angle(tbl['RA'], unit = u.hour)
dec = coord.Angle(tbl['DEC'], unit = u.degree)


    

In [67]:

    

#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)


    

    




[[  0.   -1. ]
 [ -8.3  -0.3]
 [-12.7   0.6]
 ..., 
 [  6.4  -0.5]
 [ -1.3  -1.1]
 [ -2.1  -1.2]]

In [68]:

    

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[68]:





array([ 0.   ,  0.   ,  0.   , ...,  0.   ,  0.   ,  0.106])

In [69]:

    

magI_correct = tbl['I'] - 1.10*ext
#tbl['I']
len(magI_correct)


    

    
Out[69]:





4958

In [70]:

    

#load CE periods
#ce_p = np.array(np.loadtxt('/home/glauffer/Dropbox/FURG/final_project/data/rrlyr_C_CEperiod.txt',dtype = float, usecols = 1))
ce_p = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/rrlyr_C_CEperiod.txt', guess=False, Reader=ascii.Tab)
len(ce_p['CE_p'])


    

    
Out[70]:





4958

In [71]:

    

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,magI_correct)
print('m - A= ', x[0], '*log P + ', x[1])


    

    




m - A=  -0.0222262207806 *log P +  18.8409494682

In [72]:

    

plt.plot(np.log(ce_p['CE_p']), magI_correct, 'ro', 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()
#plt.xlim(-2.0,0)


    

    
Out[72]:





<matplotlib.legend.Legend at 0x7f946d7d7cf8>






    

In [73]:

    

#calculando o modulo de distancia
u = magI_correct - x[0]*np.log(ce_p['CE_p'])-x[1]
u


    

    
Out[73]:





<Column name='I' unit=None format=None description=None>
array([-0.02111622,  0.03007214, -0.29935391, ..., -0.06607797,
       -0.17953104, -0.36972828])

In [74]:

    

#passando para distancia -> r = 10**(0.2*(u+5))
r = np.zeros(len(u))
r = 10**(0.2*u+0.2*5)#+0.2*1.10*ext)
#sendo a distancia media 50 kpc -> r_lmc = r+50
r_lmc = np.zeros(len(r))
r_lmc = r + 50e3
r_lmc


    

    
Out[74]:





<Column name='I' unit=None format=None description=None>
array([ 50009.9032275 ,  50010.13945069,  50008.71222771, ...,
        50009.70028301,  50009.20648376,  50008.43440292])

In [75]:

    

#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))


    

    




80.4891238235 -69.6773945475






    

In [76]:

    

#Calculando as coordenadas x, y, z

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 [78]:

    

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.plot(x, y, z, 'ro' ,ms=.5)


    

    
Out[78]:





[<mpl_toolkits.mplot3d.art3d.Line3D at 0x7f946f0c00f0>]






    

In [77]: