renew csst throughput

pycode.py

@@ -15,6 +15,7 @@ from astropy.io import fits
 import ctypes
 from astropy.table import Table
 import produceSED
+import astropy.constants as atcons
 # from ctypes import *
 
 # struct STAR
@@ -74,7 +75,8 @@ print(spec[500:509])
 '''
 
 from astropy.table import Table
-catalogFn = "/nfsdata/nfsdata/share/CSSTsimInputCat_TH/C9_RA300_DECm60.fits"
+# catalogFn = "/nfsdata/share/CSSTsimInputCat_TH/C9_RA300_DECm60.fits"
+catalogFn = "/home/zhangxin/CSST_SIM/star_spec/csst_spec_interp_clean/code/data/catalog/trilegal.fits"
 cat = Table.read(catalogFn)
 
 filters = ['nuv','u','g','r','i','z','y']
@@ -87,6 +89,8 @@ for fi in filters:
 for fi in filters_other:
     res[fi] = np.zeros(nrows)
 
+parallaxs = np.zeros(nrows)
+
 from mpi4py import MPI
 comm = MPI.COMM_WORLD
 rank = comm.Get_rank()
@@ -96,35 +100,46 @@ iterNum = 0
 for star in cat:
     if iterNum%10000==0:
         print(iterNum)
-    if iterNum > 1:
-        iterNum = iterNum + 1
-        break
+    # if iterNum > 10000:
+    #     iterNum = iterNum + 1
+    #     continue
     if iterNum % rank_size != rank:
         iterNum = iterNum + 1
         continue
-    specTable = np.zeros([nwv,2])
-    s=Star(star['mwmsc_logte'], star['mwmsc_logg'], 1., star['mwmsc_av'], star['mwmsc_mu0'], star['mwmsc_z'])
+    # specTable = np.zeros([nwv,2])
+    s=Star(star['mwmsc_logte'], star['mwmsc_logg'], star['mwmsc_mass'], star['mwmsc_av'], star['mwmsc_mu0'], star['mwmsc_z'])
+    # av stellarmass dm teff logg feh
+    # 0.0464   0.7512  10.8000   3.6914   4.5952   0.0122
+    #  s=StarParm(obj.param['teff'], obj.param['logg'], obj.param['stellarMass'], obj.param['av'], obj.param['DM'], obj.param['feh'])
+    # s=Star(3.6914, 4.5952, 0.7512, 0.0464, 10.8000, 0.0122) 
     #print("star[logTe], star[logg], star[Mass], star[Av], star[mu0], star[Z]: ", star['logTe'], star['logg'], star['Mass'], star['Av'], star['mu0'], star['Z'])
     d.interpSingleStar(s, spec, wave)
-    specTable[:,0] = wave[:]
-    specTable[:,1] = spec[:]
+
+    rv_c = star['mwmsc_vrad']/(atcons.c.value/1000.)
+    Doppler_factor = np.sqrt((1+rv_c)/(1-rv_c))
+    wave_RV = wave*Doppler_factor
+    # specTable[:,0] = wave[:]
+    # specTable[:,1] = spec[:]
     # print(spec[500:509])
-    spec_out = Table(np.array([wave[:], np.power(10,spec[:])]).T, names=('WAVELENGTH', 'FLUX'))
-    spec_norm_phot = produceSED.produceNormSED_photon(inputSED = spec_out,mag_norm = star['mwmsc_gsmag'], norm_filter_thr_fn= 'data/throughputs/SDSS/SLOAN_SDSS.g.fits',ws = 1000, we = 24000)
-    mags = produceSED.calculatCSSTMAG(spec = spec_norm_phot, throughput_dir = 'data/throughputs/CSST/',ws= 2000, we = 11000)
+    spec_out = Table(np.array([wave_RV, np.power(10,spec[:])]).T, names=('WAVELENGTH', 'FLUX'))
+    spec_norm_phot = produceSED.produceNormSED_photon(inputSED = spec_out,mag_norm = star['mwmsc_gsmag'], norm_filter_thr_fn= 'data/throughputs/SDSS/SLOAN_SDSS.g.fits',ws = 1050, we = 23950)
+    # spec_norm_phot = produceSED.produceNormSED_photon(inputSED = spec_out,mag_norm = 17.8360, norm_filter_thr_fn= 'data/throughputs/SDSS/SLOAN_SDSS.g.fits',ws = 1000, we = 24000)
+    mags = produceSED.calculatCSSTMAG(spec = spec_norm_phot, throughput_dir = 'data/throughputs/CSST_n/',ws= 2000, we = 11000)
 
-    mags_others = produceSED.calculatCSSTMAG(spec = spec_norm_phot, throughput_dir = 'data/throughputs/filter_transp/',ws= 1000, we = 24000, filelist=filters_other, band_instr='other')
+    mags_others = produceSED.calculatCSSTMAG(spec = spec_norm_phot, throughput_dir = 'data/throughputs/filter_transp/',ws= 1050, we = 23950, filelist=filters_other, band_instr='other')
 
     for fi in filters:
         res[fi][iterNum] = mags[fi]
     for fi in filters_other:
         res[fi][iterNum] = mags_others[fi]
+    parallaxs[iterNum] = 1/(10**(star['mwmsc_mu0']*0.2)/100)
     iterNum = iterNum + 1
-    print(mags, mags_others)
+    # print(mags, mags_others)
     
 
 
 total_res=comm.gather(res, root=0)
+total_parall = comm.gather(parallaxs, root=0)
 
 if rank == 0:
     # new_res = res
@@ -135,13 +150,15 @@ if rank == 0:
             res[fi] = res[fi] + total_res[i][fi]
         for fi in filters_other:
             res[fi] = res[fi] + total_res[i][fi]
-
+        parallaxs = parallaxs + total_parall[i]
+    cat.add_column(np.round(parallaxs,5),name='parallax')
     for fi in filters:
         cat.add_column(np.round(res[fi],5),name='interSpec_'+fi)
     for fi in filters_other:
         cat.add_column(np.round(res[fi],5),name='interSpec_'+fi)
-
-    # cat.write("C9_RA300_DECm60_calmag_.fits",overwrite=True)
+    outdir = '/nfsdata/share/CSSOSDataProductsSims/trilegalCat/'
+    outFn = "trilegal_cal_mag_n.fits"
+    cat.write(os.path.join(outdir,outFn),overwrite=True)
     # print('--------------------------')
     # print(mags['nuv']-star['mwmsc_nuvmag'], mags['u']-star['mwmsc_umag'], mags['g']-star['mwmsc_gmag'], mags['r']-star['mwmsc_rmag'], mags['i']-star['mwmsc_imag'],mags['z']-star['mwmsc_zmag'],mags['y']-star['mwmsc_ymag'])
 #exmple of handling a fits catalouge with many stars