Skip to content
GitLab
Commit 33dff68d authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

update

parent 7f51aaa2
Pipeline #7593 failed with stage
in 0 seconds
Show whitespace changes
Inline Side-by-side
csst_ifs_sim/csst_ifs_sim.py
View file @ 33dff68d
......@@ -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']))
......@@ -3422,11 +3420,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 = 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,9 +5310,9 @@ class IFSsimulator():
self.image_r = red_img.array
obsid = simnumber
bluefile = os.path.join(self.result_path,'original_sky/IFS_' + \
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_' + \
redfile = os.path.join(self.result_path, 'original_sky/IFS_' +
self.source+'_'+str(obsid)+'_R.fits')
fits.writeto(bluefile, blue_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,10 +5634,10 @@ 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 + \
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 + \
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)
......@@ -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
......@@ -5781,7 +5795,7 @@ class IFSsimulator():
else:
fn = os.path.join(self.information['dir_path'] ,
fn = os.path.join(self.information['dir_path'],
'IFS_inputdata/TianCe/orbit20160925/'+str(k+1)+'.txt')
d = np.loadtxt(fn)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment