In [1]:
%pylab inline
from astropy.cosmology import Planck13
from simqso.sqgrids import *
from simqso import sqbase
from simqso.sqrun import buildSpectraBulk,buildQsoSpectrum
random.seed(12345)
In [2]:
# cover 1000A to 20um at R=1000
wave = sqbase.fixed_R_dispersion(1000,20e4,1000)
In [3]:
# just make up a few magnitude and redshift points
nqso = 5
M = AbsMagVar(FixedSampler(linspace(-27,-25,nqso)[::-1]),restWave=1450)
z = RedshiftVar(FixedSampler(linspace(2,4,nqso)))
qsos = QsoSimPoints([M,z],cosmo=Planck13,units='luminosity')
scatter(qsos.z,qsos.absMag)
xlabel('z')
ylabel('$M_{1450}$')
ylim(-24.75,-27.25);
In [4]:
# use the canonical values for power law continuum slopes in FUV/NUV, with breakpoint at 1215A
contVar = BrokenPowerLawContinuumVar([GaussianSampler(-1.5,0.3),
GaussianSampler(-0.5,0.3)],
[1215.])
In [5]:
# add two dust components as in Lyu+Rieke 2017, but reduce the hot dust flux by a factor of 2
subDustVar = DustBlackbodyVar([ConstSampler(0.05),ConstSampler(1800.)],
name='sublimdust')
subDustVar.set_associated_var(contVar)
hotDustVar = DustBlackbodyVar([ConstSampler(0.1),ConstSampler(880.)],
name='hotdust')
hotDustVar.set_associated_var(contVar)
In [6]:
# generate lines using the Baldwin Effect emission line model from BOSS DR9
emLineVar = generateBEffEmissionLines(qsos.absMag)
In [7]:
# the default iron template from Vestergaard & Wilkes 2001 was modified to fit BOSS spectra
fescales = [(0,1540,0.5),(1540,1680,2.0),(1680,1868,1.6),(1868,2140,1.0),(2140,3500,1.0)]
feVar = FeTemplateVar(VW01FeTemplateGrid(qsos.z,wave,scales=fescales))
In [8]:
# Now add the features to the QSO grid
qsos.addVars([contVar,subDustVar,hotDustVar,emLineVar,feVar])
In [9]:
# ready to generate spectra
_,spectra = buildSpectraBulk(wave,qsos,saveSpectra=True)
In [10]:
qsos.data
Out[10]:
In [11]:
figure(figsize=(14,6))
for sp in spectra:
plot(wave,sp)
yscale('log')
In [12]:
figure(figsize=(14,8))
for i,sp in enumerate(spectra):
plot(wave,sp+i*2e-16)
xlim(3000,7000)
Out[12]:
In [13]:
# compare the resulting SEDs to the mean SED from Lyu+Rieke 2017
from astropy.table import Table
try:
lr17 = Table.read('apjaa7051t4_mrt.txt',format='ascii')
lr17wave = array(10**lr17['lambda'])
lamflam = array(spectra)*wave
figure(figsize=(8,4))
for i,z in enumerate(qsos.data['z']):
j = searchsorted(wave/(1+z),5100)
fscl = 1/lamflam[i,j]
plot(wave/(1+z)/1e4,log10(fscl*lamflam[i]))
j = searchsorted(lr17wave,0.51)
fscl = lr17['LogL-Norm'][j]
plot(lr17wave,lr17['LogL-Norm']-fscl,c='k')
axvline(0.51,c='gray',ls='--')
xlim(0.1,5)
ylim(-1,1)
xscale('log');
except IOError:
pass
Divide the spectrum into components. Offset the emission line and iron template features for clarity.
In [14]:
sp,comp = buildQsoSpectrum(wave,qsos.cosmo,qsos.getVars(SpectralFeatureVar),
qsos.data[0],save_components=True)
z1 = sp.z + 1
figure(figsize=(9,5))
plot(sp.wave/z1,sp.wave*sp.f_lambda,label='total')
for cname,cspec in comp.items():
flam = cspec.f_lambda
if cname in ['emLines','fetempl']:
flam = 0.25*comp['slopes'].f_lambda*(1+flam)
plot(cspec.wave/z1,cspec.wave*flam,label=cname)
legend(ncol=3)
ylim(1e-14,2e-12)
xscale('log')
yscale('log')
xlabel('$\lambda [\AA]$')
ylabel(r'$\nu{f}_\nu$');
In [ ]: