''' 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, 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')) 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])) sedNormFactor=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 = 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)