notebook.community

Edit and run

Leitura dos dados





In [55]:



    


%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
plt.rcParams['savefig.dpi'] = 200

cep_fu = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/cepFU_ID_RA_DEC_I.txt', guess=False, Reader=ascii.Tab)
cep_fo = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/cepFO_ID_RA_DEC_I.txt', guess=False, Reader=ascii.Tab)
rrc = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/rrc_ID_RA_DEC_I.txt', guess=False, Reader=ascii.Tab)
rrab = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/identification_onlyRRAB.dat', guess=False, Reader=ascii.Tab, header_start = 0, data_start = 1)



    











In [56]:



    


#convertendo RA e DEC
ra_cepfu_1 = coord.Angle(cep_fu['RA'], unit = u.hour)
dec_cepfu_1 = coord.Angle(cep_fu['DEC'], unit = u.degree)

ra_cepfo = coord.Angle(cep_fo['RA'], unit = u.hour)
dec_cepfo = coord.Angle(cep_fo['DEC'], unit = u.degree)

ra_rrc = coord.Angle(rrc['RA'], unit = u.hour)
dec_rrc = coord.Angle(rrc['DEC'], unit = u.degree)

ra_rrab = coord.Angle(rrab['RA'], unit = u.hour)
dec_rrab = coord.Angle(rrab['Dec'], unit = u.degree)



    











In [57]:



    


#leitura dos periodos
cepfu_p = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/cep_FU_CEperiod.txt', guess=False, Reader=ascii.Tab)
cepfo_p = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/cep_FO_CEperiod.txt', guess=False, Reader=ascii.Tab)
rrc_p = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/rrlyr_C_CEperiod.txt', guess=False, Reader=ascii.Tab)
rrab_p = ascii.read('/home/glauffer/Dropbox/FURG/final_project/data/ID_mag_Period_RA_DE.dat')



    











In [58]:



    


#remover estrelas sem banda I das CEP FU
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])

mag_fu = np.delete(cep_fu['I'], fu_exception)

#ra_cepfu = np.delete(ra_cepfu_1, fu_exception)
#dec_cepfu = np.delete(dec_cepfu_1, fu_exception)



    











In [59]:



    


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

Removendo a extincao e Relacao PL

CEP FU





In [60]:



    


#Converter RA e Dec em X e Y de acordo com o paper do Pejcha (2009)
import numpy as np
x_fu = (ra_cepfu_1.degree - 80.35) * np.cos(dec_cepfu_1.degree)
y_fu = dec_cepfu_1.degree + 69.68
ra_dec_cepfu = np.around(np.column_stack((x_fu,y_fu)), decimals=1) #arrendondamento para o primeiro decimal apos a virgula
#print(ra_dec_cepfu)



    











In [61]:



    


ext_fu = np.ones(len(ra_dec_cepfu))

for i in range(len(ra_dec_cepfu)):
    if ra_dec_cepfu[i,0]>(-4.7): # or ():
        ext_fu[i] = 0
    if ra_dec_cepfu[i,0]>4.9:
        ext_fu[i] = 0
    if ra_dec_cepfu[i,1]>(-3.4): #or 
        ext_fu[i] = 0
    if ra_dec_cepfu[i,1]>3.2:
        ext_fu[i] = 0
    else:
        try:
            ext_fu[i] = ext_map[np.where((x_map.T[0] == ra_dec_cepfu[i,0]) & (y_map.T[0] == ra_dec_cepfu[i,1]))][0][4]
        except IndexError:
            ext_fu[i] = 0
ext_fu1 = np.delete(ext_fu, fu_exception)
magI_fu = mag_fu - 1.10*ext_fu1



    











In [62]:



    


A_fu = np.column_stack((np.log(cepfu_p['CE_p']), np.ones_like(cepfu_p['CE_p'])))
x_fu, res_fu, rank_fu, s_fu = np.linalg.lstsq(A_fu, magI_fu)
print('m_FU = ', x_fu[0], '*log P + ', x_fu[1])



    


















    






m_FU =  -1.29201094747 *log P +  16.8783009955

CEP FO





In [63]:



    


#Converter RA e Dec em X e Y de acordo com o paper do Pejcha (2009)
import numpy as np
x_fo = (ra_cepfo.degree - 80.35) * np.cos(dec_cepfo.degree)
y_fo = dec_cepfo.degree + 69.68
ra_dec_cepfo = np.around(np.column_stack((x_fo,y_fo)), decimals=1) #arrendondamento para o primeiro decimal apos a virgula
#print(ra_dec_cepfo)

ext_fo = np.ones(len(ra_dec_cepfo))

for i in range(len(ra_dec_cepfo)):
    if ra_dec_cepfo[i,0]>(-4.7): # or ():
        ext_fo[i] = 0
    if ra_dec_cepfo[i,0]>4.9:
        ext_fo[i] = 0
    if ra_dec_cepfo[i,1]>(-3.4): #or 
        ext_fo[i] = 0
    if ra_dec_cepfo[i,1]>3.2:
        ext_fo[i] = 0
    else:
        try:
            ext_fo[i] = ext_map[np.where((x_map.T[0] == ra_dec_cepfo[i,0]) & (y_map.T[0] == ra_dec_cepfo[i,1]))][0][4]
        except IndexError:
            ext_fo[i] = 0

magI_fo = cep_fo['I']- 1.10*ext_fo



    











In [64]:



    


A_fo = np.column_stack((np.log(cepfo_p['CE_p']), np.ones_like(cepfo_p['CE_p'])))
x_fo, res_fo, rank_fo, s_fo = np.linalg.lstsq(A_fo, magI_fo)
print('m_FO = ', x_fo[0], '*log P + ', x_fo[1])



    


















    






m_FO =  -1.57341725022 *log P +  16.5583045317

Relacao PL das CEP





In [65]:



    


plt.plot(np.log(cepfo_p['CE_p']), magI_fo, 'y>', ms=2.5, label='Cep FO')
plt.plot(np.log(cepfu_p['CE_p']), magI_fu, 'bo', ms=2.5, label='Cep FU')
plt.plot(np.log(cepfu_p['CE_p']), x_fu[0]*np.log(cepfu_p['CE_p']) + x_fu[1], 'r-')
plt.plot(np.log(cepfo_p['CE_p']), x_fo[0]*np.log(cepfo_p['CE_p']) + x_fo[1], 'r-')
plt.gca().invert_yaxis()
#plt.legend()
plt.legend(loc=2)
plt.title('Relação PL - Cefeidas')
plt.xlabel(r"$ \log P$")
plt.ylabel(r"$\bar{m}_I$")
plt.savefig('cep_relacaoPL.png', bbox_inches='tight')
print('m_FU = ', x_fu[0], '*log P + ', x_fu[1])
print('m_FO = ', x_fo[0], '*log P + ', x_fo[1])



    


















    






m_FU =  -1.29201094747 *log P +  16.8783009955
m_FO =  -1.57341725022 *log P +  16.5583045317

RRLYR AB





In [66]:



    


#Converter RA e Dec em X e Y de acordo com o paper do Pejcha (2009)
import numpy as np
x_ab = (ra_rrab.degree - 80.35) * np.cos(dec_rrab.degree)
y_ab = dec_rrab.degree + 69.68
ra_dec_rrab = np.around(np.column_stack((x_ab,y_ab)), decimals=1) #arrendondamento para o primeiro decimal apos a virgula
#print(ra_dec_rrab)

ext_ab = np.ones(len(ra_dec_rrab))

for i in range(len(ra_dec_rrab)):
    if ra_dec_rrab[i,0]>(-4.7): # or ():
        ext_ab[i] = 0
    if ra_dec_rrab[i,0]>4.9:
        ext_ab[i] = 0
    if ra_dec_rrab[i,1]>(-3.4): #or 
        ext_ab[i] = 0
    if ra_dec_rrab[i,1]>3.2:
        ext_ab[i] = 0
    else:
        try:
            ext_ab[i] = ext_map[np.where((x_map.T[0] == ra_dec_rrab[i,0]) & (y_map.T[0] == ra_dec_rrab[i,1]))][0][4]
        except IndexError:
            ext_ab[i] = 0

magI_ab = rrab_p['col2']- 1.10*ext_ab



    











In [67]:



    


A_ab = np.column_stack((np.log(rrab_p['col7']), np.ones_like(rrab_p['col7'])))
x_ab, res_ab, rank_ab, s_ab = np.linalg.lstsq(A_ab, magI_ab)
print('m_ab = ', x_ab[0], '*log P + ', x_ab[1])



    


















    






m_ab =  -0.720619071784 *log P +  18.3639869006

RRLYR C





In [68]:



    


#Converter RA e Dec em X e Y de acordo com o paper do Pejcha (2009)
import numpy as np
x_c = (ra_rrc.degree - 80.35) * np.cos(dec_rrc.degree)
y_c = dec_rrc.degree + 69.68
ra_dec_rrc = np.around(np.column_stack((x_c,y_c)), decimals=1) #arrendondamento para o primeiro decimal apos a virgula
#print(ra_dec_rrc)

ext_c = np.ones(len(ra_dec_rrc))

for i in range(len(ra_dec_rrc)):
    if ra_dec_rrc[i,0]>(-4.7): # or ():
        ext_c[i] = 0
    if ra_dec_rrc[i,0]>4.9:
        ext_c[i] = 0
    if ra_dec_rrc[i,1]>(-3.4): #or 
        ext_c[i] = 0
    if ra_dec_rrc[i,1]>3.2:
        ext_c[i] = 0
    else:
        try:
            ext_c[i] = ext_map[np.where((x_map.T[0] == ra_dec_rrc[i,0]) & (y_map.T[0] == ra_dec_rrc[i,1]))][0][4]
        except IndexError:
            ext_c[i] = 0
#ext_c
magI_c = rrc['I']- 1.10*ext_c



    











In [69]:



    


A_c = np.column_stack((np.log(rrc_p['CE_p']), np.ones_like(rrc_p['CE_p'])))
x_c, res_c, rank_c, s_c = np.linalg.lstsq(A_c, magI_c)
print('m_c = ', x_c[0], '*log P + ', x_c[1])



    


















    






m_c =  -0.0222262207806 *log P +  18.8409494682

relacao PL das RRLYR





In [70]:



    


plt.plot(np.log(rrc_p['CE_p']), magI_c, 'g<', ms=2., label='RRLyr C')
plt.plot(np.log(rrab_p['col7']), magI_ab, 'ko', ms=2., label='RRLyr AB')
plt.plot(np.log(rrc_p['CE_p']), x_c[0]*np.log(rrc_p['CE_p']) + x_c[1], 'r-')
plt.plot(np.log(rrab_p['col7']), x_ab[0]*np.log(rrab_p['col7']) + x_ab[1], 'r-')
plt.gca().invert_yaxis()
#plt.legend()
plt.legend(loc=2)
print('m_ab = ', x_ab[0], '*log P + ', x_ab[1])
print('m_c = ', x_c[0], '*log P + ', x_c[1])



    


















    






m_ab =  -0.720619071784 *log P +  18.3639869006
m_c =  -0.0222262207806 *log P +  18.8409494682










    
























In [71]:



    


plt.plot(np.log(rrc_p['CE_p']), magI_c, 'g<', ms=2., label='RRLyr C')
plt.plot(np.log(rrc_p['CE_p']), x_c[0]*np.log(rrc_p['CE_p']) + x_c[1], 'r-')
plt.gca().invert_yaxis()
#plt.legend()
plt.xlim(-2.0,0)
plt.legend(loc=2)
plt.title('Relação PL - RRLYR C')
plt.xlabel(r"$ \log P$")
plt.ylabel(r"$\bar{m}_I$")
plt.savefig('rrc_relacaoPL.png', bbox_inches='tight')



    


















    
























In [72]:



    


plt.plot(np.log(rrab_p['col7']), magI_ab, 'ko', ms=.8, label='RRLyr AB')
plt.plot(np.log(rrab_p['col7']), x_ab[0]*np.log(rrab_p['col7']) + x_ab[1], 'r-')
plt.gca().invert_yaxis()
#plt.legend()
plt.xlim(-2.0,0)
plt.legend(loc=2)
plt.title('Relação PL - RRLYR AB')
plt.xlabel(r"$ \log P$")
plt.ylabel(r"$\bar{m}_I$")
plt.savefig('rrab_relacaoPL.png', bbox_inches='tight')

Calculando o modulo de distancia

CEP FU





In [73]:



    


#calculando o modulo de distancia
u_fu = magI_fu - x_fu[0]*np.log(cepfu_p['CE_p'])-x_fu[1]
#passando para distancia em kpc mais a distancia media-> r = 10**(0.2*(u+5))

r_fu = (10**(0.2*u_fu+0.2*5) + 50e3) / 1000

CEP FO





In [74]:



    


#calculando o modulo de distancia
u_fo = magI_fo - x_fo[0]*np.log(cepfo_p['CE_p'])-x_fo[1]
#passando para distancia -> r = 10**(0.2*(u+5))
#r = np.zeros(len(u))
r_fo = (10**(0.2*u_fo+0.2*5) + 50e3) / 1000

RR AB





In [75]:



    


#calculando o modulo de distancia
u_ab = magI_ab - x_ab[0]*np.log(rrab_p['col7'])-x_ab[1]
r_ab = (10**(0.2*u_ab+0.2*5) + 50e3) / 1000

RR C





In [76]:



    


#calculando o modulo de distancia
u_c = magI_c - x_c[0]*np.log(rrc_p['CE_p'])-x_c[1]
r_c = (10**(0.2*u_c+0.2*5) + 50e3) / 1000

Determinando a centroide





In [144]:



    


ra_m = (np.mean(ra_rrab.degree) + np.mean(ra_rrc.degree) + np.mean(ra_cepfu_1.degree) + np.mean(ra_cepfo.degree)) / 4
dec_m = (np.mean(dec_rrab.degree) + np.mean(dec_rrc.degree) + np.mean(dec_cepfu_1.degree) + np.mean(dec_cepfo.degree)) / 4
print(ra_m, dec_m)
#plt.plot(ra_rrab.degree, dec_rrab.degree, 'ko', ms=2., label='RRLyr AB')
#plt.plot(ra_rrc.degree, dec_rrc.degree, 'g<', ms=2., label='RRLyr C')
#plt.plot(ra_cepfu_1.degree, dec_cepfu_1.degree, 'y>', ms=2., label='Cep FU')
#plt.plot(ra_cepfo.degree, dec_cepfo.degree, 'bo', ms=2., label='Cep FO')
plt.plot(ra_rrab.degree, dec_rrab.degree, 'k.', ms=2., label='RRLyr AB')
plt.plot(ra_rrc.degree, dec_rrc.degree, 'g.', ms=2., label='RRLyr C')
plt.plot(ra_cepfu_1.degree, dec_cepfu_1.degree, 'y.', ms=2., label='Cep FU')
plt.plot(ra_cepfo.degree, dec_cepfo.degree, 'b.', ms=2., label='Cep FO')

plt.plot(ra_m, dec_m, 'ro', ms=5.)
plt.title('Coordenadas Equatoriais')
plt.ylabel('Declinação $(\delta)$ [graus]')
plt.xlabel('Ascensão Reta '+r"$(\alpha)$"+' [graus]')
plt.legend(prop={'size':6})
plt.savefig('equatorial.png', bbox_inches='tight')
#plt.plot(ra_rrab.degree, dec_rrab.degree, 'ko')



    


















    






80.2395185619 -69.6084739941

Plot 3d





In [78]:



    


r_lmc_0 = (50e3) / 1000
#r_lmc_0 = 50
ra_0 = (np.mean(ra_rrab.radian) + np.mean(ra_rrc.radian) + np.mean(ra_cepfu_1.radian) + np.mean(ra_cepfo.radian)) / 4
dec_0 = (np.mean(dec_rrab.radian) + np.mean(dec_rrc.radian) + np.mean(dec_cepfu_1.radian) + np.mean(dec_cepfo.radian)) / 4

ra_cepfu = np.delete(ra_cepfu_1.radian, fu_exception)
dec_cepfu = np.delete(dec_cepfu_1.radian, fu_exception)

x_fu = -r_fu * np.sin(ra_cepfu - ra_0) * np.cos(dec_cepfu)
y_fu = r_fu * np.sin(dec_cepfu)*np.cos(dec_0) - r_fu * np.sin(dec_0) * np.cos(ra_cepfu - ra_0) * np.cos(dec_cepfu)
z_fu = r_lmc_0 - r_fu * np.sin(dec_cepfu) * np.sin(dec_0) - r_fu * np.cos(dec_0)*np.cos(ra_cepfu) - ra_0 * np.cos(dec_cepfu)

x_fo = -r_fo * np.sin(ra_cepfo.radian - ra_0) * np.cos(dec_cepfo.radian)
y_fo = r_fo * np.sin(dec_cepfo.radian)*np.cos(dec_0) - r_fo * np.sin(dec_0) * np.cos(ra_cepfo.radian - ra_0) * np.cos(dec_cepfo.radian)
z_fo = r_lmc_0 - r_fo * np.sin(dec_cepfo.radian) * np.sin(dec_0) - r_fo * np.cos(dec_0)*np.cos(ra_cepfo.radian) - ra_0 * np.cos(dec_cepfo.radian)

x_ab = -r_ab * np.sin(ra_rrab.radian - ra_0) * np.cos(dec_rrab.radian)
y_ab = r_ab * np.sin(dec_rrab.radian)*np.cos(dec_0) - r_ab * np.sin(dec_0) * np.cos(ra_rrab.radian - ra_0) * np.cos(dec_rrab.radian)
z_ab = r_lmc_0 - r_ab * np.sin(dec_rrab.radian) * np.sin(dec_0) - r_ab * np.cos(dec_0)*np.cos(ra_rrab.radian) - ra_0 * np.cos(dec_rrab.radian)

x_c = -r_c * np.sin(ra_rrc.radian - ra_0) * np.cos(dec_rrc.radian)
y_c = r_c * np.sin(dec_rrc.radian)*np.cos(dec_0) - r_c * np.sin(dec_0) * np.cos(ra_rrc.radian - ra_0) * np.cos(dec_rrc.radian)
z_c = r_lmc_0 - r_c * np.sin(dec_rrc.radian) * np.sin(dec_0) - r_c * np.cos(dec_0)*np.cos(ra_rrc.radian) - ra_0 * np.cos(dec_rrc.radian)



    











In [108]:



    


fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.plot(x_ab, y_ab, z_ab, linestyle='None' ,marker='.', color='green',ms=1., label='RR AB')
#plt.plot(x_ab, y_ab,z_ab-( -0.98*x_ab + 0.29*y_ab+2.68), linestyle='None' ,marker='.',color='red',ms=1.)
plt.plot(x_c, y_c, z_c, linestyle='None' ,marker='.', color='black',ms=1., label='RR C')
plt.plot(x_fu, y_fu, z_fu, linestyle='None' ,marker='.', color='blue',ms=1., label='Cep FU')
plt.plot(x_fo, y_fo, z_fo, linestyle='None' ,marker='.', color='gold',ms=1., label='Cep FO')
ax.set_xlabel('x [kpc]')
ax.set_ylabel('y [kpc]')
ax.set_zlabel('z [kpc]')
ax.set_title('Distribuição das Estrelas na LMC')
#plt.legend()
plt.savefig('plot3d.png', bbox_inches='tight')
#ax.set_xticks([-5000, -4000, -3000, -2000, -1000, 0, 1000, 2000, 3000, 4000])
#ax.set_xticklabels([-5, -4, -3, -2, -1, 0, 1, 2, 3, 4])
#ax.set_yticks([-3000, -2000, -1000, 0, 1000, 2000, 3000])
#ax.set_yticklabels([ -3, -2, -1, 0, 1, 2, 3])
#ax.set_zticks([-2000, 0, 2000, 4000, 6000, 8000])
#ax.set_zticklabels([-2, 0, 2, 4, 6, 8])



    


















    
























In [80]:



    


#fu_plot = x_fu, y_fu, z_fu
#fo_plot = x_fo, y_fo, z_fo
#ab_plot = x_ab, y_ab, z_ab
#c_plot = x_c, y_c, z_c

np.savetxt('fu_plot.txt', np.dstack((x_fu, y_fu, z_fu))[0], fmt='%s')
np.savetxt('fo_plot.txt', np.dstack((x_fo, y_fo, z_fo))[0], fmt='%s')
np.savetxt('ab_plot.txt', np.dstack((x_ab, y_ab, z_ab))[0], fmt='%s')
np.savetxt('c_plot.txt', np.dstack((x_c, y_c, z_c))[0], fmt='%s')



    











In [81]:



    


#plt.imshow(np.dstack((ra_rrab.degree, dec_rrab.degree))[0], interpolation='Nearest')
ra_all_1 = np.insert(ra_rrab.degree, [1],ra_rrc.degree)#, ra_cepfo.degree)
ra_all_2 = np.insert(ra_all_1, [1],ra_cepfu_1.degree)
ra_all = np.insert(ra_all_2, [1], ra_cepfo.degree)
dec_all_1 = np.insert(dec_rrab.degree, [1],dec_rrc.degree)#, ra_cepfo.degree)
dec_all_2 = np.insert(dec_all_1, [1],dec_cepfu_1.degree)
dec_all = np.insert(dec_all_2, [1], dec_cepfo.degree)
#dec_all_1 = np.append(dec_rrab.degree, dec_rrc.degree, dec_cepfo.degree)
#dec_all = np.apend(dec_all_1, dec_cepfu_1.degree)
h, xedges, yedges = np.histogram2d(dec_all, ra_all, bins=[30,30])#,
                                   #range=[[np.min(ra_rrab.degree), np.max(ra_rrab.degree)],[np.min(dec_rrab.degree), np.max(dec_rrab.degree)]])
#extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
plt.imshow(h, interpolation='Nearest')#, range=extent)
plt.colorbar()
#plt.gca().invert_yaxis



    


















    
Out[81]:








<matplotlib.colorbar.Colorbar at 0x7fd296ddef98>










    
























In [82]:



    


x_all_1 = np.insert(x_c, [0],x_ab)#, ra_cepfo.degree)
x_all_2 = np.insert(x_all_1, [0], x_fo)
x_all = np.insert(x_all_2, [0], x_fu)

y_all_1 = np.insert(y_c, [0],y_ab)#, ra_cepfo.degree)
y_all_2 = np.insert(y_all_1, [0], y_fo)
y_all = np.insert(y_all_2, [0], y_fu)

z_all_1 = np.insert(z_c, [0],z_ab)#, ra_cepfo.degree)
z_all_2 = np.insert(z_all_1, [0], z_fo)
z_all = np.insert(z_all_2, [0], z_fu)

print(len(x_all), len(y_all), len(z_all))
A_all = np.column_stack((np.ones(x_all.size), x_all, y_all))
c_all, resid_all, rank_all, sigma_all = np.linalg.lstsq(A_all,z_all)
print(c_all)
#A_rrlyrab = np.column_stack((np.ones(x_ab.size), x_ab, y_ab))
#c_all, resid_all, rank_all, sigma_all = np.linalg.lstsq(A_rrlyrab,z_ab)



    


















    






25707 25707 25707
[ 2.68048019 -0.98602693  0.29295402]

















In [83]:



    


c_all
inclination = np.arccos(1 / np.sqrt(1 + c_all[0]**2 + c_all[1]**2))
theta = np.arctan( -c_all[0] / c_all[1]) + np.sign(c_all[2]) * np.pi / 2



    











In [84]:



    


print(np.rad2deg(theta), np.rad2deg(inclination))



    


















    






159.803723116 70.7033388092

















In [85]:



    


from matplotlib.patches import Ellipse
elli = Ellipse(xy=[0], width=[2], height=[1], angle=[np.rad2deg(inclination)])
fig = plt.figure()
ax = fig.add_subplot(111)#, aspect='equal')
ax.add_artist(elli)
plt.plot(x_all, y_all, '.', ms=0.5)
#plt.gca().add_artist(elli)

#plt.plot(c_all[0]*x_all - c_all[2], c_all[1]*y_all - c_all[2], '.', color='green')
plt.show(True)



    


















    








<matplotlib.figure.Figure at 0x7fd27c85f710>

















In [100]:



    


#import sklearn.linear_model as lm

#ll = lm.LinearRegression()
#ll.fit(x_all, y_all, z_all)

import scipy.linalg

A_al = np.column_stack((x_all, y_all, np.ones(x_all.size)))
C,_,_,_ = scipy.linalg.lstsq(A_al, z_all)
#A_al = np.column_stack((x_ab, y_ab, np.ones(x_ab.size)))
#C,_,_,_ = scipy.linalg.lstsq(A_al, z_ab)



    











In [101]:



    


C



    


















    
Out[101]:








array([-0.98602693,  0.29295402,  2.68048019])

















In [102]:



    


inclination = np.arccos(1 / np.sqrt(1 + C[0]**2 + C[1]**2))
theta = np.arctan( -C[0] / C[1]) + np.sign(C[2]) * np.pi / 2



    











In [103]:



    


print(np.rad2deg(theta), np.rad2deg(inclination))



    


















    






163.453026208 45.8084447759

















In [143]:



    


angles = np.linspace(0,360,21)[:-1] # A list of 20 angles between 0 and 360
 
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.plot(x_ab, y_ab, z_ab, linestyle='None' ,marker='.', color='green',ms=1.)
#plt.plot(x_ab, y_ab,z_ab-( -0.98*x_ab + 0.29*y_ab+2.68), linestyle='None' ,marker='.',color='red',ms=1.)
plt.plot(x_c, y_c, z_c, linestyle='None' ,marker='.', color='black',ms=1.)
plt.plot(x_fu, y_fu, z_fu, linestyle='None' ,marker='.', color='blue',ms=1.)
plt.plot(x_fo, y_fo, z_fo, linestyle='None' ,marker='.', color='gold',ms=1.)
#plt.plot(x_all, y_all, C[0]*x_all+C[1]*y_all+C[2], linestyle='None', marker='.', color='r', ms=1.)

ax.set_xlabel('x [kpc]')
ax.set_ylabel('y [kpc]')
ax.set_zlabel('z [kpc]')
ax.set_title('Distribuição das Estrelas na LMC')

#for ii in angles:
#    ax.view_init(0., azim=ii)
#    plt.savefig("movie"%ii+".png", bbox_inches='tight')
ax.view_init(30, 360)
plt.savefig('3d_plot_360.png',  bbox_inches='tight')



    


















    
























In [ ]:

Content source: glauffer/final_project

Similar notebooks:

notebook.community | gallery | about