From 33dff68d405f99a1d24595fe183d7a8de72c4067 Mon Sep 17 00:00:00 2001 From: yan Date: Tue, 24 Dec 2024 09:23:14 +0800 Subject: [PATCH] update --- csst_ifs_sim/csst_ifs_sim.py | 174 +++++++++++++++++++---------------- 1 file changed, 94 insertions(+), 80 deletions(-) diff --git a/csst_ifs_sim/csst_ifs_sim.py b/csst_ifs_sim/csst_ifs_sim.py index 3115c36..913226e 100644 --- a/csst_ifs_sim/csst_ifs_sim.py +++ b/csst_ifs_sim/csst_ifs_sim.py @@ -138,7 +138,7 @@ class CDM03bidir(): # read in trap information trapdata = np.loadtxt( - os.path.join(self.values['dir_path'],self.values['paralleltrapfile'])) + os.path.join(self.values['dir_path'], self.values['paralleltrapfile'])) if trapdata.ndim > 1: self.nt_p = trapdata[:, 0] self.sigma_p = trapdata[:, 1] @@ -150,7 +150,7 @@ class CDM03bidir(): self.taur_p = [trapdata[2],] trapdata = np.loadtxt( - os.path.join(self.values['dir_path'],self.values['serialtrapfile'])) + os.path.join(self.values['dir_path'], self.values['serialtrapfile'])) if trapdata.ndim > 1: self.nt_s = trapdata[:, 0] self.sigma_s = trapdata[:, 1] @@ -260,15 +260,7 @@ class CDM03bidir(): return np.asanyarray(CTIed) -################################################################################################################# -################################################################################################################# -""" -These functions can be used for logging information. - -.. Warning:: logger is not multiprocessing safe. -:version: 0.3 -""" def lg(log_filename, loggername='logger'): """ Sets up a logger. @@ -297,7 +289,7 @@ def lg(log_filename, loggername='logger'): logger.addHandler(handler) return logger -############################################################################## + def IFSinformation(): """ @@ -651,7 +643,7 @@ def ill2flux(path, E): """ # use template from sky_bkg (background_spec_hst.dat) - filename = os.path.join(path,'IFS_inputdata/refs/background_spec_hst.dat') + filename = os.path.join(path, 'IFS_inputdata/refs/background_spec_hst.dat') cat_spec = pd.read_csv(filename, sep='\\s+', header=None, comment='#') wave0 = cat_spec[0].values # A spec0 = cat_spec[2].values # erg/s/cm^2/A/arcsec^2 @@ -745,7 +737,7 @@ class StrayLight(object): self.ecliptic = self.equator.transform_to('barycentrictrueecliptic') self.pointing = transRaDec2D(radec[0], radec[1]) self.slcdll = ctypes.CDLL( - os.path.join(self.path,'IFS_inputdata/refs/libstraylight.so')) # dylib + os.path.join(self.path, 'IFS_inputdata/refs/libstraylight.so')) # dylib self.slcdll.Calculate.argtypes = [ctypes.c_double, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.POINTER( @@ -768,11 +760,12 @@ class StrayLight(object): ctypes.POINTER(ctypes.c_double)] self.slcdll.Init.argtypes = [ ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p] - self.deFn = os.path.join(self.path,"IFS_inputdata/refs/DE405") - self.PSTFn = os.path.join(self.path,"IFS_inputdata/refs/PST") - self.RFn = os.path.join(self.path,"IFS_inputdata/refs/R") - self.ZolFn = os.path.join(self.path,"IFS_inputdata/refs/Zodiacal") - self.brightStarTabFn = os.path.join(self.path,"IFS_inputdata/refs/BrightGaia_with_csst_mag") + self.deFn = os.path.join(self.path, "IFS_inputdata/refs/DE405") + self.PSTFn = os.path.join(self.path, "IFS_inputdata/refs/PST") + self.RFn = os.path.join(self.path, "IFS_inputdata/refs/R") + self.ZolFn = os.path.join(self.path, "IFS_inputdata/refs/Zodiacal") + self.brightStarTabFn = os.path.join( + self.path, "IFS_inputdata/refs/BrightGaia_with_csst_mag") self.slcdll.Init(str.encode(self.deFn), str.encode( self.PSTFn), str.encode(self.RFn), str.encode(self.ZolFn)) @@ -1878,7 +1871,7 @@ class IFSsimulator(): # interpolated zodical surface brightness at 0.5 um zodi = pd.read_csv( - os.path.join(self.information['dir_path'],'IFS_inputdata/refs/zodi_map.dat'), sep='\\s+', header=None, comment='#') + os.path.join(self.information['dir_path'], 'IFS_inputdata/refs/zodi_map.dat'), sep='\\s+', header=None, comment='#') beta_angle = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75]) lamda_angle = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180]) @@ -1892,7 +1885,7 @@ class IFSsimulator(): # read the zodical spectrum in the ecliptic cat_spec = pd.read_csv( - os.path.join(self.information['dir_path'],'IFS_inputdata/refs/solar_spec.dat'), sep='\\s+', header=None, comment='#') + os.path.join(self.information['dir_path'], 'IFS_inputdata/refs/solar_spec.dat'), sep='\\s+', header=None, comment='#') wave = cat_spec[0].values # A spec0 = cat_spec[1].values # zodi_norm = 252 # @@ -1933,9 +1926,9 @@ class IFSsimulator(): # print(self.information) crLengths = np.loadtxt( - os.path.join(self.information['dir_path'],self.information['cosmicraylengths'])) + os.path.join(self.information['dir_path'], self.information['cosmicraylengths'])) crDists = np.loadtxt( - os.path.join(self.information['dir_path'],self.information['cosmicraydistance'])) + os.path.join(self.information['dir_path'], self.information['cosmicraydistance'])) self.cr = dict(cr_u=crLengths[:, 0], cr_cdf=crLengths[:, 1], cr_cdfn=np.shape(crLengths)[0], cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0]) @@ -1991,22 +1984,23 @@ class IFSsimulator(): else: ss = '_' - - - ### create result father path; + # + # + # ##create result father path; if os.path.isdir(self.information['result_path']) is False: os.mkdir(self.information['result_path']) - self.result_path=os.path.join(self.information['result_path'], self.source+ss+result_day) + self.result_path = os.path.join( + self.information['result_path'], self.source+ss+result_day) print(self.result_path) if os.path.isdir(self.result_path) is False: os.mkdir(self.result_path) - os.mkdir(os.path.join(self.result_path,'calibration_Data')) - os.mkdir(os.path.join(self.result_path,'log_file')) - os.mkdir(os.path.join(self.result_path,'original_Calibration')) - os.mkdir(os.path.join(self.result_path,'original_sky')) - os.mkdir(os.path.join(self.result_path,'shift_rot_sky')) - os.mkdir(os.path.join(self.result_path,'sky_Data')) + os.mkdir(os.path.join(self.result_path, 'calibration_Data')) + os.mkdir(os.path.join(self.result_path, 'log_file')) + os.mkdir(os.path.join(self.result_path, 'original_Calibration')) + os.mkdir(os.path.join(self.result_path, 'original_sky')) + os.mkdir(os.path.join(self.result_path, 'shift_rot_sky')) + os.mkdir(os.path.join(self.result_path, 'sky_Data')) ############################################################################## now = datetime.utcnow() @@ -2014,7 +2008,7 @@ class IFSsimulator(): data_time = now.strftime("%Y-%m-%d-%H-%M-%S") ############################################################ - self.log = lg(os.path.join(self.result_path,'log_file/IFS_' + + self.log = lg(os.path.join(self.result_path, 'log_file/IFS_' + self.source+'_'+data_time+'_Num_'+str(simnumber)+'.log')) self.log.info('STARTING A NEW SIMULATION') @@ -2036,7 +2030,8 @@ class IFSsimulator(): self.HgArsigma = self.light_FWHM/2.35 # sigma value of the light source # load system optical and CCD efficiency data - matfn0 = os.path.join(self.information['dir_path'],'IFS_inputdata/TotalQ200923.mat') + matfn0 = os.path.join( + self.information['dir_path'], 'IFS_inputdata/TotalQ200923.mat') self.log.info('Optical optical efficiency file path is: %s' % (matfn0)) da0 = sio.loadmat(matfn0) # optical efficiency of blue channel @@ -2047,7 +2042,8 @@ class IFSsimulator(): # load all useful data; # load wavefront data; - matfn2 = os.path.join(self.information['dir_path'],'IFS_inputdata/opd/opd_638nm.mat') + matfn2 = os.path.join( + self.information['dir_path'], 'IFS_inputdata/opd/opd_638nm.mat') self.log.info('OPD0 file path is: %s' % (matfn2)) da2 = sio.loadmat(matfn2) opd0 = da2['opd'] # opd unit is in meter @@ -2055,17 +2051,20 @@ class IFSsimulator(): self.pupil = abs(opd0) > 0.0 #################### - matfn2 = os.path.join(self.information['dir_path'],'IFS_inputdata/opd/opd1.fits') + matfn2 = os.path.join( + self.information['dir_path'], 'IFS_inputdata/opd/opd1.fits') da = fits.open(matfn2) self.opd1 = da[0].data self.log.info('OPD1 file path is: %s' % (matfn2)) - matfn2 = os.path.join(self.information['dir_path'],'IFS_inputdata/opd/opd2.fits') + matfn2 = os.path.join( + self.information['dir_path'], 'IFS_inputdata/opd/opd2.fits') self.log.info('OPD2 file path is: %s' % (matfn2)) da = fits.open(matfn2) self.opd2 = da[0].data - matfn2 = os.path.join(self.information['dir_path'],'IFS_inputdata/opd/opd3.fits') + matfn2 = os.path.join( + self.information['dir_path'], 'IFS_inputdata/opd/opd3.fits') self.log.info('OPD3 file path is: %s' % (matfn2)) da = fits.open(matfn2) self.opd3 = da[0].data @@ -2452,7 +2451,6 @@ class IFSsimulator(): # ############################################################################# def applyPoissonNoise(self): - """ Returns ------- @@ -2502,7 +2500,7 @@ class IFSsimulator(): Warning:: This method does not work if the input file has exactly one line. """ cosmetics = np.loadtxt( - os.path.join(self.information['dir_path'],self.information['cosmeticsfile_b'])) + os.path.join(self.information['dir_path'], self.information['cosmeticsfile_b'])) self.log.info('cosmeticsfile_b path is: %s' % (self.information['cosmeticsfile_b'])) @@ -2526,7 +2524,7 @@ class IFSsimulator(): ############################################################################# cosmetics = np.loadtxt( - os.path.join(self.information['dir_path'],self.information['cosmeticsfile_r'])) + os.path.join(self.information['dir_path'], self.information['cosmeticsfile_r'])) self.log.info('cosmeticsfile_r path is: %s' % (self.information['cosmeticsfile_r'])) @@ -3421,12 +3419,14 @@ class IFSsimulator(): filename_b = 'CSST_IFS_B_'+self.source+'_'+exp_start_str + \ '_'+exp_end_str+'_'+str(obsid)+'_B_L0_V'+sim_ver - - self.file_b = self.result_path+'/sky_Data/'+filename_b+'.fits' + + self.file_b = os.path.join( + self.result_path, 'sky_Data/'+filename_b+'.fits') filename_r = 'CSST_IFS_R_'+self.source+'_'+exp_start_str + \ '_'+exp_end_str+'_'+str(obsid)+'_R_L0_V'+sim_ver - self.file_r = self.result_path+'/sky_Data/'+filename_r + '.fits' + self.file_r = os.path.join( + self.result_path, 'sky_Data/'+filename_r + '.fits') else: @@ -3441,11 +3441,13 @@ class IFSsimulator(): filename_b = 'CSST_IFS_B_'+self.source+'_'+exp_start_str + \ '_'+exp_end_str+'_'+str(obsid)+'_B_L0_V'+sim_ver - self.file_b = os.path.join(self.result_path,'calibration_Data/'+filename_b+'.fits') + self.file_b = os.path.join( + self.result_path, 'calibration_Data/'+filename_b+'.fits') filename_r = 'CSST_IFS_R_'+self.source+'_'+exp_start_str + \ '_'+exp_end_str+'_'+str(obsid)+'_R_L0_V'+sim_ver - self.file_r = os.path.join(self.result_path,'calibration_Data/'+filename_r+'.fits') + self.file_r = os.path.join( + self.result_path, 'calibration_Data/'+filename_r+'.fits') # create a new FITS file, using HDUList instance # #### layer 0 #### @@ -4281,19 +4283,19 @@ class IFSsimulator(): """ # read solar template - solar_template = pd.read_csv(os.path.join(self.information['dir_path'],'IFS_inputdata/refs/solar_spec.dat'), sep='\\s+', + solar_template = pd.read_csv(os.path.join(self.information['dir_path'], 'IFS_inputdata/refs/solar_spec.dat'), sep='\\s+', header=None, comment='#') template_wave = solar_template[0].values template_flux = solar_template[1].values # read earth shine surface brightness - earthshine_curve = pd.read_csv(os.path.join(self.information['dir_path'],'IFS_inputdata/refs/earthshine.dat'), + earthshine_curve = pd.read_csv(os.path.join(self.information['dir_path'], 'IFS_inputdata/refs/earthshine.dat'), header=None, comment='#') angle = earthshine_curve[0].values surface_brightness = earthshine_curve[1].values # read V-band throughtput - cat_filter_V = pd.read_csv(os.path.join(self.information['dir_path'],'IFS_inputdata/refs/filter_Bessell_V.dat'), sep='\\s+', + cat_filter_V = pd.read_csv(os.path.join(self.information['dir_path'], 'IFS_inputdata/refs/filter_Bessell_V.dat'), sep='\\s+', header=None, comment='#') filter_wave = cat_filter_V[0].values filter_response = cat_filter_V[1].values @@ -4411,21 +4413,24 @@ class IFSsimulator(): ################################################################ # ######## # #### load slicer_Qe1K.fits #### - slicerfile = os.path.join(self.information['dir_path'],'IFS_inputdata/Flatfield/slicer_QE_1K_230625.fits') + slicerfile = os.path.join( + self.information['dir_path'], 'IFS_inputdata/Flatfield/slicer_QE_1K_230625.fits') da = fits.open(slicerfile) self.log.info('hole Sim and slicer 1K QE file:%s' % (slicerfile)) slicer_Qe = da[0].data # #### load hole mask ###### # load hole mask matrix - file = os.path.join(self.information['dir_path'], 'IFS_inputdata/Hole/holemask_230612.fits') + file = os.path.join( + self.information['dir_path'], 'IFS_inputdata/Hole/holemask_230612.fits') self.log.info('hole mask file is %s' % (file)) da = fits.open(file) HoleMask = da[self.simnumber].data # load HgAr data; - matfn3 = os.path.join(self.information['dir_path'],'IFS_inputdata/HgAr.mat') + matfn3 = os.path.join( + self.information['dir_path'], 'IFS_inputdata/HgAr.mat') self.log.info('lamp flux file is : %s' % (matfn3)) da3 = sio.loadmat(matfn3) HgArst = da3['HgArst'] # opd unit is in meter @@ -4674,10 +4679,11 @@ class IFSsimulator(): data_time = now.strftime("%Y-%m-%d") # write the actual file - bluefile = os.path.join(self.result_path,'original_Calibration/IFS_'+self.source + \ - '_'+data_time+'_'+str(self.simnumber)+'_B_original.fits') - redfile = os.path.join(self.result_path,'original_Calibration/IFS_'+self.source + \ - '_'+data_time+'_'+str(self.simnumber)+'_R_original.fits') + bluefile = os.path.join(self.result_path, 'original_Calibration/IFS_' + + self.source + '_' + d ata_time + '_' + + str(self.simnumber)+'_B_original.fits') + redfile = os.path.join(self.result_path, 'original_Calibration/IFS_' + + self.source+'_'+data_time+'_'+str(self.simnumber)+'_R_original.fits') fits.writeto(bluefile, blue_img.array, overwrite=True) fits.writeto(redfile, red_img.array, overwrite=True) @@ -4697,7 +4703,8 @@ class IFSsimulator(): newd = fits.HDUList([hdu1, hdu2]) - filename = os.path.join(self.result_path,'calibration_Data/'+self.source+'_flux.fits') + filename = os.path.join( + self.result_path, 'calibration_Data/'+self.source+'_flux.fits') if self.simnumber == 1: newd.writeto(filename, overwrite=True) @@ -4738,7 +4745,8 @@ class IFSsimulator(): ############################################################################ # consider the slice optical efficiency in different slicer channel - slicer_QE_file = os.path.join(self.information['dir_path'],'IFS_inputdata/Flatfield/slicer_QE230625.fits') + slicer_QE_file = os.path.join( + self.information['dir_path'], 'IFS_inputdata/Flatfield/slicer_QE230625.fits') da = fits.open(slicer_QE_file) self.log.info('slicer_QE_file path is: %s' % (slicer_QE_file)) slicer_Qe = da[0].data @@ -4942,7 +4950,8 @@ class IFSsimulator(): # orbit parameters are not in currenct txt file if self.orbit_pars[-1, 0] < t2jd: self.orbit_file_num = self.orbit_file_num+1 - fn = os.path.join(self.information['dir_path'] , 'IFS_inputdata/TianCe/orbit20160925/'+str(self.orbit_file_num)+'.txt') + fn = os.path.join( + self.information['dir_path'], 'IFS_inputdata/TianCe/orbit20160925/'+str(self.orbit_file_num)+'.txt') self.orbit_pars = np.loadtxt(fn) self.orbit_exp_num = 0 @@ -5119,8 +5128,10 @@ class IFSsimulator(): # PSFfilename='/media/yan//IFSsim/IFSsim Data/rot_shift_Img.fits' # fits.writeto('/media/yan//IFSsim/IFSsim Data/rot_shift_subImg.fits',image0.array[50-32:50+32,50-32:50+32],overwrite=True ) - hdu1.writeto(os.path.join(self.result_path,'shift_rot_sky/rot_shift_subImg_'+str(simnumber)+'.fits'), overwrite=True) - fits.writeto(os.path.join(self.result_path,'shift_rot_sky/original_Img_'+str(simnumber)+'.fits'), Nimg, overwrite=True) + hdu1.writeto(os.path.join( + self.result_path, 'shift_rot_sky/rot_shift_subImg_'+str(simnumber)+'.fits'), overwrite=True) + fits.writeto(os.path.join( + self.result_path, 'shift_rot_sky/original_Img_'+str(simnumber)+'.fits'), Nimg, overwrite=True) ##################################################################### # ## do convolve image0 with PSF0 from primay CSST ### @@ -5285,7 +5296,6 @@ class IFSsimulator(): ######################################################## - ##################################################################### # stray light will cover 2% of input total light; if self.sky_noise: @@ -5300,10 +5310,10 @@ class IFSsimulator(): self.image_r = red_img.array obsid = simnumber - bluefile = os.path.join(self.result_path,'original_sky/IFS_' + \ - self.source+'_'+str(obsid)+'_B.fits') - redfile = os.path.join(self.result_path,'original_sky/IFS_' + \ - self.source+'_'+str(obsid)+'_R.fits') + bluefile = os.path.join(self.result_path, 'original_sky/IFS_' + + self.source+'_'+str(obsid)+'_B.fits') + redfile = os.path.join(self.result_path, 'original_sky/IFS_' + + self.source+'_'+str(obsid)+'_R.fits') fits.writeto(bluefile, blue_img.array, overwrite=True) fits.writeto(redfile, red_img.array, overwrite=True) @@ -5333,13 +5343,15 @@ class IFSsimulator(): # load HgAr data; if self.source == 'LAMP': - matfn3 = os.path.join(self.information['dir_path'],'IFS_inputdata/HgAr.mat') + matfn3 = os.path.join( + self.information['dir_path'], 'IFS_inputdata/HgAr.mat') self.log.info('lamp flux file is : %s' % (matfn3)) da3 = sio.loadmat(matfn3) HgArst = da3['HgArst'] # opd unit is in meter # load flat data if self.source == 'FLAT': - matfn3 = os.path.join(self.information['dir_path'],'IFS_inputdata/flat_light.mat') + matfn3 = os.path.join( + self.information['dir_path'], 'IFS_inputdata/flat_light.mat') self.log.info('flat flux file is : %s' % (matfn3)) da3 = sio.loadmat(matfn3) flat_light = da3['flat_light'] # opd unit is in meter @@ -5508,7 +5520,7 @@ class IFSsimulator(): # print('time=',(end-start)) # consider the slice optical efficiency in different slicer channel da = fits.open( - os.path.join(self.information['dir_path'],'IFS_inputdata/Flatfield/slicer_QE230625.fits')) + os.path.join(self.information['dir_path'], 'IFS_inputdata/Flatfield/slicer_QE230625.fits')) slicer_Qe = da[0].data img0 = img0*slicer_Qe # ######### do the slice effect ################### @@ -5622,11 +5634,11 @@ class IFSsimulator(): data_time = now.strftime("%Y-%m-%d") # write the actual file - bluefile = os.path.join(self.result_path,'original_Calibration/IFS_'+self.source + \ - '_'+data_time+'_'+str(self.simnumber)+'_B_original.fits') - - redfile = os.path.join(self.result_path,'original_Calibration/IFS_'+self.source + \ - '_'+data_time+'_'+str(self.simnumber)+'_R_original.fits') + bluefile = os.path.join(self.result_path, 'original_Calibration/IFS_'+self.source + + '_'+data_time+'_'+str(self.simnumber)+'_B_original.fits') + + redfile = os.path.join(self.result_path, 'original_Calibration/IFS_'+self.source + + '_'+data_time+'_'+str(self.simnumber)+'_R_original.fits') fits.writeto(bluefile, blue_img.array, overwrite=True) fits.writeto(redfile, red_img.array, overwrite=True) @@ -5646,7 +5658,8 @@ class IFSsimulator(): newd = fits.HDUList([hdu1, hdu2]) - filename = os.path.join(self.result_path,'calibration_Data/'+self.source+'_flux.fits') + filename = os.path.join( + self.result_path, 'calibration_Data/'+self.source+'_flux.fits') if self.simnumber == 1: newd.writeto(filename, overwrite=True) @@ -5715,7 +5728,8 @@ class IFSsimulator(): if simnumber <= 200 and simnumber > 150: self.information['exptime'] = 1200 - self.skyfilepath = os.path.join(self.information['dir_path'],self.information['sky_fitsin']) + self.skyfilepath = os.path.join( + self.information['dir_path'], self.information['sky_fitsin']) print('self.skyfilepath = ', self.skyfilepath) @@ -5732,7 +5746,7 @@ class IFSsimulator(): ########################################## df = pd.read_csv( - os.path.join(self.information['dir_path'],'IFS_inputdata/TianCe/'+starcat)) + os.path.join(self.information['dir_path'], 'IFS_inputdata/TianCe/'+starcat)) ############################################################## sn = self.simnumber-1 @@ -5759,7 +5773,7 @@ class IFSsimulator(): # fn=father_path+'/IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt'; fn = os.path.join(self.information['dir_path'], - 'IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt') + 'IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt') d = np.loadtxt(fn) for kk in range(len(d[:, 0])): @@ -5781,8 +5795,8 @@ class IFSsimulator(): else: - fn = os.path.join(self.information['dir_path'] , - 'IFS_inputdata/TianCe/orbit20160925/'+str(k+1)+'.txt') + fn = os.path.join(self.information['dir_path'], + 'IFS_inputdata/TianCe/orbit20160925/'+str(k+1)+'.txt') d = np.loadtxt(fn) self.orbit_pars = d -- GitLab