Skip to content
GitLab
Commit 471e6c0b authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

update

parent 19fd5ed7
Show whitespace changes
Inline Side-by-side
csst_ifs_sim/csst_ifs_sim.py
View file @ 471e6c0b
......@@ -138,7 +138,7 @@ class CDM03bidir():
# read in trap information
trapdata = np.loadtxt(
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(
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]
......@@ -261,10 +261,7 @@ class CDM03bidir():
return np.asanyarray(CTIed)
#################################################################################################################
#################################################################################################################
"""
These functions can be used for logging information.
......@@ -272,8 +269,6 @@ These functions can be used for logging information.
:version: 0.3
"""
def lg(log_filename, loggername='logger'):
"""
Sets up a logger.
......@@ -302,8 +297,6 @@ def lg(log_filename, loggername='logger'):
logger.addHandler(handler)
return logger
##############################################################################
def IFSinformation():
......@@ -658,7 +651,7 @@ def ill2flux(path, E):
"""
# use template from sky_bkg (background_spec_hst.dat)
filename = 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
......@@ -752,7 +745,7 @@ class StrayLight(object):
self.ecliptic = self.equator.transform_to('barycentrictrueecliptic')
self.pointing = transRaDec2D(radec[0], radec[1])
self.slcdll = ctypes.CDLL(
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(
......@@ -775,11 +768,11 @@ 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 = self.path+"IFS_inputdata/refs/DE405"
self.PSTFn = self.path+"IFS_inputdata/refs/PST"
self.RFn = self.path + "IFS_inputdata/refs/R"
self.ZolFn = self.path+"IFS_inputdata/refs/Zodiacal"
self.brightStarTabFn = 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))
......@@ -1865,9 +1858,6 @@ class IFSsimulator():
# get solar position
dt = datetime.fromisoformat(time)
jd = time2jd(dt)
#
#
#
# # jd = julian.to_jd(dt, fmt='jd')
t = Time(jd, format='jd', scale='utc')
......@@ -1888,7 +1878,7 @@ class IFSsimulator():
# interpolated zodical surface brightness at 0.5 um
zodi = pd.read_csv(
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])
......@@ -1902,7 +1892,7 @@ class IFSsimulator():
# read the zodical spectrum in the ecliptic
cat_spec = pd.read_csv(
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 #
......@@ -1922,42 +1912,6 @@ class IFSsimulator():
return wave_A, spec_erg2
##########################################################################
# def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)):
# """
# Parameters
# ----------
# image : TYPE
# DESCRIPTION.
# sigma : TYPE, optional
# DESCRIPTION. The default is (0.32, 0.32).
# Returns
# -------
# TYPE
# DESCRIPTION.
# """
# """
# Smooths a given image with a gaussian kernel with widths given as sigmas.
# This smoothing can be used to mimic charge diffusion within the CCD.
# The default values are from Table 8-2 of CCD_273_Euclid_secification_1.0.130812.pdf converted
# to sigmas (FWHM / (2sqrt(2ln2)) and rounded up to the second decimal.
# .. Note:: This method should not be called for the full image if the charge spreading
# has already been taken into account in the system PSF to avoid double counting.
# :param image: image array which is smoothed with the kernel
# :type image: ndarray
# :param sigma: widths of the gaussian kernel that approximates the charge diffusion [0.32, 0.32].
# :param sigma: tuple
# :return: smoothed image array
# :rtype: ndarray
# """
# return ndimage.filters.gaussian_filter(image, sigma)
###############################################################################
def readCosmicRayInformation(self):
......@@ -1979,9 +1933,9 @@ class IFSsimulator():
# print(self.information)
crLengths = np.loadtxt(
self.information['dir_path']+self.information['cosmicraylengths'])
os.path.join(self.information['dir_path'],self.information['cosmicraylengths']))
crDists = np.loadtxt(
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])
......@@ -2038,31 +1992,21 @@ class IFSsimulator():
else:
ss = '_'
# if self.information['dir_path'] == '/nfsdata/share/simulation-unittest/ifs_sim/':
# self.result_path = self.information['dir_path'] + \
# self.information['result_path']+'/'+self.source+ss+result_day
# else:
# home_path = os.environ['HOME']
# if home_path == '/home/yan':
# self.result_path = '../IFS_simData/'+self.source+ss+result_day
# else:
# self.result_path = '/data/ifspip/CCD_ima/'+self.source+ss+result_day
self.result_path = self.information['result_path'] + \
'/'+self.source+ss+result_day
print(self.information['result_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)
print(self.result_path)
if os.path.isdir(self.result_path) is False:
os.mkdir(self.result_path)
os.mkdir(self.result_path+'/calibration_Data')
os.mkdir(self.result_path+'/log_file')
os.mkdir(self.result_path+'/original_Calibration')
os.mkdir(self.result_path+'/original_sky')
os.mkdir(self.result_path+'/shift_rot_sky')
os.mkdir(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()
......@@ -2070,8 +2014,8 @@ class IFSsimulator():
data_time = now.strftime("%Y-%m-%d-%H-%M-%S")
############################################################
self.log = lg(self.result_path+'/log_file/IFS_' +
self.source+'_'+data_time+'_Num_'+str(simnumber)+'.log')
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')
......@@ -2092,7 +2036,7 @@ class IFSsimulator():
self.HgArsigma = self.light_FWHM/2.35 # sigma value of the light source
# load system optical and CCD efficiency data
matfn0 = 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
......@@ -2103,7 +2047,7 @@ class IFSsimulator():
# load all useful data;
# load wavefront data;
matfn2 = 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
......@@ -2111,17 +2055,17 @@ class IFSsimulator():
self.pupil = abs(opd0) > 0.0
####################
matfn2 = 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 = 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 = 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
......@@ -2424,50 +2368,50 @@ class IFSsimulator():
self.log.info('Added dark current to bule and red channel')
if self.information['dark1_b'] > 0.001 or self.information['dark1_b'] > 0.001:
if self.information['dark1_b'] < 0.0001 or self.information['dark1_b'] > 0.001:
self.log.error(
'dark1_b value error, it shoub be in [0.0001, 0.001]!')
raise ValueError(
'dark1_b value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark2_b'] > 0.001 or self.information['dark2_b'] > 0.001:
if self.information['dark2_b'] < 0.0001 or self.information['dark2_b'] > 0.001:
self.log.error(
'dark2_b value error, it shoub be in [0.0001, 0.001]!')
raise ValueError(
'dark2_b value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark3_b'] > 0.001 or self.information['dark3_b'] > 0.001:
if self.information['dark3_b'] < 0.0001 or self.information['dark3_b'] > 0.001:
self.log.error(
'dark3_b value error, it shoub be in [0.0001, 0.001]!')
raise ValueError(
'dark3_b value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark4_b'] > 0.001 or self.information['dark4_b'] > 0.001:
if self.information['dark4_b'] < 0.0001 or self.information['dark4_b'] > 0.001:
self.log.error(
'dark4_b value error, it shoub be in [0.0001, 0.001]!')
raise ValueError(
'dark4_b value error, it shoub be in [0.0001, 0.001]!')
####################################################################
if self.information['dark1_r'] > 0.001 or self.information['dark1_r'] > 0.001:
if self.information['dark1_r'] < 0.0001 or self.information['dark1_r'] > 0.001:
self.log.error(
'dark1_r value error, it shoub be in [0.0001, 0.001]!')
raise ValueError(
'dark1_r value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark2_r'] > 0.001 or self.information['dark2_r'] > 0.001:
if self.information['dark2_r'] < 0.0001 or self.information['dark2_r'] > 0.001:
self.log.error(
'dark2_r value error, it shoub be in [0.0001, 0.001]!')
raise ValueError(
'dark2_r value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark3_r'] > 0.001 or self.information['dark3_r'] > 0.001:
if self.information['dark3_r'] < 0.0001 or self.information['dark3_r'] > 0.001:
self.log.error(
'dark3_r value error, it shoub be in [0.0001, 0.001]!')
raise ValueError(
'dark3_r value error, it shoub be in [0.0001, 0.001]!')
if self.information['dark4_r'] > 0.001 or self.information['dark4_r'] > 0.001:
if self.information['dark4_r'] < 0.0001 or self.information['dark4_r'] > 0.001:
self.log.error(
'dark4_r value error, it shoub be in [0.0001, 0.001]!')
raise ValueError(
......@@ -2558,7 +2502,7 @@ class IFSsimulator():
Warning:: This method does not work if the input file has exactly one line.
"""
cosmetics = np.loadtxt(
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']))
......@@ -2582,7 +2526,7 @@ class IFSsimulator():
#############################################################################
cosmetics = np.loadtxt(
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']))
......@@ -3257,48 +3201,7 @@ class IFSsimulator():
self.log.info('Maximum and total values of the image are %i and %i, respectively' % (np.max(self.image_r),
np.sum(self.image_r)))
# ############################################################################
# def applyImageShift(self):
# """
# Returns
# -------
# None.
# """
# np.random.seed(9*self.simnumber)
# ud = np.random.random() # Choose a random rotation
# dx = 2 * (ud-0.5) * self.information['shiftmax']
# np.random.seed(99*self.simnumber)
# ud = np.random.random() # Choose a random rotation
# dy = 2 * (ud-0.5) * self.information['shiftmax']
# self.image_b = ndimage.shift(self.image_b.copy(), [
# dy+self.information['shift_b_y'], dx+self.information['shift_b_x']], order=0, mode='nearest')
# self.image_r = ndimage.shift(self.image_r.copy(), [
# dy+self.information['shift_r_y'], dx+self.information['shift_r_x']], order=0, mode='nearest')
# self.log.info('Applied image shifting to g r i channels.')
# self.information['ra'] = dx*self.information['pixel_size']
# self.information['dec'] = dy*self.information['pixel_size']
##############################################################################
# def applyImageRotate(self):
# """
# Returns
# -------
# None.
# """
# np.random.seed(10*self.simnumber)
# ud = np.random.random() # Choose a random rotation
# angle = 2 * (ud-0.5) * self.information['tel_rotmax']
# inputimg = self.image_b.copy()
# # here we choose reshape=False, the rotated image will
# rotimg = ndimage.rotate(
# inputimg, angle+self.information['rotate_b'], order=1, reshape=False)
# self.image_b = rotimg
# inputimg = self.image_r.copy()
# # here we choose reshape=False, the rotated image will
# rotimg = ndimage.rotate(
# inputimg, angle+self.information['rotate_r'], order=1, reshape=False)
# self.image_r = rotimg
# self.information['Tel_rot'] = angle
# self.log.info(
# 'Applied telescope rotation with angle (in degree)= %f.', angle)
###############################################################################
def CCDreadout(self):
......@@ -3518,6 +3421,7 @@ 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'
filename_r = 'CSST_IFS_R_'+self.source+'_'+exp_start_str + \
......@@ -3537,11 +3441,11 @@ 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+'/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 = 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 ####
......@@ -3566,7 +3470,7 @@ class IFSsimulator():
'ICRS', 'coordinate system of the object')
ofd_b.header['EQUINOX'] = (float(2000.0), '')
ofd_b.header['FITSSWV'] = (
'ifs_sim_0.8.03', 'FITS creating software version')
'csst_ifs_sim_3.0.0', 'FITS creating software version')
# ######## Object information #############
if self.source == 'SCI' or self.source == 'COMP':
ofd_b.header['OBJECT'] = (
......@@ -3933,7 +3837,7 @@ class IFSsimulator():
ofd_r.header['EQUINOX'] = (float(2000.0), '')
ofd_r.header['FITSSWV'] = (
'ifs_sim_0.8.03', 'FITS creating software version')
'csst_ifs_sim_3.0.0', 'FITS creating software version')
# ######## Object information #############
if self.source == 'SCI' or self.source == 'COMP':
......@@ -4377,19 +4281,19 @@ class IFSsimulator():
"""
# read solar template
solar_template = pd.read_csv(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(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(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
......@@ -4507,23 +4411,21 @@ class IFSsimulator():
################################################################
# ########
# #### load slicer_Qe1K.fits ####
slicerfile = 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 = 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 = 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
......@@ -4772,10 +4674,10 @@ class IFSsimulator():
data_time = now.strftime("%Y-%m-%d")
# write the actual file
bluefile = self.result_path+'/original_Calibration/IFS_'+self.source + \
'_'+data_time+'_'+str(self.simnumber)+'_B_original.fits'
redfile = 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)
......@@ -4795,7 +4697,7 @@ class IFSsimulator():
newd = fits.HDUList([hdu1, hdu2])
filename = 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)
......@@ -4836,8 +4738,7 @@ class IFSsimulator():
############################################################################
# consider the slice optical efficiency in different slicer channel
slicer_QE_file = 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
......@@ -5041,9 +4942,7 @@ 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 = 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
......@@ -5220,10 +5119,8 @@ 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(self.result_path+'/shift_rot_sky/rot_shift_subImg_' +
str(simnumber)+'.fits', overwrite=True)
fits.writeto(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 ###
......@@ -5342,22 +5239,6 @@ class IFSsimulator():
energy = energy+sum(photons.flux[idx0])
p_num = len(idx0[0])
##############################################################
# ######### code for test slice is right or not ##############
##############################################################
# photons_2=galsim.PhotonArray(N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux = photons.flux[idx0], wavelength=lam)
# if k==0:
# sliceimage = galsim.Image(100,100)
# sliceimage.scale = pixelscale
# sliceimage.setOrigin(0,0)
# #image.setCenter(int(np.mean(photons.x)), int(np.mean(photons.y)))
# rotphotons=galsim.PhotonArray(N=p_num, x=photons_2.x/pixelscale, y=photons_2.y/pixelscale, flux = photons_2.flux)
# sensor = galsim.Sensor()
# sensor.accumulate(rotphotons, sliceimage)
# ######### finish testing #############################
# #############################################################
# ####################################
# find photons for blue channel, and make the flux multiple the optical and CCD efficiency
......@@ -5404,17 +5285,6 @@ class IFSsimulator():
########################################################
# ##################### test code ################################
# ################################################################
# fits.writeto('originalImg.fits',Nimg,overwrite=True)
# fits.writeto('PSFedImg.fits',temp,overwrite=True)
# fits.writeto('rotedphotonsImg.fits',image.array,overwrite=True)
# fits.writeto('newsub_rotedphotonsImg.fits',newimage.array,overwrite=True)
# fits.writeto('slicephotonsImg.fits',sliceimage.array,overwrite=True)
# pause()
# #############################################################
# ############## finish test code ###########################
#####################################################################
# stray light will cover 2% of input total light;
......@@ -5430,10 +5300,10 @@ class IFSsimulator():
self.image_r = red_img.array
obsid = simnumber
bluefile = self.result_path+'/original_sky/IFS_' + \
self.source+'_'+str(obsid)+'_B.fits'
redfile = 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)
......@@ -5463,14 +5333,13 @@ class IFSsimulator():
# load HgAr data;
if self.source == 'LAMP':
matfn3 = 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 = 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
......@@ -5639,7 +5508,7 @@ class IFSsimulator():
# print('time=',(end-start))
# consider the slice optical efficiency in different slicer channel
da = fits.open(
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 ###################
......@@ -5698,22 +5567,6 @@ class IFSsimulator():
p_num = len(idx0[0])
# #############################################################
# ##### code for test slice is right or not ##################
#########################################################
# photons_2=galsim.PhotonArray(N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux = photons.flux[idx0], wavelength=lam)
# if k==0:
# sliceimage = galsim.Image(100,100)
# sliceimage.scale = pixelscale
# sliceimage.setOrigin(0,0)
# #image.setCenter(int(np.mean(photons.x)), int(np.mean(photons.y)))
# rotphotons=galsim.PhotonArray(N=p_num, x=photons_2.x/pixelscale, y=photons_2.y/pixelscale, flux = photons_2.flux)
# sensor = galsim.Sensor()
# sensor.accumulate(rotphotons, sliceimage)
# ######################## finish testing #######################################################
# ######################################################################################################
#####################################
# find photons for blue channel, and make the flux multiple the optical and CCD efficiency
......@@ -5769,10 +5622,11 @@ class IFSsimulator():
data_time = now.strftime("%Y-%m-%d")
# write the actual file
bluefile = self.result_path+'/original_Calibration/IFS_'+self.source + \
'_'+data_time+'_'+str(self.simnumber)+'_B_original.fits'
redfile = 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)
......@@ -5792,7 +5646,7 @@ class IFSsimulator():
newd = fits.HDUList([hdu1, hdu2])
filename = 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)
......@@ -5861,8 +5715,7 @@ class IFSsimulator():
if simnumber <= 200 and simnumber > 150:
self.information['exptime'] = 1200
self.skyfilepath = 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)
......@@ -5879,7 +5732,7 @@ class IFSsimulator():
##########################################
df = pd.read_csv(
self.information['dir_path']+'IFS_inputdata/TianCe/'+starcat)
os.path.join(self.information['dir_path'],'IFS_inputdata/TianCe/'+starcat))
##############################################################
sn = self.simnumber-1
......@@ -5905,8 +5758,8 @@ class IFSsimulator():
# fn=father_path+'/IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt';
fn = self.information['dir_path'] + \
'IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt'
fn = os.path.join(self.information['dir_path'],
'IFS_inputdata/TianCe/orbit20160925/'+str(k)+'.txt')
d = np.loadtxt(fn)
for kk in range(len(d[:, 0])):
......@@ -5928,8 +5781,8 @@ class IFSsimulator():
else:
fn = 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
......@@ -6164,20 +6017,15 @@ def runIFSsim(sourcein, configfile, dir_path, result_path, iLoop, debug, applyho
simulate[iLoop] = IFSsimulator(configfile)
simulate[iLoop].configure(sourcein, dir_path, result_path,
iLoop, debug, applyhole) # load the configfile;
iLoop, debug, applyhole)
# load the configfile;
simulate[iLoop].information['result_path'] = result_path
if applyhole == 'yes' and sourcein == 'LAMP':
simulate[iLoop].information['holemask'] = 'yes'
else:
simulate[iLoop].information['holemask'] = 'no'
# ## dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
# simulate[iLoop].information['dir_path'] = dir_path
# simulate[iLoop].information['debug'] = debug
# simulate[iLoop].information['result_path'] = result_path
simulate[iLoop].simulate(sourcein, iLoop)
return 1
......
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