In [3]:
cd ~/Dropbox/notebooks/ghost_halos/
In [4]:
import numpy as np
import matplotlib.pyplot as plt
import radprofile as rp
from astropy.io import fits
from scipy.integrate import trapz
%matplotlib inline
In [5]:
## Open function to load the profile object
def open_profile(fname):
result = rp.Profile()
ffile = fits.open(fname)
pdata = ffile[1].data
result.rleft = pdata['R'][:,0]
result.rright = pdata['R'][:,1]
try:
result.surbri = pdata['SUR_FLUX'] # photons/cm**2/pixel**2/s
result.surbri_err = pdata['SUR_FLUX_ERR']
except:
print "Could not find value:", value
print "Your options are:"
print pdata.columns.names
return result
In [6]:
SFLUX = 0.0007349181010116529 # 0.8 - 1.2 keV flux [photons cm^-2 s^-1]
fprofile = open_profile('herx1_c67_1keV_radprofile_flux.fits')
ftotflux = np.sum(fprofile.surbri * fprofile.area)
print "Total flux in profile:", ftotflux
print "This is " + np.str((ftotflux/SFLUX)*100.0)[0:4] + "% of the HETG measured flux"
In [7]:
plt.errorbar(fprofile.rmid, fprofile.surbri, yerr=fprofile.surbri_err, \
capsize=0, lw=2, color='k', marker='', ls=''),
plt.loglog()
plt.xlabel('Radius [pixels]')
plt.ylabel('Surface brightness [flux pixels$^{-2}$]')
Out[7]:
In [8]:
## Convert the profile from pixel units to arcsec
arcsec_per_pixel = 0.5 # arcsec/pixel
rmid_arcsec = fprofile.rmid * arcsec_per_pixel # arcsec
sbri_arcsec = fprofile.surbri / arcsec_per_pixel**2 # flux per arcsec^2
sbri_err_arcsec = fprofile.surbri_err / arcsec_per_pixel**2
plt.errorbar(rmid_arcsec, sbri_arcsec, yerr=sbri_err_arcsec, \
capsize=0, lw=2, color='k', marker='', ls='')
plt.loglog()
plt.xlabel('Observation Angle [arcsec]')
plt.ylabel('Surface brightness [flux arcsec$^{-2}$]')
Out[8]:
In [9]:
from astropy.modeling import models
from scipy.optimize import leastsq
In [21]:
XTRANS = 0.5 # pixels, core-wing transition radius, min x val for wing portion
def wing_model(x, amp, p, xtrans=XTRANS):
wmod = models.PowerLaw1D(amplitude=amp*SFLUX, x_0=1.0, alpha=p)
expmod = np.exp(-xtrans/x)
return wmod(x) * expmod
def core_model(x, amp, p, g):
cmod = models.Beta1D(amplitude=amp*SFLUX, x_0=0.0, alpha=p, gamma=g)
return cmod(x)
def psf_model(params, x):
a1, p1, g1, a2, p2, c0 = params
bkg = models.Const1D(c0)
result = core_model(x, a1, p1, g1) + wing_model(x, a2, p2) + bkg(x)
return result # arcsec^-2
In [22]:
def chi(params, x, y, err):
return (y - psf_model(params, x))/err
In [23]:
def plot_fit(params, x, y, err, comps=False):
plt.errorbar(x, y, yerr=err, \
capsize=0, lw=2, color='k', marker='o', markersize=3.0, ls='', label='')
plt.loglog()
rsmooth = np.logspace(-1,np.log10(max(x)),100)
plt.plot(rsmooth, psf_model(params, rsmooth), \
lw=2, alpha=0.5, color='r')
if comps:
a1, p1, g1, a2, p2, c0 = params
core = core_model(rsmooth, a1, p1, g1)
wing = wing_model(rsmooth, a2, p2)
bkg = models.Const1D(amplitude=c0)
plt.plot(rsmooth, core, 'k:', label='Core')
plt.plot(rsmooth, wing, 'k--', label='Wings')
plt.plot(rsmooth, bkg(rsmooth), '-', color='0.5', label='Background')
return
In [24]:
## Pick a good value for XTRANS
atest, ptest, xtest = 1.e-10, 2.0, 0.5
rtest = np.logspace(-1.0,1.0,100)
PLonly = models.PowerLaw1D(amplitude=atest*SFLUX, x_0=1.0, alpha=ptest)
wtest = wing_model(rtest, atest, ptest, xtrans=xtest)
ax = plt.subplot(211)
ax.plot(rtest, wtest, 'k--', label='wing model')
ax.plot(rtest, PLonly(rtest), 'r-', label='power law')
ax.axvline(3.0, ls=':', color='k')
ax.xaxis.set_ticklabels('')
plt.loglog()
plt.legend(loc='lower left', frameon=False)
ax2 = plt.subplot(212)
ax2.plot(rtest, (PLonly(rtest)-wtest)/PLonly(rtest) * 100, 'k-')
plt.loglog()
plt.axhline(10.0, ls=':', color='r')
plt.axvline(3.0, ls=':', color='k')
ax2.set_xlabel('Radius')
ax2.set_ylabel('Percent deviation')
Out[24]:
In [25]:
pinit = [0.1, 2.0, 1.0, 1.e-3, 2.0, 1.e-9]
plot_fit(pinit, rmid_arcsec, sbri_arcsec, sbri_err_arcsec, comps=True)
plt.ylim(1.e-10, 1.e-3)
Out[25]:
In [26]:
(pfit, temp) = leastsq(chi, pinit, args=(rmid_arcsec, sbri_arcsec, sbri_err_arcsec))
print pfit
In [27]:
plot_fit(pfit, rmid_arcsec, sbri_arcsec, sbri_err_arcsec, comps=True)
plt.ylim(1.e-10, 1.e-3)
plt.xlim(0.1,300)
plt.legend(loc='upper right', frameon=False)
plt.xlabel('Observation angle [arcsec]')
plt.ylabel('Surface brightness [photon flux arcsec$^{-2}$]')
print 'Reduced Chi-squared =', \
np.sum(chi(pfit,rmid_arcsec,sbri_arcsec,sbri_err_arcsec)**2) / len(pfit)
In [28]:
from astropy.io import ascii
from astropy.table import Table
In [29]:
sfilename = 'plot_radprofile_fitparams.txt'
par_table = Table([pfit], names=['params'])
ascii.write(par_table, sfilename)
In [30]:
print pfit
In [31]:
## Fraction of flux accounted for by background
bkg_totflux = pfit[-1] * np.pi * (fprofile.rright[-1]*arcsec_per_pixel)**2
print "Background contributes " + \
np.str(bkg_totflux/ftotflux * 100.0)[0:4] + \
"% of profile flux"
print "This happens to be " + np.str(bkg_totflux/SFLUX * 100.0)[0:4] + "% of the source flux"
In [32]:
xvals = np.logspace(-1,3.0,200)
yvals = core_model(xvals, pfit[0], pfit[1], pfit[2]) + \
wing_model(xvals, pfit[3], pfit[4])
int_flux = trapz( yvals*2.0*np.pi*xvals, xvals)
print "Total flux in PSF model:", int_flux
print "This is " + np.str(int_flux/SFLUX * 100.0)[0:4] + "% of source flux"
print "This is " + np.str(int_flux/ftotflux * 100.0)[0:4] + "% of the profile flux"
Source flux (measured from HETG) should be source minus bkg, so only 45.5% of 0th order point source light is accounted for -- note that pileup is an issue for this observation
50.2 + 42.6 = 92.8% of the profile flux is accounted for with my PSF model
In [ ]: