超新星のSEDから光度曲線を計算


計算式については「天体のfluxからBroad-band magnitudeを計算」を参照。
 

超新星のtemplate SED


Ia型はEric HsiaoのHPからダウンロード可能。
Peter NugentのHPからも、Type Ia 1991T-like, Type Ia 1991bg-like, Type Ib/c, Type Ib/c high-velocity, Type IIP, Type IIL, Type IInなどの各種templateが手に入る。
 

Filter Response


ここでは、Suprime-Cam Filtersをここからダウンロードして使っています。
# 使っているレスポンスは「フィルター + CCD + 主鏡反射率 + 主焦点補正光学系 + 大気透過率(secz=1.2) 」のもの。
 

計算手順


以下、HsiaoのType Ia template(snflux_1a.dat)をもとに光度曲線を計算する。

  1. まずHsiaoのtemplateでは-20日から85日までのSEDが一つのファイルにまとまっているので、あらかじめspectrum_m20d.dat, spectrum_m19d.dat,…, spectrum_m01d.dat, spectrum_p00d.dat, spectrum_p01d.dat,…, spectrum_p85d.datのように各日毎にファイルを分離しておく。
  2. 後はABMAGの計算式、
    \displaystyle  {\rm ABMAG}(R)=-2.5\log_{10}\frac{\int f_{\nu}\frac{R}{\nu}d\nu}{\int \frac{R}{\nu}d\nu}-48.6
    に従って計算するだけ。ただし、z>0の計算をする場合はスペクトルと時間が波長方向に(1+z)倍に伸びるので注意すること。

その他詳しい計算などは
K‐Corrections and Extinction Corrections for Type Ia Supernovae
Nugent et al. 2002 (ADS)
を参照。
 

コード例


distance modulusの計算にはPythonの宇宙論パッケージCosmoloPyを使用。
ここではVバンド絶対等級で-19.2のIa型超新星がSubaru/SCのB-,V-,Rc-,i’-,z’-bandでどのように観測されるかを計算している。

import numpy as np
import scipy.interpolate
import cosmolopy.distance as cosmocd

# convert wavelength to frequency
def convert_nu(filter):
    filter_lam=filter[:,0]*1.0e-8 # [cm]
    filter_nu=np.array([0.0] * len(filter_lam))
    for i in range(len(filter_lam)):
        filter_nu[i]=c/filter_lam[len(filter_lam)-1-i] # [Hz]
    return filter_nu*1.0e-15

# convert f_lambda to f_nu
def convert_fnu(filter):
    filter_lam=filter[:,0]*1.0e-8 # [cm]
    filter_responce=np.array([0.0] * len(filter_lam))
    for i in range(len(filter_lam)):
        filter_responce[i]=filter[len(filter_lam)-1-i,1]
    return filter_responce

# calculate ABmag
def abmagnitude(model_nu15,model_fnu,filter_nu15,filter_responce,z):
    flux=scipy.interpolate.interp1d(model_nu15,model_fnu,kind="linear")
    numerator=0.0
    denominator=0.0
    for i in range(len(filter_nu15)-1):
        for j in range(len(model_fnu)):
            if (filter_nu15[i]*(1.+z) = model_nu15[j-1]):
                dnu=(filter_nu15[i+1]-filter_nu15[i])
                numerator+=flux(filter_nu15[i]*(1.+z))*(1.+z)*filter_responce[i]/filter_nu15[i]*dnu
                denominator+=filter_responce[i]/filter_nu15[i]*dnu
    mag=-2.5*np.log10(numerator/denominator)-48.6
    return mag

# cosmological parameters
cosmo = {'omega_M_0':0.3, 'omega_lambda_0':0.7, 'omega_k_0':0.0, 'h':0.70}
# speed of light [cm s-1]
c = 2.989e10
# absolute V-band magnitude
MaV=-19.2

# load filter information
Bfilter=np.loadtxt('B-bandレスポンスファイル')
Vfilter=np.loadtxt('V-bandレスポンスファイル')
Rfilter=np.loadtxt('R-bandレスポンスファイル')
ifilter=np.loadtxt('i-bandレスポンスファイル')
zfilter=np.loadtxt('z-bandレスポンスファイル')

# define responce (as a function of nu)
Bfilter_nu15=convert_nu(Bfilter) # [E-15 Hz]
Bfilter_responce=convert_fnu(Bfilter)
Vfilter_nu15=convert_nu(Vfilter) # [E-15 Hz]
Vfilter_responce=convert_fnu(Vfilter)
Rfilter_nu15=convert_nu(Rfilter) # [E-15 Hz]
Rfilter_responce=convert_fnu(Rfilter)
ifilter_nu15=convert_nu(ifilter) # [E-15 Hz]
ifilter_responce=convert_fnu(ifilter)
zfilter_nu15=convert_nu(zfilter) # [E-15 Hz]
zfilter_responce=convert_fnu(zfilter)

# calculate AB magnitude from z=0 to z=1.2 in dz=0.05 interval
for z in np.arange(0.0,1.21,0.05):
    fout=open('z%03d' % (z*100.0),'w')
    if z>0:
        comoving_distance=cosmocd.comoving_distance(z, **cosmo) # comoving distance [Mpc]
        distance_modulus=5.*np.log10(comoving_distance)+5.*np.log10(1.+z)+25.
    else:
        comoving_distance=0.0
        distance_modulus=0.0
    # epoch loop
    for epoch in np.arange(-20,86,1):
        if epoch < 0:
            specdata='spectrum_m%02d.dat' % -epoch
        else:
            specdata='spectrum_p%02d.dat' % epoch
        model_lam=np.loadtxt(specdata)[:,0]*1.0e-8 # [cm]
        model_flam=np.loadtxt(specdata)[:,1]
        NBIN=len(model_lam)
        model_fnu= np.array([0.0] * NBIN)
        model_nu= np.array([0.0] * NBIN)
        for i in range(NBIN):
            model_fnu[i]=model_lam[NBIN-1-i]*model_lam[NBIN-1-i]*(model_flam[NBIN-1-i]*1.0e8)/c
            model_nu[i]=c/model_lam[NBIN-1-i]
        model_nu15=model_nu*1.0e-15
        # integration
        Bmag=abmagnitude(model_nu15,model_fnu,Bfilter_nu15,Bfilter_responce,z)+MaV+distance_modulus
        Vmag=abmagnitude(model_nu15,model_fnu,Vfilter_nu15,Vfilter_responce,z)+MaV+distance_modulus
        Rmag=abmagnitude(model_nu15,model_fnu,Rfilter_nu15,Rfilter_responce,z)+MaV+distance_modulus
        imag=abmagnitude(model_nu15,model_fnu,ifilter_nu15,ifilter_responce,z)+MaV+distance_modulus
        zmag=abmagnitude(model_nu15,model_fnu,zfilter_nu15,zfilter_responce,z)+MaV+distance_modulus
        print '%4.3f %5.1f %6.3f %6.3f %6.3f %6.3f %6.3f' % (z,epoch*(1.+z),Bmag,Vmag,Rmag,imag,zmag)

コメントを残す