In [1]:
from astroML.plotting import hist as amlhist
from astropy.coordinates import SkyCoord
from scipy.stats import ks_2samp as ks
import astropy.coordinates as coord
import matplotlib.pyplot as plt
from astropy import units as u
from astropy.io import ascii
import pandas as pd
import numpy as np
import glob
import os
%pylab inline
In [2]:
path = '/Users/apple/Downloads/VVV_phot/RRab/'
os.chdir(path)
os.system('rm RRab.dat')
tiles = sorted(glob.glob(path+'RRab*.dat'))
os.system('rm lcs/b*_*.dat')
os.system('rm lcs/b*_*.pdf')
os.system('rm o/lcs/*.dat')
os.system('rm o/lcs/*.pdf')
os.system('sort -n RRab_*.dat > RRab.dat');
In [3]:
total_ab = 0
for tile in tiles:
lcs = 0
f = open(tile,'r')
for RRab in f:
if RRab.split()[0] == '#Num':
continue
if 'b213' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b213_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b205' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b205_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b214' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b214_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b218' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b218_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b227' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b227_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b228' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b228_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b208' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b208_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b209' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b209_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b210' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b210_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b211' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b211_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b212' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b212_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b215' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b215_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b219' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b219_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b220' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b220_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b221' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b221_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b222' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b222_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b223' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b223_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b224' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b224_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b225' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b225_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
if 'b226' in tile:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b226_'+RRab.split()[0]+'.dat lcs/'+
tile.replace('/RRab/','/RRab/lcs/').split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
else:
os.system('cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b201_'+RRab.split()[0]+'.dat lcs/'+
tile.split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat')
#print 'cp ../'+tile.split('RRab_')[1].split('.dat')[0]+'/b201_'+RRab.split()[0]+'.dat lcs/'+tile.split('RRab_')[1].split('.dat')[0]+'_'+RRab.split()[0]+'.dat'
lcs = lcs + 1
print tile.split('RRab_')[1].split('.dat')[0]+':', lcs, 'RRab'
total_ab = total_ab + lcs
f.close()
print 'Total number of RRab:', total_ab, 'in', len(tiles), 'tiles'
In [4]:
#Extraccion de repetidas
overlap_names = np.genfromtxt('o/overlap.dat', usecols=(0), dtype=None)
overlap_indices = np.genfromtxt('o/overlap.dat', usecols=(17), dtype=int)
overlap_P, overlap_AKs, overlap_Ks = np.genfromtxt('o/overlap.dat', usecols=(11,9,4), unpack=True)
overlap_RA, overlap_DEC = np.genfromtxt('o/overlap.dat', usecols=(2,3), unpack=True)
hist_oP = np.empty_like(overlap_P)
hist_oAKs = np.empty_like(overlap_P)
hist_oKs = np.empty_like(overlap_P)
for i in np.unique(overlap_indices):
hist_oP[i] = overlap_P[overlap_indices == i][0] - overlap_P[overlap_indices == i][1]
hist_oAKs[i] = overlap_AKs[overlap_indices == i][0] - overlap_AKs[overlap_indices == i][1]
hist_oKs[i] = overlap_Ks[overlap_indices == i][0] - overlap_Ks[overlap_indices == i][1]
for i in np.unique(overlap_indices):
os.system('cat lcs/'+overlap_names[overlap_indices == i][0].replace('-','_')+'.dat lcs/'+overlap_names[overlap_indices == i][1].replace('-','_')+'.dat >> o/lcs/'+overlap_names[overlap_indices == i][0].replace('-','_')+'_'+overlap_names[overlap_indices == i][1].replace('-','_')+'.dat')
os.system('rm lcs/'+overlap_names[overlap_indices == i][0].replace('-','_')+'.dat lcs/'+overlap_names[overlap_indices == i][1].replace('-','_')+'.dat')
# os.system('cp o/lcs/'+overlap_names[overlap_indices == i][0].replace('-','_')+'_'+overlap_names[overlap_indices == i][1].replace('-','_')+'.dat lcs/')
In [5]:
!wc -l o/lcs/*
In [6]:
def phase(t,P):
return (t)/P - np.floor((t)/P)
In [7]:
#plt.figure(figsize=(8,4.5))
plt.figure(figsize=(12,5))
plt.subplot(121)
mjd, mag, emag = np.genfromtxt('o/lcs/b208_303363_b222_398211.dat', usecols=(0,1,2), unpack=True)
P_b208b222 = 0.640578
plt.errorbar(phase(mjd, P_b208b222), mag, emag, fmt='og', ecolor='b')
plt.errorbar(phase(mjd, P_b208b222)+1, mag, emag, fmt='og', ecolor='b')
plt.gca().invert_yaxis()
plt.ylim(14.00,13.58)
plt.xlim(-0.1,2.1)
plt.xlabel('Phase', size=14)
plt.ylabel(r'$\rm{K}_{\rm s}$', size=14)
plt.title(r'VVV J2762905.51-320926.8 ($N_{\rm{epochs}}=124$)', size=13)
plt.subplot(322)
plt.hist(hist_oP[1:33], histtype='stepfilled', lw=2.5, color='grey', alpha=0.4, bins=10)
plt.xlabel(r'$\Delta P$ (days)', size=14)
plt.yticks([0,5,10,15])
plt.xticks([-0.00004,-0.00002,0.0,0.00002,0.00004],[r'$4\times 10^{-4}$',r'$2\times 10^{-4}$',r'$0.0$',r'$2\times 10^{-4}$',r'$4\times 10^{-4}$'])
plt.grid(True)
plt.subplot(324)
plt.hist(hist_oAKs[1:33], histtype='stepfilled', lw=2.5, color='grey', alpha=0.4, bins=10)
plt.xlabel(r'$\Delta A_{\rm{K}_{\rm s}}$ (mag)', size=14)
plt.yticks([0,5,10,15])
plt.grid(True)
plt.xlim(-0.22,0.12)
plt.subplot(326)
plt.hist(hist_oKs[1:33], histtype='stepfilled', lw=2.5, color='grey', alpha=0.4, bins=10)
plt.xlabel(r'$\Delta \rm{K}_{\rm s}$ (mag)', size=14)
plt.yticks([0,5,10,15])
plt.grid(True)
plt.xlim(-0.062,0.042)
plt.tight_layout()
plt.savefig('overlap.pdf',format='pdf')
plt.savefig('/Users/apple/Dropbox/Magister/Cursos/InvestigacionB/Latex_RR/overlap.pdf',format='pdf')
In [8]:
ID = []
tile_name = []
P = np.array([])
RA = np.array([])
DEC = np.array([])
for tile in tiles:
with open(tile,'r') as f:
for ab in f:
if ab.split()[0] == '#Num':
continue
if ab.split()[0] in np.array_str(overlap_names):
continue
tile_name.append( tile.split('RRab_')[1].split('.dat')[0])
ID.append(str(ab.split()[0]))
P = np.append(P, float(ab.split()[10]))
RA = np.append(RA, float(ab.split()[1]))
DEC = np.append(DEC, float(ab.split()[2]))
ID_overlap = ['overlap']*len(overlap_names[0::2])
P_overlap = overlap_P[0::2]
tile_name_overlap = glob.glob('o/lcs/*.dat')
In [9]:
#Respaldo 1
#ID = []
#tile_name = []
#P = np.array([])
#RA = np.array([])
#DEC = np.array([])
#for tile in tiles:
# with open(tile,'r') as f:
# for ab in f:
# if ab.split()[0] == '#Num':
# continue
# tile_name.append( tile.split('RRab_')[1].split('.dat')[0])
# ID.append( str(ab.split()[0] ) )
# P = np.append(P, float(ab.split()[10]) )
# RA = np.append(RA, float(ab.split()[1]) )
# DEC = np.append(DEC, float(ab.split()[2]) )
In [10]:
out = open('target.lis','w')
for i in xrange(len(ID)):
if i == 500:
out.close()
os.system('./tff');
os.system('mv dff.dat dff_1.dat');
out = open('target.lis','w')
out.write( tile_name[i]+'-'+ID[i]+'\n'+ str(P[i])+'\n'+path+'lcs/'+tile_name[i]+'_'+ID[i]+'.dat \n' )
out.close()
os.system('./tff');
os.system('mv dff.dat dff_2.dat');
out = open('target.lis','w')
for i in xrange(len(ID_overlap)):
out.write( tile_name_overlap[i].split('lcs/')[1].split('.dat')[0].replace('_','-')+'\n'+str(overlap_P[0::2][i])+'\n'+path+tile_name_overlap[i]+'\n' )
out.close()
os.system('./tff');
os.system('mv dff.dat dff_overlap.dat');
os.system('cat dff_1.dat dff_2.dat >> dff.dat');
os.system('rm dff_1.dat dff_2.dat');
In [11]:
names_ab = []
avg_mag_ab, A1_ab, Phi1_ab, A2_ab, Phi2_ab, A3_ab, Phi3_ab, A4_ab, Phi4_ab, A5_ab, Phi5_ab, A6_ab, Phi6_ab = [np.zeros(len(ID)) for d in range(13)]
ID_ab = []
F = open('/Users/apple/Downloads/VVV_phot/RRab/dff.dat','r').readlines()
c = 0
for i in range(0,len(F),5):
ID_ab.append( F[i].split()[0].split('-')[1] )
names_ab.append( F[i].split()[0] )
avg_mag_ab[c] = float( F[i+1].split()[2] )
A1_ab[c] = float( F[i+2].split()[0] )
Phi1_ab[c] = float( F[i+2].split()[1] )
A2_ab[c] = float( F[i+2].split()[2] )
if A2_ab[c] == 0.0:
c = c + 1
continue
Phi2_ab[c] = float( F[i+2].split()[3] )
A3_ab[c] = float( F[i+2].split()[4] )
if A3_ab[c] == 0.0:
c = c + 1
continue
Phi3_ab[c] = float( F[i+2].split()[5] )
A4_ab[c] = float( F[i+2].split()[6] )
if A4_ab[c] == 0.0:
c = c + 1
continue
Phi4_ab[c] = float( F[i+2].split()[7] )
A5_ab[c] = float( F[i+2].split()[8] )
if A5_ab[c] == 0.0:
c = c + 1
continue
Phi5_ab[c] = float( F[i+2].split()[9] )
A6_ab[c] = float( F[i+3].split()[0] )
if A6_ab[c] == 0.0:
c = c + 1
continue
Phi6_ab[c] = float( F[i+3].split()[1] )
c = c + 1
In [12]:
names_ab_o = []
avg_mag_ab_o, A1_ab_o, Phi1_ab_o, A2_ab_o, Phi2_ab_o, A3_ab_o, Phi3_ab_o, A4_ab_o, Phi4_ab_o, A5_ab_o, Phi5_ab_o, A6_ab_o, Phi6_ab_o = [np.zeros(len(ID_overlap)) for d in range(13)]
ID_ab_o = []
F = open('/Users/apple/Downloads/VVV_phot/RRab/dff_overlap.dat','r').readlines()
c = 0
for i in range(0,len(F),5):
ID_ab_o.append( F[i].split()[0].split('-')[1] )
names_ab_o.append( F[i].split()[0] )
avg_mag_ab_o[c] = float( F[i+1].split()[2] )
A1_ab_o[c] = float( F[i+2].split()[0] )
Phi1_ab_o[c] = float( F[i+2].split()[1] )
A2_ab_o[c] = float( F[i+2].split()[2] )
if A2_ab_o[c] == 0.0:
c = c + 1
continue
Phi2_ab_o[c] = float( F[i+2].split()[3] )
A3_ab_o[c] = float( F[i+2].split()[4] )
if A3_ab_o[c] == 0.0:
c = c + 1
continue
Phi3_ab_o[c] = float( F[i+2].split()[5] )
A4_ab_o[c] = float( F[i+2].split()[6] )
if A4_ab_o[c] == 0.0:
c = c + 1
continue
Phi4_ab_o[c] = float( F[i+2].split()[7] )
A5_ab_o[c] = float( F[i+2].split()[8] )
if A5_ab_o[c] == 0.0:
c = c + 1
continue
Phi5_ab_o[c] = float( F[i+2].split()[9] )
A6_ab_o[c] = float( F[i+3].split()[0] )
if A6_ab_o[c] == 0.0:
c = c + 1
continue
Phi6_ab_o[c] = float( F[i+3].split()[1] )
c = c + 1
In [13]:
def fourier(x,*p):
#Fourier
f = p[0]
c = 1
N = (len(p)-1)/2
for i in range(1, 2*N, 2):
f = f + p[i]*np.sin(2*np.pi*c*x + p[i+1])
c = c + 1
return f
In [14]:
data_ab = [names_ab, avg_mag_ab, A1_ab, Phi1_ab, A2_ab, Phi2_ab, A3_ab, Phi3_ab, A4_ab, Phi4_ab, A5_ab, Phi5_ab, A6_ab, Phi6_ab]
amp = []
name_amp = []
is_plot = False
for j in range(len(names_ab)):
mjd, mag, emag = np.genfromtxt('lcs/'+data_ab[0][j].replace('-','_')+'.dat', usecols=(0,1,2), unpack=True)
popt = []
for i in range(1,len(data_ab)):
if data_ab[i][j] == 0.0:
break
popt.append(data_ab[i][j])
time = np.linspace(0,2,200)
amp.append( np.max(fourier(time, *popt)) - np.min(fourier(time, *popt)) )
name_amp.append('lcs/'+data_ab[0][j].replace('-','_')+'.jpg')
if not is_plot:
continue
plt.plot(time, fourier(time, *popt), color='r');
plt.errorbar(phase(mjd,float(P[j])), mag, yerr=emag, fmt='o',color='g', ecolor='b');
plt.errorbar(phase(mjd,float(P[j]))+1, mag, yerr=emag, fmt='o',color='g', ecolor='b');
plt.title(data_ab[0][j] + ' P = '+str(round(float(P[j]),3))+' d', fontsize=13)
plt.ylim(max(mag)+0.1,max(mag)-0.5);
plt.xlim(-0.1,2.1);
plt.xlabel('Phase', size=14)
plt.ylabel(r'$K_s$', size=14)
plt.savefig('lcs/'+data_ab[0][j].replace('-','_')+'.jpg', format='jpg', dpi=300)
plt.close()
In [15]:
if not is_plot:
sorted_amp = zip(*sorted(zip(amp, name_amp)))[0]
sorted_lcs = zip(*sorted(zip(amp, name_amp)))[1]
c = 1
for i in sorted_lcs:
os.system('cp %s lcs/%d.jpg' %(i,c))
c = c + 1
In [16]:
data_ab_o = [names_ab_o, avg_mag_ab_o, A1_ab_o, Phi1_ab_o, A2_ab_o, Phi2_ab_o, A3_ab_o, Phi3_ab_o, A4_ab_o, Phi4_ab_o, A5_ab_o, Phi5_ab_o, A6_ab_o, Phi6_ab_o]
amp_o = []
is_plot_o = False
for j in range(len(names_ab_o)):
mjd, mag, emag = np.genfromtxt(tile_name_overlap[j], usecols=(0,1,2), unpack=True)
popt = []
for i in range(1,len(data_ab_o)):
if data_ab_o[i][j] == 0.0:
break
popt.append(data_ab_o[i][j])
time = np.linspace(0,2,200)
amp_o.append( np.max(fourier(time, *popt)) - np.min(fourier(time, *popt)) )
if not is_plot_o:
continue
plt.plot(time, fourier(time, *popt), color='r');
plt.errorbar(phase(mjd,float(P_overlap[j])), mag, yerr=emag, fmt='o',color='g', ecolor='b');
plt.errorbar(phase(mjd,float(P_overlap[j]))+1, mag, yerr=emag, fmt='o',color='g', ecolor='b');
plt.title(data_ab_o[0][j] + ' P = '+str(round(float(P_overlap[j]),3))+' d', fontsize=13)
plt.ylim(max(mag)+0.1,min(mag)-0.1);
plt.xlim(-0.1,2.1);
plt.xlabel('Phase', size=14)
plt.ylabel(r'$K_s$', size=14)
plt.savefig(tile_name_overlap[j].replace('dat','pdf'), format='pdf')
plt.close()
In [17]:
data_ab_fourier = [names_ab, P, np.round(avg_mag_ab,3), np.round(A1_ab,3), np.round(Phi1_ab,3),
np.round(A2_ab/A1_ab,3), np.round(Phi2_ab-2*Phi1_ab,3), np.round(A3_ab/A1_ab,3),
np.round(Phi3_ab-3*Phi1_ab,3), np.round(A4_ab/A1_ab,3), np.round(Phi4_ab-4*Phi1_ab,3),
np.round(A5_ab/A1_ab,3), np.round(Phi5_ab-5*Phi1_ab,3), np.round(A6_ab/A1_ab,3), np.round(Phi6_ab-6*Phi1_ab,3)]
write = np.column_stack((names_ab, amp, P, np.round(avg_mag_ab,3), np.round(A1_ab,3), np.round(Phi1_ab,3),
np.round(A2_ab/A1_ab,3), np.round(Phi2_ab-2*Phi1_ab,3), np.round(A3_ab/A1_ab,3),
np.round(Phi3_ab-3*Phi1_ab,3), np.round(A4_ab/A1_ab,3), np.round(Phi4_ab-4*Phi1_ab,3),
np.round(A5_ab/A1_ab,3), np.round(Phi5_ab-5*Phi1_ab,3), np.round(A6_ab/A1_ab,3), np.round(Phi6_ab-6*Phi1_ab,3)))
np.savetxt('fourier.dat', write, newline="\n", fmt='%s',
header='names \t amp \t P \t A0 \t A1 \t Phi1 \t A21 \t Phi21 \t A31 \t Phi31 \t A41 \t Phi41 \t A51 \t Phi51 \t A61 \t Phi61')
data_ab_fourier_o = [names_ab_o, P_overlap, np.round(avg_mag_ab_o,3), np.round(A1_ab_o,3), np.round(Phi1_ab_o,3),
np.round(A2_ab_o/A1_ab_o,3), np.round(Phi2_ab_o-2*Phi1_ab_o,3), np.round(A3_ab_o/A1_ab_o,3),
np.round(Phi3_ab_o-3*Phi1_ab_o,3), np.round(A4_ab_o/A1_ab_o,3), np.round(Phi4_ab_o-4*Phi1_ab_o,3),
np.round(A5_ab_o/A1_ab_o,3), np.round(Phi5_ab_o-5*Phi1_ab_o,3), np.round(A6_ab_o/A1_ab_o,3), np.round(Phi6_ab_o-6*Phi1_ab_o,3)]
write = np.column_stack(([aux[:aux[5:].find('-')+5] for aux in names_ab_o], amp_o, P_overlap, np.round(avg_mag_ab_o,3), np.round(A1_ab_o,3), np.round(Phi1_ab_o,3),
np.round(A2_ab_o/A1_ab_o,3), np.round(Phi2_ab_o-2*Phi1_ab_o,3), np.round(A3_ab_o/A1_ab_o,3),
np.round(Phi3_ab_o-3*Phi1_ab_o,3), np.round(A4_ab_o/A1_ab_o,3), np.round(Phi4_ab_o-4*Phi1_ab_o,3),
np.round(A5_ab_o/A1_ab_o,3), np.round(Phi5_ab_o-5*Phi1_ab_o,3), np.round(A6_ab_o/A1_ab_o,3), np.round(Phi6_ab_o-6*Phi1_ab_o,3)))
np.savetxt('fourier_overlap.dat', write, newline="\n", fmt='%s',
header='names \t amp \t P \t A0 \t A1 \t Phi1 \t A21 \t Phi21 \t A31 \t Phi31 \t A41 \t Phi41 \t A51 \t Phi51 \t A61 \t Phi61')
plt.figure(figsize=(8,8))
plt.subplot(411)
plt.scatter(P, data_ab_fourier[5], color='r', marker='*', s=80, alpha=0.4)
plt.scatter(P_overlap, data_ab_fourier_o[5], color='b', marker='o', s=40, alpha=0.6)
#plt.xlabel('Period (d)', fontsize=16)
plt.ylabel(r'$R_{21}$', fontsize=16)
plt.ylim(-0.05, )
plt.grid()
plt.subplot(412)
plt.scatter(P[data_ab_fourier[6] != 0.], np.mod(data_ab_fourier[6][data_ab_fourier[6] != 0.],np.pi), color='r', marker='*', s=80, alpha=0.4)
plt.scatter(P_overlap[data_ab_fourier_o[6] != 0.], np.mod(data_ab_fourier_o[6][data_ab_fourier_o[6] != 0.],np.pi), color='b', marker='o', s=40, alpha=0.6)
#plt.xlabel('Period (d)', fontsize=16)
plt.ylabel(r'$\phi_{21}$', fontsize=16)
plt.ylim(-0.05, )
plt.grid()
plt.subplot(413)
plt.scatter(np.asfarray(P)[data_ab_fourier[7] != 0.], data_ab_fourier[7][data_ab_fourier[7] != 0.] , color='r', marker='*', s=80, alpha=0.4)
plt.scatter(np.asfarray(P_overlap)[data_ab_fourier_o[7] != 0.], data_ab_fourier_o[7][data_ab_fourier_o[7] != 0.] , color='b', marker='o', s=40, alpha=0.6)
#plt.xlabel('Period (d)', fontsize=16)
plt.ylabel(r'$R_{31}$', fontsize=16)
plt.ylim(-0.05, )
plt.grid()
plt.subplot(414)
plt.scatter(np.asfarray(P)[data_ab_fourier[8] != 0.], np.mod(data_ab_fourier[8][data_ab_fourier[8] != 0.],np.pi), color='r', marker='*', s=80, alpha=0.4)
plt.scatter(np.asfarray(P_overlap)[data_ab_fourier_o[8] != 0.], np.mod(data_ab_fourier_o[8][data_ab_fourier_o[8] != 0.],np.pi), color='b', marker='o', s=40, alpha=0.6)
plt.xlabel(r'$P$ (days)', fontsize=16)
plt.ylabel(r'$\phi_{31}$', fontsize=16)
plt.ylim(-0.05, )
plt.grid()
plt.tight_layout()
plt.savefig('/Users/apple/Downloads/VVV_phot/RRab/Fourier.pdf', format='pdf')
plt.savefig('/Users/apple/Dropbox/Magister/Cursos/InvestigacionB/Latex_RR/Fourier.pdf', format='pdf')
plt.show()
In [18]:
def OoI(x):
return 0.32*(-2.627-22.046*x-30.876*x**2)**(2/3.)
def OoII(x):
return ((0.064-2.481*x+10.345*x**3)/2.6)**(2/3.)
In [19]:
P_all = np.concatenate((P, P_overlap))
plt.figure(figsize=(8,5))
plt.subplot(211)
plt.scatter(P, amp, marker='*', color='r', label='RRab', s=70., alpha=0.5)
plt.scatter(P_overlap, amp_o, marker='o', color='b', label='RRab', s=40., alpha=0.6)
plt.ylabel(r'$\Delta K_{\rm{s}}$ (mag)', fontsize=14)
plt.xlim(0.35,0.95)
plt.ylim(0.12,0.51)
plt.grid()
plt.subplot(212)
amlhist(np.asfarray(P_all), bins='blocks', histtype='stepfilled', lw=2.5, color='grey', alpha=0.4)
plt.xlabel(r'$\log{P}$ (days)', fontsize=15)
plt.ylabel(r'# RRab', fontsize=14)
plt.xlim(0.35,0.95)
plt.ylim(0,)
plt.grid()
plt.tight_layout()
plt.savefig('/Users/apple/Downloads/VVV_phot/RRab/Bailey.pdf', format='pdf')
plt.savefig('/Users/apple/Dropbox/Magister/Cursos/InvestigacionB/Latex_RR/Bailey.pdf', format='pdf')
plt.show()
In [20]:
P_all = np.concatenate((P, P_overlap))
plt.figure(figsize=(8,5))
plt.scatter(np.log10(P), amp, marker='*', color='r', s=70., alpha=0.5)
plt.scatter(np.log10(P_overlap), amp_o, marker='o', color='b', s=40., alpha=0.6)
#Lineas de OoI
x = np.linspace(-0.24,-0.05)
plt.plot(x, OoII(x), color='k', ls='dashed', lw=3, alpha=0.8, label='Oo II')
#Lineas de OoII
x = np.linspace(-0.38,-0.18)
plt.plot(x, OoI(x), color='k', ls='-', lw=4, alpha=0.5, label='Oo I')
x = np.linspace(-0.3,-0.1)
plt.plot(x-0.05, 0.32*(-7.007*x-0.343)**(2/3.), color='g', ls='dashed', lw=3, label='Oo division')
plt.ylabel(r'$\Delta K_{\rm{s}}$ (mag)', fontsize=14)
plt.xlabel(r'$\log{P}$ (days)', fontsize=14)
plt.xlim(-0.48,0)
plt.ylim(0.14,0.51)
plt.legend()
plt.savefig('/Users/apple/Downloads/VVV_phot/RRab/Bailey_Oo.pdf', format='pdf')
plt.savefig('/Users/apple/Dropbox/Magister/Cursos/InvestigacionB/Latex_RR/Bailey_Oo.pdf', format='pdf')
plt.show()
In [21]:
#HASP percentage
print 'N_HASP / N_RRL = ', len(P_all[P_all < 0.45]), '/', len(P_all), '*100 = ', len(P_all[P_all < 0.45])/float(len(P_all))*100
In [22]:
out = open('coo_beam_p.dat','w')
for i in np.linspace(0.0,11.0,300):
for j in np.linspace(-8.0,-10.8,80):
out.write( str(i)+'\t'+str(j)+'\t2\t1\n' )
out = open('coo_beam_n.dat','w')
for i in np.linspace(0.0,-11.0,300):
for j in np.linspace(-8.0,-10.8,80):
out.write( str(i)+'\t'+str(j)+'\t2\t1\n' )
In [23]:
def HMS(hours, pos):
minutes = np.abs ((hours-int(hours) ) * 60.)
seconds = 0
hours = int(hours)
if hours == 0:
if minutes == 0:
if seconds == 0:
return r"$%d^\circ%02d'$" % (hours, minutes)
return "%d''" % (seconds)
return "%d'%02d''" % (minutes, seconds)
return r"$%d^\circ%02d'$" % (hours, minutes)
plt.figure(figsize=(16,1.6))
plt.gca().xaxis.set_major_formatter(plt.FuncFormatter(HMS))
plt.gca().yaxis.set_major_formatter(plt.FuncFormatter(HMS))
l, b, AKs = np.genfromtxt('/Users/apple/Downloads/VVV_phot/RRab/coo_beam_n_out.dat', usecols=(0,1,5), unpack=True)
l_2, b_2, AKs_2 = np.genfromtxt('/Users/apple/Downloads/VVV_phot/RRab/coo_beam_p_out.dat', usecols=(0,1,5), unpack=True)
AB = SkyCoord(ra=RA, dec=DEC, unit=(u.degree, u.degree), frame='icrs')
AB_o = SkyCoord(ra=overlap_RA, dec=overlap_DEC, unit=(u.degree, u.degree), frame='icrs')
gal_ab = AB.galactic
gal_ab_o = AB_o.galactic
l_ab = np.array([])
l_ab_o = np.array([])
for i,ell in enumerate(gal_ab.l.deg):
if ell > 50.:
l_ab = np.append(l_ab, ell - 360.)
else:
l_ab = np.append(l_ab, ell)
for i,ell in enumerate(gal_ab_o.l.deg):
if ell > 50.:
l_ab_o = np.append(l_ab_o, ell - 360.)
else:
l_ab_o = np.append(l_ab_o, ell)
b_ab = gal_ab.b.deg
b_ab_o = gal_ab_o.b.deg
a = plt.hexbin(l, b, AKs, gridsize=80, cmap=plt.get_cmap('gist_yarg'), vmax=0.2)
a = plt.hexbin(l_2, b_2, AKs_2, gridsize=80, cmap=plt.get_cmap('gist_yarg'), vmax=0.2)
plt.scatter(l_ab, b_ab, marker='*', color='r', label='RRab', s=10.)
plt.scatter(l_ab_o, b_ab_o, marker='o', color='b', label='RRab', s=10.)
cb = plt.colorbar(a, pad=0.01, shrink=1.0, ticks=[0, 0.05, 0.10, 0.15, 0.2], aspect=5)
cb.set_label(r'$A_{\rm{K}_{\rms}}$', fontsize=19)
plt.xticks(np.arange(-10,11,2.5))
plt.yticks([-8,-9,-10])
plt.xlabel(r'$\ell$', fontsize=19)
plt.ylabel('$b$', fontsize=19)
plt.xlim(11,-10.2)
plt.ylim(-10.35,-7.9)
plt.tight_layout()
plt.savefig('/Users/apple/Downloads/VVV_phot/RRab/lb.pdf', format='pdf')
plt.savefig('/Users/apple/Dropbox/Magister/Cursos/InvestigacionB/Latex_RR/lb.pdf', format='pdf')
plt.show()
In [24]:
plt.figure(figsize=(8,4))
amlhist(np.concatenate((avg_mag_ab,avg_mag_ab_o)), bins='blocks', histtype='stepfilled', lw=2.5, color='grey', alpha=0.4)
plt.xlabel(r'$\langle$'+r'K$_{\rm{s}}$'+r'$\rangle$', fontsize=14)
plt.ylabel('# RRab', fontsize=14)
plt.ylim(0,)
plt.savefig('/Users/apple/Downloads/VVV_phot/RRab/Ks_hist.pdf', format='pdf')
plt.savefig('/Users/apple/Dropbox/Magister/Cursos/InvestigacionB/Latex_RR/Ks_hist.pdf', format='pdf')
plt.show()
In [25]:
#Hacer un catalogo de RRab con magnitudes en J
stilts = 'java -jar -Xmx4096M /Users/apple/Downloads/stilts.jar '
os.system('rm /Users/apple/Downloads/VVV_phot/J/RRab_J.dat')
os.system('rm /Users/apple/Downloads/VVV_phot/RRab/fourier_J.dat')
os.system('rm /Users/apple/Downloads/VVV_phot/RRab/fourier_J_overlap.dat')
J_files = glob.glob('/Users/apple/Downloads/VVV_phot/J/RRab*.dat')
stilts_in = ''
for i,name in enumerate(J_files):
stilts_in = stilts_in + 'in'+str(i+1)+'='+name+' ifmt'+str(i+1)+'=ascii '
os.system(stilts+' tcatn nin='+str(len(J_files))+' '+stilts_in+' ofmt=ascii out=/Users/apple/Downloads/VVV_phot/J/RRab_J.dat');
os.system(stilts+' tmatch2 in1=/Users/apple/Downloads/VVV_phot/J/RRab_J.dat ifmt1=ascii '+
'in2=/Users/apple/Downloads/VVV_phot/RRab/fourier.dat ifmt2=ascii '+
'matcher=exact values1="Num" values2="names" ofmt=ascii out=/Users/apple/Downloads/VVV_phot/RRab/fourier_J.dat');
os.system(stilts+' tmatch2 in1=/Users/apple/Downloads/VVV_phot/J/J_overlap.dat ifmt1=ascii '+
'in2=/Users/apple/Downloads/VVV_phot/RRab/fourier_overlap.dat ifmt2=ascii '+
'matcher=exact values1="Num" values2="names" ofmt=ascii out=/Users/apple/Downloads/VVV_phot/RRab/fourier_J_overlap.dat');
In [26]:
os.chdir('/Users/apple/Downloads/VVV_phot/RRab/')
ra, dec, J, eJ, sep, amp, P, Ks = np.genfromtxt('fourier_J.dat', unpack=True, usecols=(1,2,3,4,5,7,8,9))
num = np.genfromtxt('fourier_J.dat', dtype=None, usecols=(0))
ra_o, dec_o, J_o, eJ_o, sep_o, amp_o, P_o, Ks_o = np.genfromtxt('fourier_J_overlap.dat', unpack=True, usecols=(1,2,3,4,5,7,8,9))
num_o = np.genfromtxt('fourier_J_overlap.dat', dtype=None, usecols=(0))
max_sep = 1
num = num[sep < max_sep]
num_o = num_o[sep_o < max_sep]
ra = ra[sep < max_sep]
ra_o = ra_o[sep_o < max_sep]
dec = dec[sep < max_sep]
dec_o = dec_o[sep_o < max_sep]
J = J[sep < max_sep]
J_o = J_o[sep_o < max_sep]
eJ = eJ[sep < max_sep]
eJ_o = eJ_o[sep_o < max_sep]
amp = amp[sep < max_sep]
amp_o = amp_o[sep_o < max_sep]
P = P[sep < max_sep]
P_o = P_o[sep_o < max_sep]
Ks = Ks[sep < max_sep]
Ks_o = Ks_o[sep_o < max_sep]
In [27]:
def amp_J(amp_Ks):
return 2.6 * ( amp_Ks )**(1.5)
def M_Ks(P):
return - 0.6365 - 2.347 * np.log10( P ) + 0.1747 * log_Z
def M_J(P):
return - 0.2361 - 1.830 * np.log10( P ) + 0.1886 * log_Z
def amp_J_OoII(P):
return 0.064 - 2.481 * np.log10( P ) + 10.345 * ( np.log10( P ) )**3.
def mean_J_mag(Ks):
return 0.98077573567*Ks+0.552882100368
Fe_H = -1.02
log_Z = Fe_H - 1.765
In [45]:
E_JKs = (mean_J_mag(Ks) - Ks) - (M_J(P) - M_Ks(P))
E_JKs_o = (mean_J_mag(Ks_o) - Ks_o) - (M_J(P_o) - M_Ks(P_o))
#plt.figure(figsize=(5,5))
#amlhist(np.concatenate((E_JKs,E_JKs_o)), histtype='stepfilled', lw=2.5, color='grey', alpha=0.4)
#plt.xlabel('E(J-Ks)', size=14)
#plt.ylabel('# RRab', size=14)
#plt.show()
In [29]:
A_Ks = 0.689*E_JKs
A_Ks_o = 0.689*E_JKs_o
Ks0 = Ks - A_Ks
Ks0_o = Ks_o - A_Ks_o
logR = 1 + 0.2*( Ks0 - M_Ks(P) )
logR_o = 1 + 0.2*( Ks0_o - M_Ks(P_o) )
distance = 10**logR
distance_o = 10**logR_o
In [30]:
plt.figure(figsize=(6,6))
amlhist(np.concatenate((distance, distance_o))/1000, bins='blocks', histtype='stepfilled', lw=2.5, color='grey', alpha=0.4);
plt.axvline(8,33, 0, 250, color='red', lw=3, alpha=0.5);
plt.xlabel(r'$d$ (kpc)', size=14);
plt.ylabel('# RRab', fontsize=14)
plt.savefig('/Users/apple/Dropbox/Magister/Cursos/InvestigacionB/Latex_RR/distances.pdf', format='pdf');
In [31]:
def RC_dist(JKs, Ks):
A = 0.73
JKs0 = 0.68
M_Ks = -1.55
rhs = Ks - A*(JKs - JKs0) - M_Ks
return 1/1000.*10**((rhs+5)/5.)
In [32]:
JKs1, Ks1 = np.genfromtxt('../RC/distances/gran01.dat', usecols=(2,3), unpack=True)
JKs2, Ks2 = np.genfromtxt('../RC/distances/gran02.dat', usecols=(2,3), unpack=True)
JKs3, Ks3 = np.genfromtxt('../RC/distances/gran03.dat', usecols=(2,3), unpack=True)
dist1 = RC_dist(JKs1, Ks1)
dist2 = RC_dist(JKs2, Ks2)
dist3 = RC_dist(JKs3, Ks3)
In [33]:
AB = SkyCoord(ra=ra, dec=dec ,unit=(u.degree, u.degree), frame='icrs')
AB_o = SkyCoord(ra=ra_o, dec=dec_o ,unit=(u.degree, u.degree), frame='icrs')
gal_ab = AB.galactic
gal_ab_o = AB_o.galactic
l_ab = np.array([])
l_ab_o = np.array([])
for i,ell in enumerate(gal_ab.l.deg):
if ell > 50.:
l_ab = np.append(l_ab, ell - 360.)
else:
l_ab = np.append(l_ab, ell)
for i,ell in enumerate(gal_ab_o.l.deg):
if ell > 50.:
l_ab_o = np.append(l_ab_o, ell - 360.)
else:
l_ab_o = np.append(l_ab_o, ell)
b_ab = gal_ab.b.deg
b_ab_o = gal_ab_o.b.deg
distance_all = np.concatenate((distance, distance_o))
l_all = np.concatenate((l_ab, l_ab_o))
b_all = np.concatenate((b_ab, b_ab_o))
ra_all = np.concatenate((ra, ra_o))
dec_all = np.concatenate((dec, dec_o))
bins = 32
plt.figure(figsize=(14,5))
plt.subplot(131)
median_1 = np.median(distance_all[l_all > 3.5]/1000)
amlhist(distance_all[l_all > 3.5]/1000, bins=bins, histtype='stepfilled', lw=2.5, color='grey', alpha=0.4);
plt.xlabel(r'$d$ (kpc)', size=16)
plt.ylabel('# RRab', size=16)
plt.title(r'$\ell \geq 3.5^\circ$, '+r'$N_{\rm{RRab}} = %3d $' %len(distance_all[l_all > 3.5]), size=16)
plt.axvline(median_1, 0, 250, color='red', lw=2, alpha=0.5)
plt.xlim(0.,25)
plt.ylim(0,45)
print 'Mediana para l > 3.5:', median_1
###
plt.twinx()
amlhist(dist1, bins=bins, histtype='step', lw=2.5, color='black', alpha=0.8, normed=True);
plt.ylim(0.045, 0.14)
plt.xlim(0.,25)
###
plt.subplot(132)
median_2 = np.median(distance_all[np.abs(l_all) < 3.5]/1000)
amlhist(distance_all[np.abs(l_all) < 3.5]/1000, bins=bins, histtype='stepfilled', lw=2.5, color='grey', alpha=0.4);
plt.xlabel(r'$d$ (kpc)', size=16)
plt.title(r'$3.5^\circ \geq \ell \geq -3.5^\circ$, '+r'$N_{\rm{RRab}} = %3d $' %len(distance_all[np.abs(l_all) < 3.5]), size=16)
plt.axvline(median_2, 0, 250, color='red', lw=2, alpha=0.5)
plt.xlim(0.,25)
plt.ylim(0,45)
print 'Mediana para |l| < 3.5:', median_2
###
plt.twinx()
amlhist(dist2, bins=bins, histtype='step', lw=2.5, color='black', alpha=0.8, normed=True);
plt.ylim(0.045, 0.14)
plt.xlim(0.,25)
###
plt.subplot(133)
median_3 = np.median(distance_all[l_all < -3.5]/1000)
amlhist(distance_all[l_all < -3.5]/1000, bins=bins, histtype='stepfilled', lw=2.5, color='grey', alpha=0.4);
plt.xlabel(r'$d$ (kpc)', size=16)
plt.title(r'$-3.5^\circ \geq \ell$, '+r'$N_{\rm{RRab}} = %3d $' %len(distance_all[l_all < -3.5]), size=16)
plt.axvline(median_2, 0, 250, color='red', lw=2, alpha=0.5)
plt.xlim(0.,25)
plt.ylim(0,45)
print 'Mediana para l < -3.5:', median_3
###
plt.twinx()
amlhist(dist3, bins=bins, histtype='step', lw=2.5, color='black', alpha=0.8, normed=True);
plt.ylim(0.045, 0.14)
plt.xlim(0.,25)
plt.ylabel('# RC stars (normalized)', size=16)
###
plt.tight_layout()
plt.savefig(path+'xshape.pdf', format='pdf')
plt.savefig('/Users/apple/Dropbox/Magister/Cursos/InvestigacionB/Latex_RR/xshape.pdf', format='pdf');
In [58]:
bins = 32
print 'bins of %1.3f kpc \n' % (20/32.)
n_rc_1, _ = np.histogram(dist1, normed=True, bins=bins, range=(0,20))
n_rrl_1, _ = np.histogram(distance_all[l_all > 3.5]/1000, normed=True, bins=bins, range=(0,20))
n_rc_2, _ = np.histogram(dist2, normed=True, bins=bins, range=(0,20))
n_rrl_2, _ = np.histogram(distance_all[np.abs(l_all) < 3.5]/1000, normed=True, bins=bins, range=(0,20))
n_rc_3, _ = np.histogram(dist3, normed=True, bins=bins, range=(0,20))
n_rrl_3, _ = np.histogram(distance_all[l_all < -3.5]/1000, normed=True, bins=bins, range=(0,20))
print 'Normalizado'
print 'KS test RC1 - RRL1: p-value = %1.6f, D = %2.5f, Confidence = %2.1f %%' % (ks(n_rc_1, n_rrl_1)[1], ks(n_rc_1, n_rrl_1)[0], 100.*(1-ks(n_rc_1, n_rrl_1)[1]))
print 'KS test RC2 - RRL2: p-value = %1.6f, D = %2.5f, Confidence = %2.1f %%' % (ks(n_rc_2, n_rrl_2)[1], ks(n_rc_2, n_rrl_2)[0], 100.*(1-ks(n_rc_2, n_rrl_2)[1]))
print 'KS test RC3 - RRL3: p-value = %1.6f, D = %2.5f, Confidence = %2.1f %%' % (ks(n_rc_3, n_rrl_3)[1], ks(n_rc_3, n_rrl_3)[0], 100.*(1-ks(n_rc_3, n_rrl_3)[1]))
print
print 'Sin normalizar'
print 'KS test RC1 - RRL1: p-value = %1.6f, D = %2.5f' % (ks(dist1, distance_all[l_all > 3.5]/1000)[1], ks(dist1, distance_all[l_all > 3.5]/1000)[0])
print 'KS test RC2 - RRL2: p-value = %1.6f, D = %2.5f' % (ks(dist2, distance_all[np.abs(l_all) < 3.5]/1000)[1], ks(dist2, distance_all[np.abs(l_all) < 3.5]/1000)[0])
print 'KS test RC3 - RRL3: p-value = %1.6f, D = %2.5f' % (ks(dist3, distance_all[l_all < -3.5]/1000)[1], ks(dist3, distance_all[l_all < -3.5]/1000)[0])
In [35]:
data1 = dist1.T
np.savetxt('RC1.dat', data1, delimiter='\t', fmt='%2.3f', header='dist (kpc)')
data2 = dist2.T
np.savetxt('RC2.dat', data2, delimiter='\t', fmt='%2.3f', header='dist (kpc)')
data3 = dist3.T
np.savetxt('RC3.dat', data3, delimiter='\t', fmt='%2.3f', header='dist (kpc)')
data4 = (distance_all[l_all > 3.5]/1000).T
np.savetxt('RR1.dat', data4, delimiter='\t', fmt='%2.3f', header='dist (kpc)')
data5 = (distance_all[np.abs(l_all) < 3.5]/1000).T
np.savetxt('RR2.dat', data5, delimiter='\t', fmt='%2.3f', header='dist (kpc)')
data6 = (distance_all[l_all < -3.5]/1000).T
np.savetxt('RR3.dat', data6, delimiter='\t', fmt='%2.3f', header='dist (kpc)')
In [36]:
#KS tests
a = np.array([15.7, 16.1, 15.9, 16.2, 15.9, 16.0, 15.8, 16.1, 16.3, 16.5, 15.5])
b = np.array([15.4, 16.0, 15.6, 15.7, 16.6, 16.3, 16.4, 16.8, 15.2, 16.9, 15.1])
print ks(a,b)
In [ ]:
In [ ]:
In [37]:
from matplotlib.transforms import Affine2D
import mpl_toolkits.axisartist.floating_axes as floating_axes
import mpl_toolkits.axisartist.angle_helper as angle_helper
from matplotlib.projections import PolarAxes
from mpl_toolkits.axisartist.grid_finder import FixedLocator, MaxNLocator, DictFormatter
def setup_axes_l(fig, rect, dist_min, dist_max):
# rotate a bit for better orientation
tr_rotate = Affine2D().translate(0, 0)
# scale degree to radians
tr_scale = Affine2D().scale(np.pi/180., 1.)
tr = tr_rotate + tr_scale + PolarAxes.PolarTransform()
#Number of ticks in l
grid_locator1 = MaxNLocator(10)
tick_formatter1 = angle_helper.FormatterDMS()
#Number of ticks in distance
grid_locator2 = MaxNLocator(25)
grid_helper = floating_axes.GridHelperCurveLinear(tr,
extremes=(13.0, -13.0, dist_min, dist_max),
grid_locator1=grid_locator1,
grid_locator2=grid_locator2,
tick_formatter1=tick_formatter1,
tick_formatter2=None,
)
ax1 = floating_axes.FloatingSubplot(fig, rect, grid_helper=grid_helper)
fig.add_subplot(ax1)
# adjust axis
ax1.axis["right"].set_axis_direction("bottom")
ax1.axis["right"].label.set_axis_direction("bottom")
ax1.axis["right"].label.set_text(r"$d$ [kpc]")
ax1.axis["right"].toggle(ticklabels=True, label=True)
ax1.axis["right"].major_ticklabels.set_axis_direction("bottom")
ax1.axis["right"].label.set_fontsize(15)
ax1.axis["top"].set_axis_direction("bottom")
ax1.axis["top"].label.set_axis_direction("right")
ax1.axis["top"].toggle(ticklabels=True, label=True)
ax1.axis["top"].major_ticklabels.set_axis_direction("right")
ax1.axis["top"].label.set_text(r"$\ell$")
ax1.axis["top"].label.set_fontsize(25)
ax1.axis["left"].set_axis_direction("top")
ax1.axis["left"].toggle(ticklabels=False, label=False)
ax1.axis["bottom"].set_visible(False)
# create a parasite axes whose transData in RA, cz
aux_ax = ax1.get_aux_axes(tr)
aux_ax.patch = ax1.patch # for aux_ax to have a clip path as in ax
ax1.patch.zorder=0.9 # but this has a side effect that the patch is
# drawn twice, and possibly over some other
# artists. So, we decrease the zorder a bit to
# prevent this.
return ax1, aux_ax
In [38]:
fig = plt.figure(1, figsize=(9,5.5))
ax, aux_ax = setup_axes_l(fig, 111, 2.5, 22.5)
#aux_ax.plot(np.linspace(16,-16,100), np.linspace(8.33,8.33,100), color='k', ls='dashed');
aux_ax.scatter(l_ab, distance/1000, marker='*', s=20., color='r', alpha=0.5);
aux_ax.scatter(l_ab_o, distance_o/1000, marker='o', s=20, color='b', alpha=0.8);
ax.grid(True)
plt.tight_layout()
plt.savefig('/Users/apple/Dropbox/Magister/Cursos/InvestigacionB/Latex_RR/cone.pdf', format='pdf');
In [39]:
d_kpc = distance_all/1000
X = np.empty_like(d_kpc)
Y = np.empty_like(d_kpc)
Z = np.empty_like(d_kpc)
for i in range(len(d_kpc)):
gal_coo = SkyCoord(ra=ra_all[i], dec=dec_all[i], distance=d_kpc[i], unit=(u.degree, u.degree, u.kpc))
X[i] = gal_coo.galactocentric.x.value
Y[i] = gal_coo.galactocentric.y.value
Z[i] = gal_coo.galactocentric.z.value
import plotly.plotly as py from plotly.graph_objs import * py.sign_in('fgran','vyvc5yqufy')
GC = Scatter3d(x=[0], y=[0], z=[0], name='Galactic Center', mode='markers', marker=Marker(color='rgb(255,0,0)', symbol='circle', line=Line(width=2)), visible=True) RRab = Scatter3d(x=X, y=Y, z=Z, name='RRab outer bulge', mode='markers', marker=Marker(color='rgb(104,104,104)', symbol='.', size=3, line=Line(width=1)), visible=True)
data = Data([GC, RRab]) xaxis = XAxis(title='X', autorange=True, range=[-6,15], domain=[-6,16], nticks=15) yaxis = YAxis(title='Y', autorange=False, range=[-4,4], domain=[-4,4], nticks=5) zaxis = ZAxis(title='Z', autorange=True, range=[-5.5,0.5], domain=[0.5,5.5], nticks=5) layout = Layout(autosize=False, width=600, height=500, title = 'VVV RRab stars in the Outer Bulge', scene=Scene(xaxis=xaxis, yaxis=yaxis, zaxis=zaxis), margin=Margin(l=50, r=50, b=50, t=50, pad=4)) figure = Figure(data=data, layout=layout) py.iplot(figure, filename='RRab_outer_bulge')
In [40]:
J_b201, Ks_b201 = np.genfromtxt('/Users/apple/Dropbox/Respaldo/vars/colors/b201_allcol.dat', usecols=(12,18), unpack=True)
In [41]:
#CMD
plt.figure(figsize=(8,8))
plt.scatter(J_b201-Ks_b201, Ks_b201, marker='.', color='gray', alpha=0.3, s=5, rasterized=True)
plt.text(0.58, 12.5, 'Bulge RGB', rotation=80, size=14.5)
plt.text(0.3, 12.5, 'Disk MS', rotation=100, size=14.5)
plt.scatter(J-Ks, Ks, marker='*', color='r', alpha=0.5, s=20, rasterized=True)
plt.scatter(J_o-Ks_o, Ks_o, marker='o', color='b', alpha=0.8, s=20, rasterized=True)
plt.xlim(-0.5, 1.5)
plt.ylim(18.5, 11)
plt.xlabel(r'$J-K_{\rm s}$', size=16)
plt.ylabel(r'$K_{\rm s}$', size=16)
plt.savefig('/Users/apple/Dropbox/Magister/Cursos/InvestigacionB/Latex_RR/cmd.pdf', format='pdf', dpi=300)
In [42]:
data = np.vstack((ra_all, dec_all, l_all, b_all, X, Y, Z, np.concatenate((J, J_o)),
np.concatenate((Ks, Ks_o)), np.concatenate((P,P_o)),
np.concatenate((amp, amp_o)), distance_all/1000)).T
In [43]:
np.savetxt('VVV_RRL_Outer_Bulge.dat', data, delimiter='\t',
fmt='%3.5f %2.5f %2.5f %2.5f %2.2f %2.2f %2.2f %2.3f %2.3f %1.6f %1.2f %2.2f',
header='RA Dec L B X Y Z J Ks P Amp dist(kpc)')
In [44]:
#np.savetxt(sys.stdout, data, newline=' \\\\\n', delimiter='&',
# fmt='%3.5f %2.5f %2.5f %2.5f %2.2f %2.2f %2.2f %2.3f %2.3f %1.6f %1.2f %2.2f',
# header='RA Dec L B X Y Z J Ks P Amp dist(kpc)')
In [ ]: