In [87]:
%matplotlib inline
import numpy as np
import sep
from astropy.io import fits
import matplotlib.pyplot as plt
from scipy import interpolate
from sorting_routines import *
In [89]:
fname_001 = "sim_cube_F090W_487_001.slp.fits"
fits.info(fname_001)
hdulist = fits.open(fname_001)
data_001 = hdulist[0].data[0]
data_001 = data_001.byteswap().newbyteorder()
print np.size(data_001)
for i in range(2048):
for j in range(2048):
if(np.isnan(data_001[i,j])):
data_001[i,j] = 0
In [ ]:
bkg = sep.Background(data_001)
In [90]:
#print hdulist[0].header
m_ab_zeropoint = hdulist[0].header['ABMAG']
print m_ab_zeropoint
In [62]:
plt.imshow(np.log10( data_001 - data_001.min() + 1.0e-3*data_001.max() ),origin="lower")
Out[62]:
In [91]:
plt.imshow(bkg,origin="lower")
Out[91]:
In [92]:
#subtract background
data_sub = data_001 - bkg
In [93]:
plt.imshow(np.log10( data_sub - data_sub.min() + 1.0e-3*data_sub.max() ),origin="lower")
Out[93]:
In [94]:
objects = sep.extract(data_sub, 3.0, err=bkg.globalrms)
objects_w_bkg = sep.extract(data_001, 3.0, err=bkg.globalrms)
In [95]:
print len(objects)
In [96]:
pixel_radius = 35
print pixel_radius
flux, fluxerr, flag = sep.sum_circle(data_sub, objects['x'], objects['y'], pixel_radius, err=bkg.globalrms, gain=1.0)
In [70]:
flux_w_bkg, fluxerr_w_bkg, flag_w_bkg = sep.sum_circle(data_001, objects['x'], objects['y'], pixel_radius, err=bkg.globalrms, gain=1.0)
In [97]:
print flux
In [98]:
print flux_w_bkg
In [99]:
mag_sep = -2.5*np.log10(flux) + m_ab_zeropoint
idx = sorted_index(-1.*mag_sep)
mag_sep = mag_sep[idx]
print mag_sep
In [100]:
fname = "stars_f090w.cat"
fp = open(fname,"r")
fl = fp.readlines()
fp.close()
In [102]:
nexp = 4
nobj = len(fl)/nexp
j = 0
mag_cat = np.zeros([nexp,nobj])
f_real_cat = np.zeros([nexp,nobj])
f_true_cat = np.zeros([nexp,nobj])
for k in range(nexp):
for i in range(nobj):
mag_cat[k,i] = fl[j].split()[5]
f_real_cat[k,i] = fl[j].split()[8]
f_true_cat[k,i] = fl[j].split()[7]
j += 1
print mag_cat[0,:]
In [104]:
print mag_sep - mag_cat[0,:]
In [105]:
#dmag_aperture_correction = 0.204462939361 ## 0.3"
dmag_aperture_correction = 0.0274911626418 ## 35 pixel aperture correction
In [106]:
k = 0
print mag_sep - mag_cat[k,:] - dmag_aperture_correction
In [107]:
k = 0
dmag_real = -2.5*np.log10(f_real_cat[k,:]/f_true_cat[k,:])
print dmag_real
In [108]:
k = 0
dmag_real = -2.5*np.log10(f_real_cat[k,:]/f_true_cat[k,:])
mag_error = mag_sep - mag_cat[k,:] - dmag_aperture_correction - dmag_real
print mag_error
In [109]:
plt.scatter(mag_cat[0,:],mag_error)
plt.xlabel(r'$M [AB]$')
plt.ylabel(r'$\Delta M [AB]$')
print np.mean(mag_error)
print np.std(mag_error)
print 10**(-0.4*(np.mean(mag_error)+np.std(mag_error))) - 1.
print 10**(-0.4*(np.mean(mag_error))) - 1.
print 10**(-0.4*(np.mean(mag_error)-np.std(mag_error))) - 1.
In [113]:
mag_sep_w_bkg = -2.5*np.log10(flux_w_bkg) + m_ab_zeropoint
idx = sorted_index(-1.*mag_sep_w_bkg)
mag_sep_w_bkg = mag_sep_w_bkg[idx]
print mag_sep_w_bkg
In [114]:
k = 0
dmag_real = -2.5*np.log10(f_real_cat[k,:]/f_true_cat[k,:])
mag_error_w_bkg = mag_sep_w_bkg - mag_cat[k,:] - dmag_aperture_correction - dmag_real
In [115]:
plt.scatter(mag_cat[0,:],mag_error_w_bkg)
plt.xlabel(r'$M [AB]$')
plt.ylabel(r'$\Delta M [AB]$')
Out[115]:
In [ ]: