In [1]:

    

%pylab inline
from simqso.sqgrids import *
from simqso import sqbase
from simqso.sqmodels import QLF_McGreer_2013


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

# set up a luminosity-redshift grid
M = AbsMagVar(UniformSampler(-30,-25),restWave=1450)
z = RedshiftVar(UniformSampler(1,5))
MzGrid = QsoSimGrid([M,z],(4,3),2,seed=12345)
scatter(MzGrid.z,MzGrid.absMag,)
xlabel('z')
ylabel('absolute mag');


    

In [3]:

    

# set up a flux-redshift grid with Lya EW as a third dimension
m = AppMagVar(UniformSampler(17,22),'SDSS-i')
lya = GaussianLineEqWidthVar(UniformSampler(0,2),'LyaEW',1215.67,10.)
mzlyaGrid = QsoSimGrid([m,z,lya],(4,3,2),2,seed=12345)
scatter(mzlyaGrid.z,mzlyaGrid.appMag,c=mzlyaGrid.LyaEW)
cb = colorbar()
xlabel('z')
ylabel('apparent mag')
cb.set_label('log EW(Lya)')


    

In [4]:

    

# Generate m,z points by sampling from a model QLF
cosmo = QLF_McGreer_2013.cosmo
kcorr = sqbase.ContinuumKCorr('SDSS-i',1450)
qlfGrid =  generateQlfPoints(QLF_McGreer_2013,
                             (17,22),(4.5,5.3),kcorr,
                             skyArea=100,
                             qlfseed=12345,gridseed=67890)


    

    




integration returned  43  objects

In [5]:

    

scatter(qlfGrid.z,qlfGrid.appMag)
xlabel('z')
ylabel('apparent mag')


    

    
Out[5]:





Text(0,0.5,'apparent mag')






    

In [6]:

    

# test the canonical values for power law continuum slopes in FUV/NUV
contVar = BrokenPowerLawContinuumVar([GaussianSampler(-1.5,0.3),
                                      GaussianSampler(-0.5,0.3)],[1215.],
                                     seed=12345)
contVals = contVar(1000)
_ = hist(contVals[:,0],alpha=0.8)
_ = hist(contVals[:,1],alpha=0.8)


    

In [7]:

    

# Add gaussian emission lines
g = GaussianSampler
lya = [g(1215.7,0.5),g(80,20),g(7,2)]
civ = [g(1550.0,1.5),g(20,10),g(10,2)]
emLinesVar = GaussianEmissionLinesTemplateVar([lya,civ],seed=12345)
emLinesVals = emLinesVar(1000)
figure(figsize=(14,4))
for k in range(3):
    subplot(1,3,k+1)
    for j in range(2):
        hist(emLinesVals[:,j,k],alpha=0.8)


    

In [8]:

    

# add gaussian lines to the qlf grid
qlfGrid.addVar(emLinesVar)
print(qlfGrid)


    

    




 absMag   appMag    z      emLines [2,3]   
-------- ------- ------- ------------------
-24.7238 21.9073  5.1932 1216.23 .. 13.1207
-24.6332 21.7689 4.65953 1215.99 .. 13.4276
-24.6492 21.7199 4.58749  1214.51 .. 14.337
-24.9536 21.4204 4.59819 1215.31 .. 12.2846
-25.0347 21.4418 4.82577 1216.15 .. 10.2813
-25.5588 20.9271 4.84734 1216.23 .. 7.64866
-24.8302 21.8213 5.24391 1215.87 .. 10.3208
-25.1854  21.321 4.89462  1215.32 .. 9.9414
-24.7136 21.8302 4.98235 1216.81 .. 8.40218
-25.7508 20.7557 4.89495   1215.5 .. 9.6769
     ...     ...     ...                ...
  -24.85 21.7669 5.15819 1216.36 .. 12.7096
-24.5976 21.7357  4.5111 1215.95 .. 6.94498
-25.7891  20.663 4.77046 1215.91 .. 11.5107
-25.0814 21.3814 4.79489 1215.69 .. 11.9493
-26.4142 20.0719 4.84799 1215.88 .. 8.00809
-26.5823 19.7568 4.52347 1216.76 .. 11.0919
-25.2844 21.3279 5.14694 1215.16 .. 13.0083
-24.8029 21.7325 4.96231 1215.64 .. 7.98605
-25.6632 20.9112 5.05508 1215.34 .. 6.04667
-24.8754 21.5792 4.77637 1215.18 .. 11.1137
-25.1863 21.3843 5.04601 1214.93 .. 8.78715
Length = 43 rows

In [9]:

    

scatter(qlfGrid.z,qlfGrid.appMag,c=qlfGrid.emLines[:,0,1])
cb = colorbar()
xlabel('z')
ylabel('apparent mag')
cb.set_label('EW(Lya)')


    

In [10]:

    

# Now use the BOSS DR9 emission line template, including the Baldwin Effect
lineTemplate = generateBEffEmissionLines(qlfGrid.absMag,seed=12345)
lines = lineTemplate()
lines.shape


    

    
Out[10]:





(43, 62, 3)

In [11]:

    

scatter(qlfGrid.z,qlfGrid.appMag,c=lines[:,13,1])
cb = colorbar()
xlabel('z')
ylabel('apparent mag')
cb.set_label('EW(CIV)')


    

In [12]:

    

vdblines = generateVdBCompositeEmLines(minEW=0)
vdblines = vdblines(qlfGrid.nObj)


    

    




using the following lines from VdB template:  Ly{epsilon},CIII,Ly{beta},ArI,FeIII,CIII*,Ly{alpha},NV,SiII,OI,CII,SiIV,CIV,HeII,OIII],NIV,NIII],FeII,SiII,AlIII,SiIII],CIII],FeIII,FeIII,FeII,FeII,FeII,CII],[NeIV],[OII],FeII,AlII],MgII,FeII,OIII,HeI,FeI,FeII,[NeV],[NeV],FeII,[FeVII],[OII],[FeVII],FeII,FeII,[NeIII],HeI,[NeIII],[FeV],H{delta},FeII,FeII,[FeII],[FeII],H{gamma},[OIII],FeII,FeII,HeII,H{beta},FeII,[OIII],[OIII],[FeVII],[FeVI],[NI],[FeVII],FeII,[FeXIV],[ClIII],[FeVII],HeI,[FeVII],[OI],[OI],[NII],H{alpha},[NII],[SII],[SII],HeI,[ArIII],[OII],[NiIII]

In [13]:

    

wave = arange(3000,1e4,0.5)
z = 2.0
def quick_spec(emlines):
    lineWave,eqWidth,sigma = emlines.T * (1+z)
    spec = np.zeros_like(wave)
    A = eqWidth/(np.sqrt(2*np.pi)*sigma)
    twosig2 = 2*sigma**2
    nsig = np.sqrt(-2*np.log(1e-3/A))
    nsig = (nsig*np.array([[-1.],[1]])).T
    for i in range(lineWave.shape[0]):
        i1,i2 = np.searchsorted(wave,lineWave[i]+nsig[i]*sigma[i])
        if i2 != i1:
            lineprofile = A[i]*np.exp(-(wave[i1:i2]-lineWave[i])**2
                                       / twosig2[i])
            spec[i1:i2] += lineprofile
    return spec


    

In [14]:

    

figure(figsize=(12,4))
plot(wave/(1+z),quick_spec(lines[0]),label='BOSS DR9')
plot(wave/(1+z),quick_spec(vdblines[0]),label='Vanden Berk composite')
legend()


    

    
Out[14]:





<matplotlib.legend.Legend at 0x1c1804a860>






    

In [ ]: