produce one spec by parameters

produceSED.py

@@ -352,7 +352,7 @@ def getABMAG(spec, bandpass_fn, sw=None, ew = None, bpass_type = 'fits'):
     # 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))
+    x = np.linspace(sWave,eWave, 1000)
     y_spec = spectrumi(x)*thri(x)
     interFlux = np.trapz(y_spec, x)
 
@@ -471,6 +471,30 @@ def produceNormSED_photon(inputSED = None,mag_norm = 24.0, norm_filter_thr_fn= '
     norm_spec_phot = Table(Table(np.array([lamb,all_sed*sedNormFactor]).T,  names=(['WAVELENGTH', 'FLUX'])))
     return norm_spec_phot
 
+def produceNormSED_spec(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]))
+    # print("###################",sedNormFactor, mean(orig_spec_phot['FLUX']))
+    # sedNormFactor = getNormFactorForSpecWithABMAG(ABMag = mag_norm, spectrum = orig_spec_phot, norm_thr=normThr, sWave=2000,eWave=11000)
+                                                    
+    
+    norm_spec = Table(Table(np.array([lamb,y*sedNormFactor]).T,  names=(['WAVELENGTH', 'FLUX'])))
+    return norm_spec
+
 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}

pycode.py

@@ -16,6 +16,7 @@ import ctypes
 from astropy.table import Table
 import produceSED
 import astropy.constants as atcons
+import sys
 # from ctypes import *
 
 # struct STAR
@@ -75,8 +76,19 @@ print(spec[500:509])
 '''
 
 from astropy.table import Table
-# 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"
+
+if len(sys.argv) < 2:
+	print('usage:\n'+sys.argv[0]+' survey_result.txt')
+	sys.exit(0)
+
+catalogFn = str(sys.argv[1])
+
+# catalogFn = "/nfsdata/share/CSSTsimInputCat_TH/C9_RA240_DECp30.fits"
+# catalogFn = "/home/zhangxin/CSST_SIM/star_spec/csst_spec_interp_clean/code/data/catalog/trilegal.fits"
+outdir = '/nfsdata/share/CSSOSDataProductsSims/trilegalCat/'
+# outFn = "C9_RA240_DECp30_refCat_1.fits"
+
+outFn = catalogFn.strip().split('/')[-1][0:-5] + "_refCat_1.fits"
 cat = Table.read(catalogFn)
 
 filters = ['nuv','u','g','r','i','z','y']
@@ -95,12 +107,12 @@ from mpi4py import MPI
 comm = MPI.COMM_WORLD
 rank = comm.Get_rank()
 rank_size = comm.Get_size()
-print("------------------",rank_size, rank)
+print("------------------",rank_size, rank, nrows)
 iterNum = 0    
 for star in cat:
     if iterNum%10000==0:
         print(iterNum)
-    # if iterNum > 10000:
+    # if iterNum > 100:
     #     iterNum = iterNum + 1
     #     continue
     if iterNum % rank_size != rank:
@@ -135,7 +147,31 @@ for star in cat:
     parallaxs[iterNum] = 1/(10**(star['mwmsc_mu0']*0.2)/100)
     iterNum = iterNum + 1
     # print(mags, mags_others)
-    
+
+# send_data1  = res
+# send_data2 = parallaxs
+# if rank == 0:
+#     total_res = comm.gather(send_data1, root=0)
+#     total_parall = comm.gather(send_data2, root=0)
+
+#     for i in np.arange(1, rank_size, 1):
+#         for fi in filters:
+#             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)
+#     outdir = '/nfsdata/share/CSSOSDataProductsSims/trilegalCat/'
+#     outFn = "C9_RA300_DECm60_refCat.fits"
+#     cat.write(os.path.join(outdir,outFn),overwrite=True)
+
+# else:
+#     comm.gather(send_data1, root=0)
+#     comm.gather(send_data2, root=0)
 
 
 total_res=comm.gather(res, root=0)
@@ -156,8 +192,8 @@ if rank == 0:
         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)
-    outdir = '/nfsdata/share/CSSOSDataProductsSims/trilegalCat/'
-    outFn = "trilegal_cal_mag_n.fits"
+    # outdir = '/nfsdata/share/CSSOSDataProductsSims/trilegalCat/'
+    # outFn = "C9_RA240_DECp30_refCat_1.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'])

pycode_cal_sigle.py

+'''
+Author: Zhang Xin zhangx@bao.ac.cn
+Date: 2024-01-02 13:34:39
+LastEditors: Zhang Xin zhangx@bao.ac.cn
+LastEditTime: 2024-03-25 13:51:33
+FilePath: /csst-simulation/Users/zhangxin/Work/SlitlessSim/CSST_SIM/Star_spec/csst_spec_interp_clean/code/pycode.py
+Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE
+'''
+##运行下面在mac/linux下执行 
+# cc -fPIC -shared main_singlestar.c -I/usr/include -I/usr/include/cfitsio -lcfitsio -lm -o test.dylib
+
+import os
+import numpy as np
+from astropy.io import fits
+import ctypes
+from astropy.table import Table
+import produceSED
+import astropy.constants as atcons
+import sys
+# from ctypes import *
+
+# struct STAR
+# {
+#   float logte;
+#   float logg;
+#   //float logL;
+#   float Mass;
+#   float Av;
+#   float mu0;
+#   float Z;
+#   //float mbolmag;
+# };
+
+class Star(ctypes.Structure):
+    _fields_ = [
+               ('logte',ctypes.c_float),
+               ('logg',ctypes.c_float),
+               ('Mass',ctypes.c_float),
+               ('Av', ctypes.c_float),
+               ('mu0', ctypes.c_float),
+               ('Z', ctypes.c_float)]
+
+#CHANGE
+#topdir='/run/media/chen/1TS/dupe/gitlab/csst_spec_interp/'
+topdir='/home/zhangxin/CSST_SIM/star_spec/csst_spec_interp_clean/code/'
+code_name='test.so'
+
+d = ctypes.CDLL(os.path.join(topdir,code_name))
+d.loadSpecLibs.argtypes=[ctypes.c_char_p, ctypes.c_char_p]
+d.loadExts.argtypes=[ctypes.c_char_p]
+
+
+#######################################################################################
+#CHANGE
+#nwv = d.loadSpecLibs(str.encode(os.path.join(topdir,'file_CK04.par')),str.encode(topdir))
+#d.loadExts(str.encode(os.path.join(topdir,"spec/Ext_odonnell94_R3.1_CK04W.fits")))
+#TO
+nwv = d.loadSpecLibs(str.encode(os.path.join(topdir,'file_BT-Settl_CSST_wl1000-24000_R1000.par')),str.encode(topdir))
+d.loadExts(str.encode(os.path.join(topdir,"spec/Ext_odonnell94_R3.1_CSST_wl1000-24000_R1000.fits")))
+########################################################################################
+
+print("Done loading Exts")
+spec = (ctypes.c_float*nwv)()
+wave = (ctypes.c_float*nwv)()
+d.interpSingleStar.argtypes=[ctypes.Structure, ctypes.POINTER(ctypes.c_float)]
+
+
+'''
+#example for a single star
+#s=Star(3.620752, 4.701155, 0.599979, 0.067540, 11.200001, 0.008501)
+s=Star(3.9739766, 8.108173, 0.6525841, 0.077022046, 9.05, 0.024376422)
+#s=Star(3.9739766, 4.99, 0.6525841, 0.077022046, 9.05, 0.024376422)
+d.interpSingleStar(s, spec)
+spec_ = spec[:]
+print(spec[500:509])
+'''
+
+from astropy.table import Table
+
+if len(sys.argv) < 10:
+	print('need 9 parameters: logte logg mass av mu0 Z rv normMag normFilter')
+	sys.exit(0)
+
+logte = float(sys.argv[1])
+logg = float(sys.argv[2])
+mass = float(sys.argv[3])
+av = float(sys.argv[4])
+mu0 = float(sys.argv[5]) 
+Z = float(sys.argv[6]) 
+rv  = float(sys.argv[7])
+normMag = float(sys.argv[8])
+normFilter = str(sys.argv[9])
+
+filters = ['nuv','u','g','r','i','z','y']
+filters_other = ['2MASS_H','2MASS_J','2MASS_Ks','GAIA_GAIA3.G','GAIA_GAIA3.Gbp','GAIA_GAIA3.Grp','GALEX_GALEX.FUV','GALEX_GALEX.NUV','LSST_u','LSST_g','LSST_r','LSST_i','LSST_z','LSST_y','PAN-STARRS_PS1.g','PAN-STARRS_PS1.r','PAN-STARRS_PS1.i','PAN-STARRS_PS1.z','PAN-STARRS_PS1.y','SLOAN_SDSS.u','SLOAN_SDSS.g','SLOAN_SDSS.r','SLOAN_SDSS.i','SLOAN_SDSS.z']
+
+
+s=Star(logte, logg, mass, av, mu0, Z)
+
+d.interpSingleStar(s, spec, wave)
+
+rv_c = rv/(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_RV, np.power(10,spec[:])]).T, names=('WAVELENGTH', 'FLUX'))
+
+outSpec = produceSED.produceNormSED_spec(inputSED = spec_out,mag_norm = normMag, norm_filter_thr_fn= normFilter,ws = 1050, we = 23950)
+
+writeTxt = True
+
+outfn = 'outspec'
+if writeTxt:
+    fn = outfn + ".txt"
+    with open(fn, 'w') as f:
+        f.writelines([str(w) + '    ' + str(f) + '\n' for w,f in zip(spec_out['WAVELENGTH'],spec_out['FLUX'])])
+else:
+    outSpec.write(outfn+'.fits')
+
+# spec_norm_phot = produceSED.produceNormSED_photon(inputSED = spec_out,mag_norm = magg, 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= 1050, we = 23950, filelist=filters_other, band_instr='other')
+
+
+
+