Skip to content
GitLab
Commit cf67d264 authored by Zhang Xin's avatar Zhang Xin
Browse files

init

parents
Show whitespace changes
Inline Side-by-side
Showing
with 565 additions and 0 deletions
+565 -0
produceSED_1.py 0 → 100644
+ 565
- 0
View file @ cf67d264
'''
Author: xin zhangxinbjfu@gmail.com
Date: 2021-08-31 14:58:40
LastEditors: xin zhangxinbjfu@gmail.com
LastEditTime: 2022-11-21 16:09:25
FilePath: /src/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/produceSED_1.py
Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE
'''
from tkinter.font import names
from pylab import *
import h5py as h5
from astropy.table import Table
from scipy import interpolate
import astropy.constants as cons
from scipy.interpolate import InterpolatedUnivariateSpline, UnivariateSpline, interp1d
import healpy as hp
from datatable import dt,f
import numpy as np
from astropy.cosmology import FlatLambdaCDM
import os
import math
def lyman_forest(wavelen, flux, z):
"""
Compute the Lyman forest mean absorption of an input spectrum,
according to D_A and D_B evolution from Madau (1995).
The waveeln and flux are in observed frame
"""
if z<=0:
flux0 = flux
else:
nw = 200
istep = np.linspace(0,nw-1,nw)
w1a, w2a = 1050.0*(1.0+z), 1170.0*(1.0+z)
w1b, w2b = 920.0*(1.0+z), 1015.0*(1.0+z)
wstepa = (w2a-w1a)/float(nw)
wstepb = (w2b-w1b)/float(nw)
wtempa = w1a + istep*wstepa
ptaua = np.exp(-3.6e-03*(wtempa/1216.0)**3.46)
wtempb = w1b + istep*wstepb
ptaub = np.exp(-1.7e-3*(wtempb/1026.0)**3.46\
-1.2e-3*(wtempb/972.50)**3.46\
-9.3e-4*(wtempb/950.00)**3.46)
da = (1.0/(120.0*(1.0+z)))*np.trapz(ptaua, wtempa)
db = (1.0/(95.0*(1.0+z)))*np.trapz(ptaub, wtempb)
if da>1.0: da=1.0
if db>1.0: db=1.0
if da<0.0: da=0.0
if db<0.0: db=0.0
flux0 = flux.copy()
id0 = wavelen<=1026.0*(1.0+z)
id1 = np.logical_and(wavelen<1216.0*(1.0+z),wavelen>=1026.0*(1.0+z))
flux0[id0] = db*flux[id0]
flux0[id1] = da*flux[id1]
return flux0
def __red(alan,al0,ga,c1,c2,c3,c4):
fun1 = lambda x: c3/(((x-(al0**2/x))**2)+ga*ga)
fun2 = lambda x,cc: cc*(0.539*((x-5.9)**2)+0.0564*((x-5.9)**3))
fun = lambda x,cc: c1+c2*x+fun1(x)+fun2(x,cc)
ala = alan*1.0e-04 #wavelength in microns
p = 1.0/ala
if np.size(p)>1:
p1, p2 = p[p>=5.9], p[p<5.9]
ff = np.append(fun(p1,c4), fun(p2,0.0))
elif np.size(p)==1:
if p<5.9: c4 = 0.0
ff = fun(p, c4)
else:
return
return ff
def reddening(sw, sf, av=0.0, model=0):
"""
calculate the intrinsic extinction of a given template
Parameters:
sw: array
wavelength
sf: array
flux
av: float or array
model: int
Five models will be used:
1: Allen (1976) for the Milky Way
2: Seaton (1979) fit by Fitzpatrick (1986) for the Milky Way
3: Fitzpatrick (1986) for the Large Magellanic Cloud (LMC)
4: Prevot et al (1984) and Bouchet (1985) for the Small Magellanic Cloud (SMC)
5: Calzetti et al (2000) for starburst galaxies
6: Reddy et al (2015) for star forming galaxies
Return:
reddening-corrected flux or observed flux
"""
if model==0 or av==0.0:
flux=sf
elif model==1: # Allen (1976) for the Milky Way
lambda0 = np.array([1000, 1110, 1250, 1430, 1670, \
2000, 2220, 2500, 2850, 3330, \
3650, 4000, 4400, 5000, 5530, \
6700, 9000, 10000, 20000, 100000], dtype=float)
kR = np.array([4.20, 3.70, 3.30, 3.00, 2.70, \
2.80, 2.90, 2.30, 1.97, 1.69, \
1.58, 1.45, 1.32, 1.13, 1.00, \
0.74, 0.46, 0.38, 0.11, 0.00],dtype=float)
ext0 = InterpolatedUnivariateSpline(lambda0, kR, k=1)
A_lambda = av*ext0(sw)
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
elif model==2: # Seaton (1979) fit by Fitzpatrick (1986) for the Milky Way
Rv=3.1
al0, ga, c1, c2, c3, c4 = 4.595, 1.051, -0.38, 0.74, 3.96, 0.26
ff11 = __red(1100.0,al0,ga,c1,c2,c3,c4)
ff12 = __red(1200.0,al0,ga,c1,c2,c3,c4)
slope=(ff12-ff11)/100.0
lambda0 = np.array([3650, 4000, 4400, 5000, 5530, \
6700, 9000, 10000, 20000, 100000], dtype=float)
kR = np.array([1.58, 1.45, 1.32, 1.13, 1.00, \
0.74, 0.46, 0.38, 0.11, 0.00],dtype=float)
fun = interp1d(lambda0, kR, kind='linear')
sw0 = sw[sw<1200.0]
A_lambda0 = (ff11+(sw0-1100.0)*slope)/Rv+1.0
sw1 = sw[np.logical_and(sw>=1200.0, sw<=3650.0)]
ff = __red(sw1,al0,ga,c1,c2,c3,c4)
A_lambda1 = ff/Rv+1.0
sw2 = sw[np.logical_and(sw>3650.0, sw<=100000.0)]
A_lambda2 = fun(sw2)
A_lambda3 = sw[sw>100000.0]*0.0
A_lambda = av*np.hstack([A_lambda0,A_lambda1,A_lambda2,A_lambda3])
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
elif model==3: # Fitzpatrick (1986) for the Large Magellanic Cloud (LMC)
Rv=3.1
al0, ga, c1, c2, c3, c4 = 4.608, 0.994, -0.69, 0.89, 2.55, 0.50
ff11 = __red(1100.0,al0,ga,c1,c2,c3,c4)
ff12 = __red(1200.0,al0,ga,c1,c2,c3,c4)
slope=(ff12-ff11)/100.0
lambda0 = np.array([3330, 3650, 4000, 4400, 5000, 5530, \
6700, 9000, 10000, 20000, 100000], dtype=float)
kR = np.array([1.682, 1.58, 1.45, 1.32, 1.13, 1.00, \
0.74, 0.46, 0.38, 0.11, 0.00],dtype=float)
fun = interp1d(lambda0, kR, kind='linear')
sw0 = sw[sw<1200.0]
A_lambda0 = (ff11+(sw0-1100.0)*slope)/Rv+1.0
sw1 = sw[np.logical_and(sw>=1200.0, sw<=3330.0)]
ff = __red(sw1,al0,ga,c1,c2,c3,c4)
A_lambda1 = ff/Rv+1.0
sw2 = sw[np.logical_and(sw>3330.0, sw<=100000.0)]
A_lambda2 = fun(sw2)
A_lambda3 = sw[sw>100000.0]*0.0
A_lambda = av*np.hstack([A_lambda0,A_lambda1,A_lambda2,A_lambda3])
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
elif model==4: # Prevot et al (1984) and Bouchet (1985) for the Small Magellanic Cloud (SMC)
Rv = 2.72
lambda0 = np.array([1275, 1330, 1385, 1435, 1490, 1545, \
1595, 1647, 1700, 1755, 1810, 1860, \
1910, 2000, 2115, 2220, 2335, 2445, \
2550, 2665, 2778, 2890, 2995, 3105, \
3704, 4255, 5291, 12500, 16500, 22000], dtype=float)
kR = np.array([13.54, 12.52, 11.51, 10.80, 9.84, 9.28, \
9.06, 8.49, 8.01, 7.71, 7.17, 6.90, 6.76, \
6.38, 5.85, 5.30, 4.53, 4.24, 3.91, 3.49, \
3.15, 3.00, 2.65, 2.29, 1.81, 1.00, 0.00, \
-2.02, -2.36, -2.47],dtype=float)
kR = kR/Rv+1.0
ext0 = InterpolatedUnivariateSpline(lambda0, kR, k=1)
A_lambda = av*ext0(sw)
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
elif model==5: # Calzetti et al (2000) for starburst galaxies
Rv = 4.05
sw = sw*1.0e-04 #wavelength in microns
fun1 = lambda x: 2.659*(-2.156+1.509/x-0.198/x**2+0.011/x**3)+Rv
fun2 = lambda x: 2.659*(-1.857+1.040/x)+Rv
ff11, ff12 = fun1(0.11), fun1(0.12)
slope1=(ff12-ff11)/0.01
ff99, ff100 = fun2(2.19), fun2(2.2)
slope2=(ff100-ff99)/0.01
sw0 = sw[sw<0.12]
sw1 = sw[np.logical_and(sw>=0.12, sw<=0.63)]
sw2 = sw[np.logical_and(sw>0.63, sw<=2.2)]
sw3 = sw[sw>2.2]
k_lambda0 = ff11+(sw0-0.11)*slope1
k_lambda1, k_lambda2 = fun1(sw1), fun2(sw2)
k_lambda3 = ff99+(sw3-2.19)*slope2
A_lambda = av*np.hstack([k_lambda0,k_lambda1,k_lambda2,k_lambda3])/Rv
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
elif model==6: # Reddy et al (2015) for satr forming galaxies
Rv = 2.505
sw = sw*1.0e-04
fun1 = lambda x: -5.726+4.004/x-0.525/x**2+0.029/x**3+Rv
fun2 = lambda x: -2.672-0.010/x+1.532/x**2-0.412/x**3+Rv
ff11, ff12 = fun1(0.14), fun1(0.15)
slope1=(ff12-ff11)/0.01
ff99, ff100 = fun2(2.84), fun2(2.85)
slope2=(ff100-ff99)/0.01
sw0 = sw[sw<0.15]
sw1 = sw[np.logical_and(sw>=0.15, sw<0.60)]
sw2 = sw[np.logical_and(sw>=0.60, sw<2.85)]
sw3 = sw[sw>=2.85]
k_lambda0 = ff11+(sw0-0.14)*slope1
k_lambda1, k_lambda2 = fun1(sw1), fun2(sw2)
k_lambda3 = ff99+(sw3-2.84)*slope2
A_lambda = av*np.hstack([k_lambda0,k_lambda1,k_lambda2,k_lambda3])/Rv
A_lambda[A_lambda<0.0] = 0.0
flux = sf*10**(-0.4*A_lambda)
else:
print("!!! Please select a proper reddening model")
sys.exit(0)
return flux
def getObservedSED(sedCat, redshift=0.0, av=0.0, redden=0):
z = redshift + 1.0
sw, sf = sedCat[:,0], sedCat[:,1]
# reddening
sf = reddening(sw, sf, av=av, model=redden)
sw, sf = sw*z, sf*(z**3)
# lyman forest correction
sf = lyman_forest(sw, sf, redshift)
isedObs = (sw.copy(), sf.copy())
return isedObs
def getSEDData(sedDir='',sedType = 0):
sedListF = open(sedDir + 'galaxy.list')
sedIter = 1
l=''
while sedIter<=sedType:
l = sedListF.readline()
sedIter +=1
sedfn = l.split()[0]
sedData = loadtxt(sedDir + sedfn)
return sedData
def tag_sed(starSpecfile, model_tag, teff=5000, logg=2, feh=0):
h5file = h5.File(starSpecfile,'r')
model_tag_str = model_tag.decode("utf-8").strip()
teff_grid = np.unique(h5file["teff"][model_tag_str])
logg_grid = np.unique(h5file["logg"][model_tag_str])
feh_grid = np.unique(h5file["feh"][model_tag_str])
close_teff = teff_grid[np.argmin(np.abs(teff_grid - teff))]
close_logg = logg_grid[np.argmin(np.abs(logg_grid - logg))]
if model_tag_str == "koester" or model_tag_str == "MC":
close_feh = 99
else:
close_feh = feh_grid[np.argmin(np.abs(feh_grid - feh))]
path = model_tag_str + f"_teff_{close_teff:.1f}_logg_{close_logg:.2f}_feh_{close_feh:.1f}"
wave = np.array(h5file["wave"][model_tag_str][()]).ravel()
flux = np.array(h5file["sed"][path][()]).ravel()
return path, wave, flux
def getABMagAverageVal(ABmag=20.,norm_thr=None, sWave=6840, eWave=8250):
"""
norm_thr: astropy.table, 2 colum, 'WAVELENGTH', 'SENSITIVITY'
Return:
the integerate flux of AB magnitude in the norm_thr(the throughtput of band), photos s-1 m-2 A-1
"""
inverseLambda = norm_thr['SENSITIVITY']/norm_thr['WAVELENGTH']
norm_thr_i = interpolate.interp1d(norm_thr['WAVELENGTH'], inverseLambda)
x = np.linspace(sWave,eWave, int(eWave)-int(sWave)+1)
y = norm_thr_i(x)
AverageLamdaInverse = np.trapz(y,x)/(eWave-sWave)
norm = 54798696332.52474 * pow(10.0, -0.4 * ABmag) * AverageLamdaInverse
# print('AverageLamdaInverse = ', AverageLamdaInverse)
# print('norm = ', norm)
return norm
def getNormFactorForSpecWithABMAG(ABMag=20., spectrum=None, norm_thr=None, sWave=6840, eWave=8250):
"""
Use AB magnitude system (zero point, fv = 3631 janskys) in the normal band(norm_thr) normalize the spectrum by inpute ABMag
Parameters
----------
spectrum: astropy.table, 2 colum, 'WAVELENGTH', 'FLUX'
norm_thr: astropy.table, 2 colum, 'WAVELENGTH', 'SENSITIVITY'
sWave: the start of norm_thr
eWave: the end of norm_thr
Return:
the normalization factor flux of AB system(fix a band and magnitude ) /the flux of inpute spectrum(fix a band)
"""
spectrumi = interpolate.interp1d(spectrum['WAVELENGTH'], spectrum['FLUX'])
norm_thri = interpolate.interp1d(norm_thr['WAVELENGTH'], norm_thr['SENSITIVITY'])
x = np.linspace(sWave,eWave, int(eWave)-int(sWave)+1)
y_spec = spectrumi(x)
y_thr = norm_thri(x)
y = y_spec*y_thr
specAve = np.trapz(y,x)/(eWave-sWave)
norm = getABMagAverageVal(ABmag=ABMag, norm_thr=norm_thr, sWave=sWave, eWave=eWave)
if specAve == 0:
return 0
# print('specAve = ', specAve)
return norm / specAve
def getABMAG(spec, bandpass_fn):
throughtput = Table.read(bandpass_fn)
thr_ids = throughtput['SENSITIVITY'] > 0.01
sWave = np.floor(throughtput[thr_ids][0][0])
eWave = np.ceil(throughtput[thr_ids][-1][0])
# sWave=2000
# eWave = 18000
# print('in getABMAG', sWave, eWave)
spectrumi = interpolate.interp1d(spec['WAVELENGTH'], spec['FLUX'])
thri = interpolate.interp1d(throughtput['WAVELENGTH'],throughtput['SENSITIVITY'])
x = np.linspace(sWave,eWave, (int(eWave)-int(sWave)+1))
y_spec = spectrumi(x)*thri(x)
interFlux = np.trapz(y_spec, x)
ABMAG_zero = getABMagAverageVal(
ABmag=0,
norm_thr=throughtput,
sWave=sWave,
eWave=eWave)
flux_ave = interFlux / (eWave-sWave)
ABMAG_spec = -2.5 * np.log10(flux_ave/ABMAG_zero)
return ABMAG_spec
def getABMAGANDErr(spec, bandpass_fn, readout = 5, sky = 0.2, dark = 0.02, t = 150, aper = 2, frame = 1, noisepix_num = 22, flux_ratio = 1.0):
throughtput = Table.read(bandpass_fn)
thr_ids = throughtput['SENSITIVITY'] > 0.01
sWave = np.floor(throughtput[thr_ids][0][0])
eWave = np.ceil(throughtput[thr_ids][-1][0])
# sWave=2000
# eWave = 18000
# print('in getABMAG', sWave, eWave)
spectrumi = interpolate.interp1d(spec['WAVELENGTH'], spec['FLUX'])
thri = interpolate.interp1d(throughtput['WAVELENGTH'],throughtput['SENSITIVITY'])
x = np.linspace(sWave,eWave, (int(eWave)-int(sWave)+1))
y_spec = spectrumi(x)*thri(x)
interFlux = np.trapz(y_spec, x)
ABMAG_zero = getABMagAverageVal(
ABmag=0,
norm_thr=throughtput,
sWave=sWave,
eWave=eWave)
flux_ave = interFlux / (eWave-sWave)
ABMAG_spec = -2.5 * np.log10(flux_ave/ABMAG_zero)
totalFlux = interFlux * t * frame * math.pi*(aper*aper/4.)
noise_var = totalFlux*flux_ratio + (sky + dark) * t * frame * noisepix_num + readout * readout * noisepix_num * frame
mag_err = 1.087*(noise_var/(totalFlux*flux_ratio))
return ABMAG_spec, mag_err
def getAveEnerge(spec, bandpass_fn):
throughtput = Table.read(bandpass_fn)
thr_ids = throughtput['SENSITIVITY'] > 0.01
sWave = np.floor(throughtput[thr_ids][0][0])
eWave = np.ceil(throughtput[thr_ids][-1][0])
# sWave=2000
# eWave = 18000
# print('in getABMAG', sWave, eWave)
spectrumi = interpolate.interp1d(spec['WAVELENGTH'], spec['FLUX'])
thri = interpolate.interp1d(throughtput['WAVELENGTH'],throughtput['SENSITIVITY'])
x = np.linspace(sWave,eWave, (int(eWave)-int(sWave)+1))
y_spec = spectrumi(x)
interFlux = np.trapz(y_spec, x)
return interFlux/(eWave-sWave)
def produceSourceSED(sedSoureType = 0,mag_norm = 24.0, norm_filter_thr_fn= 'g.fits',gal_sed_lib_dir = 'Galaxy/', z=0.1, av = 0.1, redden = 0, gal_sedType =1,star_sed_lib_fn='SpecLib.hdf5', lib_tag = 'phoe',teff = 5000, logg = 2, feh = 0):
if sedSoureType == 0:
tag = lib_tag.encode('UTF-8')
_, wave, flux = tag_sed(star_sed_lib_fn, tag, teff=teff, logg=logg, feh=feh)
elif sedSoureType==1:
sedData = getSEDData(gal_sed_lib_dir, sedType = gal_sedType)
sed_data = getObservedSED(
sedCat=sedData,
redshift=z,
av=av,
redden=redden)
wave = sed_data[0]
flux = sed_data[1]
speci = interpolate.interp1d(wave, flux)
lamb = np.arange(2000, 18001 + 0.5, 0.5)
y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
orig_spec_phot = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
normThr = Table.read(norm_filter_thr_fn)
# orig_spec = Table(np.array([wave,flux]).T, names=(['WAVELENGTH', 'FLUX']))
norm_thr_rang_ids = normThr['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag_norm, spectrum=orig_spec_phot,
norm_thr=normThr,
sWave=np.floor(normThr[norm_thr_rang_ids][0][0]),
eWave=np.ceil(normThr[norm_thr_rang_ids][-1][0]))
norm_spec = Table(Table(np.array([wave,flux*sedNormFactor]).T, names=(['WAVELENGTH', 'FLUX'])))
norm_spec_phot = Table(Table(np.array([lamb,all_sed*sedNormFactor]).T, names=(['WAVELENGTH', 'FLUX'])))
return norm_spec, norm_spec_phot
def calculatCSSTMAG_ERR(spec = None, throughput_dir = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/CSST/', t = 150,frame = 1, noisepix_num = 22, flux_ratio = 1.0):
fil = ['nuv','u','g','r','i','z','y']
skybg = {'nuv': 0.00261,'u':0.01823,'g':0.15897,'r':0.20705,'i':0.21433,'z':0.12658,'y':0.037}
resMag = {}
for fi in fil:
mag,err = getABMAGANDErr(spec, throughput_dir+fi+'.Throughput.fits', readout = 5, sky = skybg[fi], dark = 0.02, t = t, aper = 2, frame = frame, noisepix_num = noisepix_num, flux_ratio = flux_ratio)
resMag[fi] = [mag,err]
return resMag
def calculatCSSTMAG(spec = None, throughput_dir = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/CSST/'):
fil = ['nuv','u','g','r','i','z','y']
resMag = {}
for fi in fil:
mag = getABMAG(spec, throughput_dir+fi+'.Throughput.fits')
resMag[fi] = mag
return resMag
def calculatCSSTFilEnergy(spec = None, throughput_dir = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/CSST/'):
fil = ['nuv','u','g','r','i','z','y']
resMag = {}
for fi in fil:
ene = getAveEnerge(spec, throughput_dir+fi+'.Throughput.fits')
resMag[fi] = ene
return resMag
def produceGalSED_C6( gal_id_s = '03593100052300144566', gal_z = 1.6927,mag_norm = 24.0, norm_filter_thr_fn= 'g.fits',galaxy_cat_dir = "/Volumes/EAGET/C6_data/inputData/Catalog_C6_20221212/cat2CSSTSim_bundle/",sedlib_dir = "/Volumes/EAGET/C6_data/inputData/Catalog_C6_20221212/sedlibs/"):
healPix_id = int(gal_id_s[0:6])
galcat_file = galaxy_cat_dir + "galaxies_C6_bundle" + gal_id_s[6:12] + '.h5'
g_id = int(gal_id_s[12:])
gals_cat = h5.File(galcat_file, 'r')['galaxies']
coeff = gals_cat[str(healPix_id)]['coeff'][:][g_id]
pcs = h5.File(os.path.join(sedlib_dir, "pcs.h5"), "r")
lamb = h5.File(os.path.join(sedlib_dir, "lamb.h5"), "r")
lamb_gal = lamb['lamb'][()]
pcs = pcs['pcs'][()]
cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
factor = 10**(-.4 * cosmo.distmod(gal_z).value)
flux = np.matmul(pcs, coeff) * factor
flux[flux < 0] = 0.
sedcat = np.vstack((lamb_gal, flux)).T
sed_data = getObservedSED(
sedCat=sedcat,
redshift=gal_z,
av=0.0,
redden=0.0
)
wave, flux = sed_data[0], sed_data[1]
speci = interpolate.interp1d(wave, flux)
lamb = np.arange(2000, 18001 + 0.5, 0.5)
y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
orig_spec_phot = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
normThr = Table.read(norm_filter_thr_fn)
# orig_spec = Table(np.array([wave,flux]).T, names=(['WAVELENGTH', 'FLUX']))
norm_thr_rang_ids = normThr['SENSITIVITY'] > 0.001
sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag_norm, spectrum=orig_spec_phot,
norm_thr=normThr,
sWave=np.floor(normThr[norm_thr_rang_ids][0][0]),
eWave=np.ceil(normThr[norm_thr_rang_ids][-1][0]))
norm_spec = Table(Table(np.array([wave,flux*sedNormFactor]).T, names=(['WAVELENGTH', 'FLUX'])))
norm_spec_phot = Table(Table(np.array([lamb,all_sed*sedNormFactor]).T, names=(['WAVELENGTH', 'FLUX'])))
return norm_spec, norm_spec_phot
fileDir = os.getcwd()
#normlization filter star
# star_norm_fn = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/SDSS/SLOAN_SDSS.g.fits'
star_norm_fn = os.path.join(fileDir, "data/throughputs/SDSS/SLOAN_SDSS.g.fits")
#恒星模板库
star_sed_lib = "/Volumes/ExtremeSSD/SimData/Catalog_20210126/SpecLib.hdf5"
#输入参数，星等，得到两个光谱，spec是能量，spec_p是光子
spec, spec_photo = produceSourceSED(sedSoureType = 0,mag_norm = 20., norm_filter_thr_fn = star_norm_fn, star_sed_lib_fn=star_sed_lib, lib_tag = 'MM',teff = 3800., logg = 0. , feh = -1.)
#星系星表文件
galaxy_cat_dir = "/Volumes/EAGET/C6_data/inputData/Catalog_C6_20221212/cat2CSSTSim_bundle/"
#星系光谱模板，PCA
sedlib_dir = "/Volumes/EAGET/C6_data/inputData/Catalog_C6_20221212/sedlibs/"
#normlization filter galaxy
# gal_norm_fn = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/LSST/lsst_throuput_g.fits'
gal_norm_fn = os.path.join(fileDir, "data/throughputs/LSST/lsst_throuput_g.fits")
gal_spec, gal_spec_photo = produceGalSED_C6(gal_id_s = '03593100052300144566', gal_z = 1.6927,mag_norm = 24.0, norm_filter_thr_fn= gal_norm_fn)
#根据上面计算的光谱计算csst星等，噪声不需要就不用管了，t, frame, noisepix_num, flux_ratio,都是为了估计噪声
# mags = calculatCSSTMAG_ERR(spec = spec_photo,throughput_dir = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/CSST/', t = 150,frame = 2, noisepix_num = 22, flux_ratio = 1.0)
# # #打印csst星等
# for k in list(mags.keys()):
# print(k, mags[k][0])
# gal_norm_fn = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/LSST/lsst_throuput_g.fits'
# gal_sed_dir = "/Volumes/Extreme SSD/SimData/Templates/Galaxy/"
# # spec1, spec1_p = produceSourceSED(sedSoureType = 1,mag_norm = 22.075, norm_filter_thr_fn= gal_norm_fn,gal_sed_lib_dir = gal_sed_dir, z=0.1, av = 0.1, redden = 0)
# spec1, spec1_p = produceSourceSED(sedSoureType = 1,mag_norm = 22.075, norm_filter_thr_fn= gal_norm_fn,gal_sed_lib_dir = gal_sed_dir, z=0.35, av = 0.35, redden = 3.0000,gal_sedType=22)
# fil = ['nuv','u','g','r','i','z','y']
# throughput_dir = '/Users/zhangxin/Work/SlitlessSim/sed/produceSED_bycatfile/data/throughputs/CSST/'
# for fi in fil:
# throughput_fn = throughput_dir + fi + '.throughput.fits'
# mag = getABMAG(spec_p, throughput_dir+fi+'.Throughput.fits')
# print(fi,mag)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment