''' 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 # 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, sw=None, ew = None, bpass_type = 'fits'): if bpass_type == 'fits': throughtput = Table.read(bandpass_fn) else: th_dat = np.loadtxt(bandpass_fn) throughtput = Table(th_dat, names=['WAVELENGTH','SENSITIVITY']) if sw is None or ew is None: thr_ids = throughtput['SENSITIVITY'] > 0.01 sWave = np.floor(throughtput[thr_ids][0][0]) eWave = np.ceil(throughtput[thr_ids][-1][0]) else: sWave=sw eWave = ew if throughtput['WAVELENGTH'][0]>sWave: throughtput.insert_row(0,[sWave,0.0]) if throughtput['WAVELENGTH'][-1] 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 produceNormSED_photon(inputSED = None,mag_norm = 24.0, norm_filter_thr_fn= 'g.fits',ws = 2000, we = 18000): wave = inputSED['WAVELENGTH'] flux = inputSED['FLUX'] speci = interpolate.interp1d(wave, flux) lamb = np.arange(ws, we + 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_phot = Table(Table(np.array([lamb,all_sed*sedNormFactor]).T, names=(['WAVELENGTH', 'FLUX']))) return 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/',ws = None, we = None, filelist = None, band_instr = 'CSST'): fil = ['nuv','u','g','r','i','z','y'] if filelist is not None: fil = filelist resMag = {} for fi in fil: if band_instr == 'CSST': mag = getABMAG(spec, throughput_dir+fi+'.Throughput.fits',sw = ws, ew = we, bpass_type = 'fits') else: mag = getABMAG(spec, throughput_dir+fi+'.dat',sw = ws, ew = we, bpass_type = 'dat') 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 = -99, 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=1 if gal_z != 0: 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, 11001 + 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')) sedNormFactor=1.0 if mag_norm != -99: 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 = '03490800052300010462', gal_z = 0.3764,mag_norm = -99, 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)