produceSED_1.py

'''
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)