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

update

parent d851b140
Pipeline #7089 failed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
csst_mci_sim/csst_mci_sim.py
View file @ a2bcbefb
...@@ -105,10 +105,10 @@ def get_file_extension(filename): ...@@ -105,10 +105,10 @@ def get_file_extension(filename):
return extension return extension
########################################################## ##########################################################
""" """
Charge Transfer Inefficiency Charge Transfer Inefficiency
============================ ============================
This file contains a simple class to run a CDM03 CTI model developed by Alex Short (ESA). This file contains a simple class to run a CDM03 CTI model developed by Alex Short (ESA).
This now contains both the official CDM03 and a new version that allows different trap This now contains both the official CDM03 and a new version that allows different trap
...@@ -237,7 +237,7 @@ class CDM03bidir(): ...@@ -237,7 +237,7 @@ class CDM03bidir():
:return: image that has been run through the CDM03 model :return: image that has been run through the CDM03 model
:rtype: ndarray """"" :rtype: ndarray """""
#
# return data # return data
iflip = iquadrant / 2 iflip = iquadrant / 2
...@@ -293,7 +293,7 @@ def ill2flux(E, path): ...@@ -293,7 +293,7 @@ def ill2flux(E, path):
# use template from sky_bkg (background_spec_hst.dat) # use template from sky_bkg (background_spec_hst.dat)
filename = path+'MCI_inputData/refs/background_spec_hst.dat' filename = path+'MCI_inputData/refs/background_spec_hst.dat'
cat_spec = pd.read_csv(filename, sep='\s+', header=None, comment='#') cat_spec = pd.read_csv(filename, sep='\\s+', header=None, comment='#')
wave0 = cat_spec[0].values # A wave0 = cat_spec[0].values # A
spec0 = cat_spec[2].values # erg/s/cm^2/A/arcsec^2 spec0 = cat_spec[2].values # erg/s/cm^2/A/arcsec^2
...@@ -350,9 +350,9 @@ def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj): ...@@ -350,9 +350,9 @@ def earth_angle(time_jd, x_sat, y_sat, z_sat, ra_obj, dec_obj):
def CCDnonLinearityModel(data, beta=6e-7): def CCDnonLinearityModel(data, beta=6e-7):
""" """
##
The non-linearity is modelled based on the results presented. The non-linearity is modelled based on the results presented.
:param data: data to which the non-linearity model is being applied to :param data: data to which the non-linearity model is being applied to
:type data: ndarray :type data: ndarray
...@@ -522,7 +522,7 @@ def jd2mjd(jd): ...@@ -522,7 +522,7 @@ def jd2mjd(jd):
def dt2hmd(dt): def dt2hmd(dt):
## dt is datetime # dt is datetime
hour = dt.hour hour = dt.hour
minute = dt.minute minute = dt.minute
second = dt.second second = dt.second
...@@ -590,7 +590,7 @@ def deg2HMS(ra0, dec0): ...@@ -590,7 +590,7 @@ def deg2HMS(ra0, dec0):
def distortField(ra, dec, ch): def distortField(ra, dec, ch):
# % ra ,dec are the idea position in arcsec , and the center position is ra=0, dec=0 # % ra ,dec are the idea position in arcsec , and the center position is ra=0, dec=0
'''MCI geometry distortion effect in channel G R and I''' '''MCI geometry distortion effect in channel G R and I'''
distortA = dict() distortA = dict()
...@@ -739,9 +739,9 @@ def centroid(data): ...@@ -739,9 +739,9 @@ def centroid(data):
Returns Returns
------- -------
cx: the centroid column number, in horizontal direction definet in python image show cx: the centroid column number, in horizontal direction definet in python image show
cy: the centroid row number , in vertical direction cy: the centroid row number , in vertical direction
''' '''
### ###
from scipy import ndimage from scipy import ndimage
...@@ -758,7 +758,7 @@ PSF interpolation for MCI_sim ...@@ -758,7 +758,7 @@ PSF interpolation for MCI_sim
''' '''
###find neighbors-KDtree ### # find neighbors-KDtree ###
def findNeighbors(tx, ty, px, py, dn=5): def findNeighbors(tx, ty, px, py, dn=5):
""" """
...@@ -782,7 +782,7 @@ def findNeighbors(tx, ty, px, py, dn=5): ...@@ -782,7 +782,7 @@ def findNeighbors(tx, ty, px, py, dn=5):
return indexq return indexq
############################################################################### ###############################################################################
###PSF-IDW### # PSF-IDW###
def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, dn=5, IDWindex=3, OnlyNeighbors=True): def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, dn=5, IDWindex=3, OnlyNeighbors=True):
...@@ -807,10 +807,10 @@ def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, dn=5, IDWindex=3, OnlyNeighbo ...@@ -807,10 +807,10 @@ def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, dn=5, IDWindex=3, OnlyNeighbo
ngy, ngx = PSFMat[:, :, 0, 0].shape # PSF data size ngy, ngx = PSFMat[:, :, 0, 0].shape # PSF data size
#################################### ####################################
####### my code ###### # my code ######
psfWeight = np.zeros([dn]) psfWeight = np.zeros([dn])
if OnlyNeighbors == True: if OnlyNeighbors:
neigh = findNeighbors(px, py, cen_col, cen_row, dn) neigh = findNeighbors(px, py, cen_col, cen_row, dn)
...@@ -929,10 +929,10 @@ class MCIsimulator(): ...@@ -929,10 +929,10 @@ class MCIsimulator():
# self.config.readfp(open(self.configfile)) # self.config.readfp(open(self.configfile))
self.config.read_file(open(self.configfile)) self.config.read_file(open(self.configfile))
#######
########################################################################### ###########################################################################
#######
#######
def processConfigs(self): def processConfigs(self):
""" """
Processes configuration information and save the information to a dictionary self.information. Processes configuration information and save the information to a dictionary self.information.
...@@ -998,7 +998,7 @@ class MCIsimulator(): ...@@ -998,7 +998,7 @@ class MCIsimulator():
self.section, 'save_cosmicrays') self.section, 'save_cosmicrays')
self.lensing = self.config.getboolean(self.section, 'lensing') self.lensing = self.config.getboolean(self.section, 'lensing')
###############################################3################################# # #################################
self.booleans = dict(cosmicRays=self.cosmicRays, self.booleans = dict(cosmicRays=self.cosmicRays,
darknoise=self.darknoise, darknoise=self.darknoise,
cosmetics=self.cosmetics, cosmetics=self.cosmetics,
...@@ -1026,17 +1026,17 @@ class MCIsimulator(): ...@@ -1026,17 +1026,17 @@ class MCIsimulator():
return return
############################################################################## # ##############################################################################
############################################################################### # #############################################################################
def _createEmpty(self): def _createEmpty(self):
""" """
Creates and empty array of a given x and y size full of zeros. Creates and empty array of a given x and y size full of zeros.
add g r i channel images; add g r i channel images;
Creates lensing parameters; Creates lensing parameters;
""" """
self.image_g = np.zeros( self.image_g = np.zeros(
(self.information['ysize'], self.information['xsize']), dtype=float) (self.information['ysize'], self.information['xsize']), dtype=float)
self.image_r = np.zeros( self.image_r = np.zeros(
...@@ -1062,7 +1062,7 @@ class MCIsimulator(): ...@@ -1062,7 +1062,7 @@ class MCIsimulator():
self.ghostMax = self.information['ghostratio'] self.ghostMax = self.information['ghostratio']
self.log.info('Maximum in the ghost model %e' % self.ghostMax) self.log.info('Maximum in the ghost model %e' % self.ghostMax)
return return
############################################### # ##############################################
def readCosmicRayInformation(self): def readCosmicRayInformation(self):
""" """
...@@ -1103,7 +1103,7 @@ class MCIsimulator(): ...@@ -1103,7 +1103,7 @@ class MCIsimulator():
self.source = sourcein self.source = sourcein
##print('print information:', self.information) # print('print information:', self.information)
########################################################################### ###########################################################################
...@@ -1128,7 +1128,7 @@ class MCIsimulator(): ...@@ -1128,7 +1128,7 @@ class MCIsimulator():
self.result_path = self.information['result_path'] + \ self.result_path = self.information['result_path'] + \
self.source+"_"+result_day self.source+"_"+result_day
if os.path.isdir(self.result_path) == False: if os.path.isdir(self.result_path) is False:
os.mkdir(self.result_path) os.mkdir(self.result_path)
os.mkdir(self.result_path+'/cali_Data') os.mkdir(self.result_path+'/cali_Data')
os.mkdir(self.result_path+'/log_Data') os.mkdir(self.result_path+'/log_Data')
...@@ -1207,7 +1207,7 @@ class MCIsimulator(): ...@@ -1207,7 +1207,7 @@ class MCIsimulator():
Parameters: Parameters:
px,py : target position px,py : target position
filternmae : the selected filter name filternmae : the selected filter name
ch : the MCI channle , g or r or i ch : the MCI channle , g or r or i
Returns: Returns:
PSF: PSF array at given 7 wavelength and at target position PSF: PSF array at given 7 wavelength and at target position
""" """
...@@ -1305,8 +1305,8 @@ class MCIsimulator(): ...@@ -1305,8 +1305,8 @@ class MCIsimulator():
############################################################################### ###############################################################################
def cal_sky_noise(self): def cal_sky_noise(self):
####### add earthshine ####### # ###### add earthshine #######
####### add earthshine ####### # ###### add earthshine #######
# self.earthshine_wave # A # self.earthshine_wave # A
# self.earthshine_flux # erg/s/cm^2/A/arcsec^2 # self.earthshine_flux # erg/s/cm^2/A/arcsec^2
...@@ -1383,7 +1383,7 @@ class MCIsimulator(): ...@@ -1383,7 +1383,7 @@ class MCIsimulator():
theta = rotTelPos - rotSkyPos theta = rotTelPos - rotSkyPos
############ load star data catlog ##################### # ########### load star data catlog #####################
# 'selection_MACSJ0744.9 3927_20230517_concat.csv' # 'selection_MACSJ0744.9 3927_20230517_concat.csv'
starcat = self.information['starcat'] starcat = self.information['starcat']
...@@ -1440,7 +1440,7 @@ class MCIsimulator(): ...@@ -1440,7 +1440,7 @@ class MCIsimulator():
channel = ['g', 'r', 'i'] channel = ['g', 'r', 'i']
####################################################################### #######################################################################
#################### cal_PSF_array ################################# # ################### cal_PSF_array #################################
if self.save_starpsf: if self.save_starpsf:
...@@ -1476,7 +1476,7 @@ class MCIsimulator(): ...@@ -1476,7 +1476,7 @@ class MCIsimulator():
fsx, fsy = cal_pos(center_ra, center_dec, fsx, fsy = cal_pos(center_ra, center_dec,
rotTelPos, rotSkyPos, galRa, galDec) rotTelPos, rotSkyPos, galRa, galDec)
############# do field distottion ########## # ############ do field distottion ##########
if self.distortion: if self.distortion:
...@@ -1556,8 +1556,8 @@ class MCIsimulator(): ...@@ -1556,8 +1556,8 @@ class MCIsimulator():
del PSF_i del PSF_i
self.log.info('finish save star PSF data......') self.log.info('finish save star PSF data......')
########## finish save PSF fits ########################################## # ######### finish save PSF fits ##########################################
############################################################################ # ###########################################################################
self.log.info('load star catlog successfully') self.log.info('load star catlog successfully')
...@@ -1720,7 +1720,7 @@ class MCIsimulator(): ...@@ -1720,7 +1720,7 @@ class MCIsimulator():
st_magi = [] st_magi = []
st_magz = [] st_magz = []
#################### generate star image ########## # ################### generate star image ##########
for j in tqdm(range(nsrcs)): # nsrcs length for j in tqdm(range(nsrcs)): # nsrcs length
...@@ -1728,7 +1728,7 @@ class MCIsimulator(): ...@@ -1728,7 +1728,7 @@ class MCIsimulator():
self.log.info('begin iteration........') self.log.info('begin iteration........')
# 33 # 33
# if self.debug: # if self.debug:
#self.log.info('j = %i' % (j)) # self.log.info('j = %i' % (j))
starRa = 3600.0*(newRa[j]) # ra of star, arcsecond starRa = 3600.0*(newRa[j]) # ra of star, arcsecond
starDec = 3600.0*(newDec[j]) # dec of star, arcsecond starDec = 3600.0*(newDec[j]) # dec of star, arcsecond
...@@ -1760,7 +1760,7 @@ class MCIsimulator(): ...@@ -1760,7 +1760,7 @@ class MCIsimulator():
star_input[ch][int(nlayccd)-1, 0] = fsx # ra star_input[ch][int(nlayccd)-1, 0] = fsx # ra
star_input[ch][int(nlayccd)-1, 1] = fsy # dec star_input[ch][int(nlayccd)-1, 1] = fsy # dec
############# do field distottion ########## # ############ do field distottion ##########
if self.distortion: if self.distortion:
...@@ -1791,7 +1791,7 @@ class MCIsimulator(): ...@@ -1791,7 +1791,7 @@ class MCIsimulator():
####################################################################### #######################################################################
####### use SED_code to generate star SED ####### # ###### use SED_code to generate star SED #######
if get_file_extension(starcat) == '.fits': if get_file_extension(starcat) == '.fits':
umag = self.star['gmag'][j] umag = self.star['gmag'][j]
gmag = self.star['gmag'][j] gmag = self.star['gmag'][j]
...@@ -1810,7 +1810,7 @@ class MCIsimulator(): ...@@ -1810,7 +1810,7 @@ class MCIsimulator():
flag = (umag > 0 and gmag > 0 and rmag > flag = (umag > 0 and gmag > 0 and rmag >
0 and imag > 0 and zmag > 0) 0 and imag > 0 and zmag > 0)
if flag == False: if flag is False:
continue continue
wave = np.linspace(2500, 11000, 8501) wave = np.linspace(2500, 11000, 8501)
...@@ -1832,7 +1832,7 @@ class MCIsimulator(): ...@@ -1832,7 +1832,7 @@ class MCIsimulator():
gx = dict() gx = dict()
gy = dict() gy = dict()
#self.log.info('magu =%f, magg=%f, magr=%f, magi=%f,magz=%f' % (umag,gmag,rmag,imag,zmag)) # self.log.info('magu =%f, magg=%f, magr=%f, magi=%f,magz=%f' % (umag,gmag,rmag,imag,zmag))
for i in range(3): for i in range(3):
ch = channel[i] ch = channel[i]
gy[ch] = obj[ch][1]/pixelscale+self.information['ysize'] / \ gy[ch] = obj[ch][1]/pixelscale+self.information['ysize'] / \
...@@ -1840,7 +1840,7 @@ class MCIsimulator(): ...@@ -1840,7 +1840,7 @@ class MCIsimulator():
gx[ch] = obj[ch][0]/pixelscale+self.information['xsize'] / \ gx[ch] = obj[ch][0]/pixelscale+self.information['xsize'] / \
2-0.5 # col number on CCD image 2-0.5 # col number on CCD image
#self.log.info('Channel in =%s, PosX(in pixel)=%f, PosY(in pixel)=%f' % (channel[i], gx[ch], gy[ch])) # self.log.info('Channel in =%s, PosX(in pixel)=%f, PosY(in pixel)=%f' % (channel[i], gx[ch], gy[ch]))
if self.debug: if self.debug:
...@@ -1872,8 +1872,8 @@ class MCIsimulator(): ...@@ -1872,8 +1872,8 @@ class MCIsimulator():
st_magi.append(imag) st_magi.append(imag)
st_magz.append(zmag) st_magz.append(zmag)
## self.log.info('begin cal PSF') # # self.log.info('begin cal PSF')
###################################################################### # #####################################################################
psf = dict() psf = dict()
for i in range(3): for i in range(3):
...@@ -1913,7 +1913,7 @@ class MCIsimulator(): ...@@ -1913,7 +1913,7 @@ class MCIsimulator():
# self.log.info('stamp_img.array.max()= %i' % (stamp_img.array.max())) # self.log.info('stamp_img.array.max()= %i' % (stamp_img.array.max()))
###print('stamp_img.array.max(): ', stamp_img.array.max() ) # ##print('stamp_img.array.max(): ', stamp_img.array.max() )
################################################# #################################################
...@@ -1922,10 +1922,10 @@ class MCIsimulator(): ...@@ -1922,10 +1922,10 @@ class MCIsimulator():
cx0, cy0 = centroid(stamp_img.array) cx0, cy0 = centroid(stamp_img.array)
##self.log.info('appfat= %s'%(self.appFatt)) # #self.log.info('appfat= %s'%(self.appFatt))
if self.appFatt: if self.appFatt:
###self.log.info('begin appfat ........') # ##self.log.info('begin appfat ........')
# apply treering and bright fatter and diffusion; # apply treering and bright fatter and diffusion;
SimpleTreeRing = galsim.SiliconSensor().simple_treerings( SimpleTreeRing = galsim.SiliconSensor().simple_treerings(
...@@ -1974,7 +1974,7 @@ class MCIsimulator(): ...@@ -1974,7 +1974,7 @@ class MCIsimulator():
# print(icx,icy) # print(icx,icy)
############################################################### ###############################################################
#### add ghost image #### # ### add ghost image ####
if self.ghosts: if self.ghosts:
ghostphotons = galsim.PhotonArray( ghostphotons = galsim.PhotonArray(
len(photons.x), x=photons.x, y=photons.y, flux=photons.flux) len(photons.x), x=photons.x, y=photons.y, flux=photons.flux)
...@@ -2046,7 +2046,7 @@ class MCIsimulator(): ...@@ -2046,7 +2046,7 @@ class MCIsimulator():
# get solar position # get solar position
dt = datetime.fromisoformat(time) dt = datetime.fromisoformat(time)
###jd = julian.to_jd(dt, fmt='jd') # ##jd = julian.to_jd(dt, fmt='jd')
jd = time2jd(dt) jd = time2jd(dt)
t = Time(jd, format='jd', scale='utc') t = Time(jd, format='jd', scale='utc')
...@@ -2067,7 +2067,7 @@ class MCIsimulator(): ...@@ -2067,7 +2067,7 @@ class MCIsimulator():
# interpolated zodical surface brightness at 0.5 um # interpolated zodical surface brightness at 0.5 um
zodi = pd.read_csv( zodi = pd.read_csv(
self.information['dir_path']+'MCI_inputData/refs/zodi_map.dat', sep='\s+', header=None, comment='#') self.information['dir_path']+'MCI_inputData/refs/zodi_map.dat', sep='\\s+', header=None, comment='#')
beta_angle = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75]) 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, lamda_angle = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40, 45,
60, 75, 90, 105, 120, 135, 150, 165, 180]) 60, 75, 90, 105, 120, 135, 150, 165, 180])
...@@ -2077,7 +2077,7 @@ class MCIsimulator(): ...@@ -2077,7 +2077,7 @@ class MCIsimulator():
# read the zodical spectrum in the ecliptic # read the zodical spectrum in the ecliptic
cat_spec = pd.read_csv( cat_spec = pd.read_csv(
self.information['dir_path']+'MCI_inputData/refs/solar_spec.dat', sep='\s+', header=None, comment='#') self.information['dir_path']+'MCI_inputData/refs/solar_spec.dat', sep='\\s+', header=None, comment='#')
wave = cat_spec[0].values # A wave = cat_spec[0].values # A
spec0 = cat_spec[1].values # 10^-8 W m^�? sr^�? μm^�? spec0 = cat_spec[1].values # 10^-8 W m^�? sr^�? μm^�?
zodi_norm = 252 # 10^-8 W m^�? sr^�? μm^�? zodi_norm = 252 # 10^-8 W m^�? sr^�? μm^�?
...@@ -2115,7 +2115,7 @@ class MCIsimulator(): ...@@ -2115,7 +2115,7 @@ class MCIsimulator():
dec = -34.155241212932374 dec = -34.155241212932374
# halo_id = 229600100382 # halo_id = 229600100382
if self.sim_star == False: if self.sim_star is False:
self.information['star_ra'] = ra self.information['star_ra'] = ra
self.information['star_dec'] = dec self.information['star_dec'] = dec
...@@ -2138,7 +2138,7 @@ class MCIsimulator(): ...@@ -2138,7 +2138,7 @@ class MCIsimulator():
self.information['ra_pnt0'] = center_ra self.information['ra_pnt0'] = center_ra
self.information['dec_pnt0'] = center_dec self.information['dec_pnt0'] = center_dec
#################################################################### ####################################################################
################################################################################## # ###################################################################
time_jd = time2jd(self.dt) time_jd = time2jd(self.dt)
...@@ -2186,7 +2186,7 @@ class MCIsimulator(): ...@@ -2186,7 +2186,7 @@ class MCIsimulator():
################################################################################## ##################################################################################
self.cal_sky_noise() self.cal_sky_noise()
############ load galaxy data with SED ############################ # ########### load galaxy data with SED ############################
self.information['target'] = deg2HMS( self.information['target'] = deg2HMS(
self.information['star_ra'], self.information['star_dec']) self.information['star_ra'], self.information['star_dec'])
...@@ -2309,7 +2309,7 @@ class MCIsimulator(): ...@@ -2309,7 +2309,7 @@ class MCIsimulator():
srcs_cat = fits.open(filename) srcs_cat = fits.open(filename)
#### load galaxy SED fitsfile ### # ### load galaxy SED fitsfile ###
if self.lensing: if self.lensing:
...@@ -2326,7 +2326,7 @@ class MCIsimulator(): ...@@ -2326,7 +2326,7 @@ class MCIsimulator():
self.log.info('galaxy_Input SED path is: %s' % (filename)) self.log.info('galaxy_Input SED path is: %s' % (filename))
###### do Tiance Effect ### # ##### do Tiance Effect ###
ra_list = [] ra_list = []
dec_list = [] dec_list = []
...@@ -2359,7 +2359,7 @@ class MCIsimulator(): ...@@ -2359,7 +2359,7 @@ class MCIsimulator():
################################################################### ###################################################################
for k1 in tqdm(range(len(newRa))): for k1 in tqdm(range(len(newRa))):
#print('k1=', k1) # print('k1=', k1)
# *(srcs_cat[k1].header['new_ra'] -self.information['gal_ra'] +self.information['star_ra']) # ra of galaxies, arcsecond # *(srcs_cat[k1].header['new_ra'] -self.information['gal_ra'] +self.information['star_ra']) # ra of galaxies, arcsecond
galRa = newRa[k1]*3600 galRa = newRa[k1]*3600
...@@ -2368,7 +2368,7 @@ class MCIsimulator(): ...@@ -2368,7 +2368,7 @@ class MCIsimulator():
########################################################################### ###########################################################################
############################################################################ ############################################################################
### test code #### # ## test code ####
################################################################# #################################################################
fsx, fsy = cal_pos(center_ra, center_dec, fsx, fsy = cal_pos(center_ra, center_dec,
...@@ -2382,7 +2382,7 @@ class MCIsimulator(): ...@@ -2382,7 +2382,7 @@ class MCIsimulator():
gal_flux = srcs_sed[k1+1].data gal_flux = srcs_sed[k1+1].data
################################ ################################
### rotate the lensed_images_g ### # ## rotate the lensed_images_g ###
if abs(theta.deg) > 0: if abs(theta.deg) > 0:
# here we choose reshape=False, the rotated image will # here we choose reshape=False, the rotated image will
lensed_images_g = ndimage.rotate( lensed_images_g = ndimage.rotate(
...@@ -2401,7 +2401,7 @@ class MCIsimulator(): ...@@ -2401,7 +2401,7 @@ class MCIsimulator():
galaxy_input[ch][int(nlayccd)-1, 1] = fsy # dec galaxy_input[ch][int(nlayccd)-1, 1] = fsy # dec
# galaxy_input[ch][int(nlayccd)-1, 2] =mag[ch] # mag g # galaxy_input[ch][int(nlayccd)-1, 2] =mag[ch] # mag g
############# do field distottion ########## # ############ do field distottion ##########
if self.distortion: if self.distortion:
...@@ -2436,10 +2436,10 @@ class MCIsimulator(): ...@@ -2436,10 +2436,10 @@ class MCIsimulator():
# galaxy_output[ch][int(nlayccd)-1, 2] =mag[ch] # mag g # galaxy_output[ch][int(nlayccd)-1, 2] =mag[ch] # mag g
####################################################################### #######################################################################
#self.log.info('Galaxy number=%i, Ra(in arcsec)=%f, Dec(in arcsec)=%f' % (nlayccd, fpx, fpy)) # self.log.info('Galaxy number=%i, Ra(in arcsec)=%f, Dec(in arcsec)=%f' % (nlayccd, fpx, fpy))
#self.log.info('Galaxy number=%i, Mag_g=%f, Mag_r=%f,Mag_i=%f' % (nlayccd, mag['g'],mag['r'],mag['i'])) # self.log.info('Galaxy number=%i, Mag_g=%f, Mag_r=%f,Mag_i=%f' % (nlayccd, mag['g'],mag['r'],mag['i']))
#print('Flux sum =', gal_flux.sum()) # print('Flux sum =', gal_flux.sum())
# cal_SED_photons(self, flux_arr): # cal_SED_photons(self, flux_arr):
intscales, PSF_eff_weight = self.cal_SED_photons(gal_flux) intscales, PSF_eff_weight = self.cal_SED_photons(gal_flux)
...@@ -2449,7 +2449,7 @@ class MCIsimulator(): ...@@ -2449,7 +2449,7 @@ class MCIsimulator():
if intscales['g'].max()/ls_imgsum < 0.1 and intscales['r'].max()/ls_imgsum < 0.1 and intscales['i'].max()/ls_imgsum < 0.1: if intscales['g'].max()/ls_imgsum < 0.1 and intscales['r'].max()/ls_imgsum < 0.1 and intscales['i'].max()/ls_imgsum < 0.1:
######### #########
continue continue
######## test code ########## # ####### test code ##########
if self.debug: if self.debug:
if nlayccd > 1: if nlayccd > 1:
...@@ -2520,7 +2520,7 @@ class MCIsimulator(): ...@@ -2520,7 +2520,7 @@ class MCIsimulator():
gal_rescale.append(srcs_cat[k1+1].header['RESCALE']) gal_rescale.append(srcs_cat[k1+1].header['RESCALE'])
gal_szs.append(srcs_cat[k1+1].header['SZS']) gal_szs.append(srcs_cat[k1+1].header['SZS'])
#self.log.info('Galaxy number=%i, Mag_g=%f, Mag_r=%f,Mag_i=%f' % (nlayccd, mag['g'],mag['r'],mag['i'])) # self.log.info('Galaxy number=%i, Mag_g=%f, Mag_r=%f,Mag_i=%f' % (nlayccd, mag['g'],mag['r'],mag['i']))
if nlayccd % 1000 == 0: if nlayccd % 1000 == 0:
self.log.info('Galaxy number on CCD is = %i' % (nlayccd)) self.log.info('Galaxy number on CCD is = %i' % (nlayccd))
...@@ -2564,7 +2564,7 @@ class MCIsimulator(): ...@@ -2564,7 +2564,7 @@ class MCIsimulator():
stamp_img, max_flux=1.0) stamp_img, max_flux=1.0)
cx0, cy0 = centroid(stamp_img.array) cx0, cy0 = centroid(stamp_img.array)
########### apply fat effect to galaxy image ##### # ########## apply fat effect to galaxy image #####
if self.appFatt: if self.appFatt:
# apply treering and bright fatter and diffusion; # apply treering and bright fatter and diffusion;
...@@ -2663,12 +2663,12 @@ class MCIsimulator(): ...@@ -2663,12 +2663,12 @@ class MCIsimulator():
fullimg['i'] = final_image['i'].array fullimg['i'] = final_image['i'].array
return fullimg return fullimg
#############################################
#################################
############################################################################### ###############################################################################
###########
############################################################################### ###############################################################################
################
def generatePRNU(self, ave=1.0, sigma=0.01): def generatePRNU(self, ave=1.0, sigma=0.01):
""" """
...@@ -2744,8 +2744,8 @@ class MCIsimulator(): ...@@ -2744,8 +2744,8 @@ class MCIsimulator():
Because the pixel-to-pixel non-uniformity effect (i.e. multiplicative) flat fielding takes place Because the pixel-to-pixel non-uniformity effect (i.e. multiplicative) flat fielding takes place
before CTI and other effects, the flat field file must be the same size as the pixels that see before CTI and other effects, the flat field file must be the same size as the pixels that see
the sky. the sky."""
"""
# generate flatfiledfile, with # generate flatfiledfile, with
flatsigma = self.information['flatsigma'] flatsigma = self.information['flatsigma']
...@@ -2790,7 +2790,7 @@ class MCIsimulator(): ...@@ -2790,7 +2790,7 @@ class MCIsimulator():
self.information['coveringfraction']*self.information['exptime']/300.0, limit=None) self.information['coveringfraction']*self.information['exptime']/300.0, limit=None)
# paste the information # paste the information
#self.image += CCD_cr # self.image += CCD_cr
self.image_g += CCD_cr_g self.image_g += CCD_cr_g
self.image_r += CCD_cr_r self.image_r += CCD_cr_r
self.image_i += CCD_cr_i self.image_i += CCD_cr_i
...@@ -2806,12 +2806,12 @@ class MCIsimulator(): ...@@ -2806,12 +2806,12 @@ class MCIsimulator():
self.information['simnumber'])+'.fits', np.int32(CCD_cr_i), overwrite=True) self.information['simnumber'])+'.fits', np.int32(CCD_cr_i), overwrite=True)
# count the covering factor # count the covering factor
#area_cr = np.count_nonzero(self.cosmicMap) # area_cr = np.count_nonzero(self.cosmicMap)
area_cr_g = np.count_nonzero(CCD_cr_g) area_cr_g = np.count_nonzero(CCD_cr_g)
area_cr_r = np.count_nonzero(CCD_cr_r) area_cr_r = np.count_nonzero(CCD_cr_r)
area_cr_i = np.count_nonzero(CCD_cr_i) area_cr_i = np.count_nonzero(CCD_cr_i)
#self.log.info('The cosmic ray covering factor is %i pixels ' % area_cr) # self.log.info('The cosmic ray covering factor is %i pixels ' % area_cr)
self.log.info( self.log.info(
'The cosmic ray in C1 channel covering factor is %i pixels ' % area_cr_g) 'The cosmic ray in C1 channel covering factor is %i pixels ' % area_cr_g)
self.log.info( self.log.info(
...@@ -3009,7 +3009,7 @@ class MCIsimulator(): ...@@ -3009,7 +3009,7 @@ class MCIsimulator():
# cosmetics_i[yc,xc]=val # cosmetics_i[yc,xc]=val
self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val)) self.log.info('x=%i, y=%i, value=%f' % (xc, yc, val))
####### save cosmetics image code #### # ###### save cosmetics image code ####
# if self.simnumber<2: # if self.simnumber<2:
# #save cosmic ray image map # #save cosmic ray image map
...@@ -3171,7 +3171,7 @@ class MCIsimulator(): ...@@ -3171,7 +3171,7 @@ class MCIsimulator():
ymin = dict() ymin = dict()
ymax = dict() ymax = dict()
# calculate the non-over subimg of the final matrix of the image # calculate the non-over subimg of the final matrix of the image
for k in range(1, 17): for k in range(1, 17):
xmin[k] = 0 xmin[k] = 0
xmax[k] = 4616+self.overscan xmax[k] = 4616+self.overscan
...@@ -3192,7 +3192,7 @@ class MCIsimulator(): ...@@ -3192,7 +3192,7 @@ class MCIsimulator():
# end for # end for
return return
#################################################################################### ##############################################################################
def applyBias(self): def applyBias(self):
""" """
...@@ -3271,10 +3271,10 @@ class MCIsimulator(): ...@@ -3271,10 +3271,10 @@ class MCIsimulator():
hdulist_g.writeto(file_g, overwrite=True) hdulist_g.writeto(file_g, overwrite=True)
return return
###########
########
#################################################################################### ##############################################################################
######
def applyBleeding(self, img, direction='not_horizon'): def applyBleeding(self, img, direction='not_horizon'):
""" """
...@@ -3304,14 +3304,14 @@ class MCIsimulator(): ...@@ -3304,14 +3304,14 @@ class MCIsimulator():
overload = value - self.information['fullwellcapacity'] overload = value - self.information['fullwellcapacity']
if overload > 0.: if overload > 0.:
overload /= 2. overload /= 2.
#self.image[j, i] -= overload # self.image[j, i] -= overload
data[i, j] -= overload data[i, j] -= overload
sum += overload sum += overload
elif sum > 0.: elif sum > 0.:
if -overload > sum: if -overload > sum:
overload = -sum overload = -sum
#self.image[j, i] -= overload # self.image[j, i] -= overload
data[i, j] -= overload data[i, j] -= overload
sum += overload sum += overload
###################################################### ######################################################
...@@ -3323,14 +3323,14 @@ class MCIsimulator(): ...@@ -3323,14 +3323,14 @@ class MCIsimulator():
# second round - from top to bottom (bleeding was half'd already, so now full) # second round - from top to bottom (bleeding was half'd already, so now full)
overload = value - self.information['fullwellcapacity'] overload = value - self.information['fullwellcapacity']
if overload > 0.: if overload > 0.:
#self.image[-j-1, i] -= overload # self.image[-j-1, i] -= overload
data[i, -j-1] -= overload data[i, -j-1] -= overload
sum += overload sum += overload
elif sum > 0.: elif sum > 0.:
if -overload > sum: if -overload > sum:
overload = -sum overload = -sum
#self.image[-j-1, i] -= overload # self.image[-j-1, i] -= overload
data[i, -j-1,] -= overload data[i, -j-1,] -= overload
sum += overload sum += overload
...@@ -3344,14 +3344,14 @@ class MCIsimulator(): ...@@ -3344,14 +3344,14 @@ class MCIsimulator():
overload = value - self.information['fullwellcapacity'] overload = value - self.information['fullwellcapacity']
if overload > 0.: if overload > 0.:
overload /= 2. overload /= 2.
#self.image[j, i] -= overload # self.image[j, i] -= overload
data[j, i] -= overload data[j, i] -= overload
sum += overload sum += overload
elif sum > 0.: elif sum > 0.:
if -overload > sum: if -overload > sum:
overload = -sum overload = -sum
#self.image[j, i] -= overload # self.image[j, i] -= overload
data[j, i] -= overload data[j, i] -= overload
sum += overload sum += overload
################################ ################################
...@@ -3361,18 +3361,18 @@ class MCIsimulator(): ...@@ -3361,18 +3361,18 @@ class MCIsimulator():
# second round - from top to bottom (bleeding was half'd already, so now full) # second round - from top to bottom (bleeding was half'd already, so now full)
overload = value - self.information['fullwellcapacity'] overload = value - self.information['fullwellcapacity']
if overload > 0.: if overload > 0.:
#self.image[-j-1, i] -= overload # self.image[-j-1, i] -= overload
data[-j-1, i] -= overload data[-j-1, i] -= overload
sum += overload sum += overload
elif sum > 0.: elif sum > 0.:
if -overload > sum: if -overload > sum:
overload = -sum overload = -sum
#self.image[-j-1, i] -= overload # self.image[-j-1, i] -= overload
data[-j-1, i] -= overload data[-j-1, i] -= overload
sum += overload sum += overload
#####print('Applying column bleeding finished.......') # ####print('Applying column bleeding finished.......')
return data return data
############################################################################ ############################################################################
...@@ -3427,7 +3427,7 @@ class MCIsimulator(): ...@@ -3427,7 +3427,7 @@ class MCIsimulator():
''' '''
### xsize is column, ysize is row # ## xsize is column, ysize is row
# horizon direction , columnn # horizon direction , columnn
self.prescan = self.information['prescan'] self.prescan = self.information['prescan']
...@@ -3571,7 +3571,7 @@ class MCIsimulator(): ...@@ -3571,7 +3571,7 @@ class MCIsimulator():
# exposure end time is t2 ; # exposure end time is t2 ;
t2 = self.dt+timedelta(seconds=self.information['exptime']) t2 = self.dt+timedelta(seconds=self.information['exptime'])
exp_endtime = t2.strftime("%Y%m%d%H%M%S") exp_endtime = t2.strftime("%Y%m%d%H%M%S")
### data read time is the exposure end time ### # ## data read time is the exposure end time ###
# data read end time is t3, and this is the write file time; # data read end time is t3, and this is the write file time;
t3 = self.dt+timedelta(seconds=self.information['exptime'])+timedelta( t3 = self.dt+timedelta(seconds=self.information['exptime'])+timedelta(
seconds=self.information['readouttime']) seconds=self.information['readouttime'])
...@@ -3610,7 +3610,7 @@ class MCIsimulator(): ...@@ -3610,7 +3610,7 @@ class MCIsimulator():
ofd_g.header['FITSSWV'] = ( ofd_g.header['FITSSWV'] = (
'mci_sim_0.8.03', 'FITS creating software version') 'mci_sim_0.8.03', 'FITS creating software version')
############################################################################## ##############################################################################
######### Object information ######################################### # ######## Object information #########################################
ofd_g.header['OBJECT'] = ( ofd_g.header['OBJECT'] = (
self.information['name_obj'][:30], 'object name') self.information['name_obj'][:30], 'object name')
ofd_g.header['TARGET'] = ( ofd_g.header['TARGET'] = (
...@@ -3622,7 +3622,7 @@ class MCIsimulator(): ...@@ -3622,7 +3622,7 @@ class MCIsimulator():
ofd_g.header['DEC_OBJ'] = ( ofd_g.header['DEC_OBJ'] = (
float(self.information['dec_obj']), 'object Dec (deg)') float(self.information['dec_obj']), 'object Dec (deg)')
######## Telescope information ############### # ####### Telescope information ###############
ofd_g.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') ofd_g.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version')
ofd_g.header['DATE-OBS'] = (data_time[:21], ofd_g.header['DATE-OBS'] = (data_time[:21],
...@@ -3760,7 +3760,7 @@ class MCIsimulator(): ...@@ -3760,7 +3760,7 @@ class MCIsimulator():
ofd_g.header['CHECKSUM'] = (0, '0') ofd_g.header['CHECKSUM'] = (0, '0')
################################################ ################################################
########## finish header for 0 layer #################### # ######### finish header for 0 layer ####################
# new image HDU # new image HDU
...@@ -3769,7 +3769,7 @@ class MCIsimulator(): ...@@ -3769,7 +3769,7 @@ class MCIsimulator():
else: else:
hdu_g = fits.ImageHDU(data=np.uint16(self.image_g)) hdu_g = fits.ImageHDU(data=np.uint16(self.image_g))
# #
### MCI header for image layer ### # ## MCI header for image layer ###
# need review thie header # need review thie header
hdu_g.header['XTENSION'] = ('IMAGE', 'extension type') hdu_g.header['XTENSION'] = ('IMAGE', 'extension type')
hdu_g.header['BITPIX'] = (np.int16(16), 'array data type') hdu_g.header['BITPIX'] = (np.int16(16), 'array data type')
...@@ -3789,7 +3789,7 @@ class MCIsimulator(): ...@@ -3789,7 +3789,7 @@ class MCIsimulator():
hdu_g.header['BUNIT'] = ('ADU', 'physical unit of array values') hdu_g.header['BUNIT'] = ('ADU', 'physical unit of array values')
###### file information ############### # ##### file information ###############
temp = t3.utcnow() temp = t3.utcnow()
data_time = temp.strftime("%Y-%m-%dT%H:%M:%S.%f") data_time = temp.strftime("%Y-%m-%dT%H:%M:%S.%f")
...@@ -3804,7 +3804,7 @@ class MCIsimulator(): ...@@ -3804,7 +3804,7 @@ class MCIsimulator():
######################################## ########################################
#### instrument status and object information ###### # ### instrument status and object information ######
hdu_g.header['TELESCOP'] = ('CSST', 'telescope name') hdu_g.header['TELESCOP'] = ('CSST', 'telescope name')
hdu_g.header['INSTRUME'] = ('MCI', 'instrument name') hdu_g.header['INSTRUME'] = ('MCI', 'instrument name')
hdu_g.header['CHANNEL'] = ('C1', 'channel number') hdu_g.header['CHANNEL'] = ('C1', 'channel number')
...@@ -3828,9 +3828,9 @@ class MCIsimulator(): ...@@ -3828,9 +3828,9 @@ class MCIsimulator():
hdu_g.header['DEC_OBJ'] = ( hdu_g.header['DEC_OBJ'] = (
float(self.information['dec_obj']), 'object Dec (deg)') float(self.information['dec_obj']), 'object Dec (deg)')
##### detector and Filter information ##### # #### detector and Filter information #####
hdu_g.header['FILTER'] = (self.filter_g[:6], 'filter band') hdu_g.header['FILTER'] = (self.filter_g[:6], 'filter band')
##### Detector information #### # #### Detector information ####
hdu_g.header['DETSN'] = ( hdu_g.header['DETSN'] = (
'E2V-CCD-290-0000000', 'detector serial number') 'E2V-CCD-290-0000000', 'detector serial number')
hdu_g.header['DETNAME'] = ('uv', 'detector name') hdu_g.header['DETNAME'] = ('uv', 'detector name')
...@@ -3871,7 +3871,7 @@ class MCIsimulator(): ...@@ -3871,7 +3871,7 @@ class MCIsimulator():
hdu_g.header['BIN_X'] = (np.int16(1), 'bin number in X') hdu_g.header['BIN_X'] = (np.int16(1), 'bin number in X')
hdu_g.header['BIN_Y'] = (np.int16(1), 'bin number in Y') hdu_g.header['BIN_Y'] = (np.int16(1), 'bin number in Y')
#### World coordinate system information ##### # ### World coordinate system information #####
hdu_g.header['WCSAXES'] = (np.int16(2), 'number of WCS axes') hdu_g.header['WCSAXES'] = (np.int16(2), 'number of WCS axes')
hdu_g.header['CRPIX1'] = ( hdu_g.header['CRPIX1'] = (
...@@ -3896,7 +3896,7 @@ class MCIsimulator(): ...@@ -3896,7 +3896,7 @@ class MCIsimulator():
hdu_g.header['CD2_2'] = ( hdu_g.header['CD2_2'] = (
float(self.information['CD2_2']), 'transformation matrix element (deg/pix)') float(self.information['CD2_2']), 'transformation matrix element (deg/pix)')
###### Readout information ####### # ##### Readout information #######
hdu_g.header['GAINLVL'] = (np.float32( hdu_g.header['GAINLVL'] = (np.float32(
1.5), 'gain level (e-/ADU)') 1.5), 'gain level (e-/ADU)')
hdu_g.header['GAIN01'] = (np.float32( hdu_g.header['GAIN01'] = (np.float32(
...@@ -4001,7 +4001,7 @@ class MCIsimulator(): ...@@ -4001,7 +4001,7 @@ class MCIsimulator():
hdu_g.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)') hdu_g.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)')
###################################### ######################################
#### exposure and shutter information ##### # ### exposure and shutter information #####
hdu_g.header['EXPSTART'] = (np.float64( hdu_g.header['EXPSTART'] = (np.float64(
time2mjd(self.dt)), 'exposure start time (MJD)') time2mjd(self.dt)), 'exposure start time (MJD)')
hdu_g.header['EXPEND'] = (float(time2mjd(t2)), hdu_g.header['EXPEND'] = (float(time2mjd(t2)),
...@@ -4012,7 +4012,7 @@ class MCIsimulator(): ...@@ -4012,7 +4012,7 @@ class MCIsimulator():
self.information['exptime']), 'dark current time (s)') self.information['exptime']), 'dark current time (s)')
hdu_g.header['SHTSTAT'] = (bool(True), 'shutter status') hdu_g.header['SHTSTAT'] = (bool(True), 'shutter status')
###### satellite and its allitide information ############## # ##### satellite and its allitide information ##############
hdu_g.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') hdu_g.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version')
hdu_g.header['SATESWV'] = ( hdu_g.header['SATESWV'] = (
'softwave-1.0', 'satellite software version') 'softwave-1.0', 'satellite software version')
...@@ -4210,7 +4210,7 @@ class MCIsimulator(): ...@@ -4210,7 +4210,7 @@ class MCIsimulator():
ofd_r.header['EQUINOX'] = (float(2000.0), '') ofd_r.header['EQUINOX'] = (float(2000.0), '')
ofd_r.header['FITSSWV'] = ( ofd_r.header['FITSSWV'] = (
'mci_sim_0.8.03', 'FITS creating software version') 'mci_sim_0.8.03', 'FITS creating software version')
######### Object information ######################################### # ######## Object information #########################################
ofd_r.header['OBJECT'] = ( ofd_r.header['OBJECT'] = (
self.information['name_obj'][:30], 'object name') self.information['name_obj'][:30], 'object name')
ofd_r.header['TARGET'] = ( ofd_r.header['TARGET'] = (
...@@ -4222,7 +4222,7 @@ class MCIsimulator(): ...@@ -4222,7 +4222,7 @@ class MCIsimulator():
ofd_r.header['DEC_OBJ'] = ( ofd_r.header['DEC_OBJ'] = (
float(self.information['dec_obj']), 'object Dec (deg)') float(self.information['dec_obj']), 'object Dec (deg)')
######## Telescope information ############### # ####### Telescope information ###############
ofd_r.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') ofd_r.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version')
ofd_r.header['DATE-OBS'] = (data_time[:21], ofd_r.header['DATE-OBS'] = (data_time[:21],
...@@ -4312,15 +4312,15 @@ class MCIsimulator(): ...@@ -4312,15 +4312,15 @@ class MCIsimulator():
ofd_r.header['CHECKSUM'] = (0, '0') ofd_r.header['CHECKSUM'] = (0, '0')
################################################ ################################################
########## finish header for 0 layer #################### # ######### finish header for 0 layer ####################
########## finish header for 0 layer #################### # ######### finish header for 0 layer ####################
# finish header for 0 layer # finish header for 0 layer
if HeaderTest == 'yes': if HeaderTest == 'yes':
hdu_r = fits.ImageHDU() hdu_r = fits.ImageHDU()
else: else:
hdu_r = fits.ImageHDU(data=np.uint16(self.image_r)) hdu_r = fits.ImageHDU(data=np.uint16(self.image_r))
### MCI header for image layer ### # ## MCI header for image layer ###
# need review thie header # need review thie header
hdu_r.header['XTENSION'] = ('IMAGE', 'extension type') hdu_r.header['XTENSION'] = ('IMAGE', 'extension type')
hdu_r.header['BITPIX'] = (np.int16(16), 'array data type') hdu_r.header['BITPIX'] = (np.int16(16), 'array data type')
...@@ -4355,7 +4355,7 @@ class MCIsimulator(): ...@@ -4355,7 +4355,7 @@ class MCIsimulator():
######################################## ########################################
#### instrument status and object information ###### # ### instrument status and object information ######
hdu_r.header['TELESCOP'] = ('CSST', 'telescope name') hdu_r.header['TELESCOP'] = ('CSST', 'telescope name')
hdu_r.header['INSTRUME'] = ('MCI', 'instrument name') hdu_r.header['INSTRUME'] = ('MCI', 'instrument name')
hdu_r.header['CHANNEL'] = ('C2', 'channel number') hdu_r.header['CHANNEL'] = ('C2', 'channel number')
...@@ -4379,7 +4379,7 @@ class MCIsimulator(): ...@@ -4379,7 +4379,7 @@ class MCIsimulator():
hdu_r.header['DEC_OBJ'] = ( hdu_r.header['DEC_OBJ'] = (
float(self.information['dec_obj']), 'object Dec (deg)') float(self.information['dec_obj']), 'object Dec (deg)')
##### detector and Filter information ##### # #### detector and Filter information #####
hdu_r.header['FILTER'] = (self.filter_r[:6], 'filter band') hdu_r.header['FILTER'] = (self.filter_r[:6], 'filter band')
##### Detector information #### ##### Detector information ####
hdu_r.header['DETSN'] = ( hdu_r.header['DETSN'] = (
...@@ -4422,7 +4422,7 @@ class MCIsimulator(): ...@@ -4422,7 +4422,7 @@ class MCIsimulator():
hdu_r.header['BIN_X'] = (np.int16(1), 'bin number in X') hdu_r.header['BIN_X'] = (np.int16(1), 'bin number in X')
hdu_r.header['BIN_Y'] = (np.int16(1), 'bin number in Y') hdu_r.header['BIN_Y'] = (np.int16(1), 'bin number in Y')
#### World coordinate system information ##### # ### World coordinate system information #####
hdu_r.header['WCSAXES'] = (np.int16(2), 'number of WCS axes') hdu_r.header['WCSAXES'] = (np.int16(2), 'number of WCS axes')
hdu_r.header['CRPIX1'] = ( hdu_r.header['CRPIX1'] = (
...@@ -4447,7 +4447,7 @@ class MCIsimulator(): ...@@ -4447,7 +4447,7 @@ class MCIsimulator():
hdu_r.header['CD2_2'] = ( hdu_r.header['CD2_2'] = (
float(self.information['CD2_2']), 'transformation matrix element (deg/pix)') float(self.information['CD2_2']), 'transformation matrix element (deg/pix)')
###### Readout information ####### # ##### Readout information #######
hdu_r.header['GAINLVL'] = (np.float32( hdu_r.header['GAINLVL'] = (np.float32(
1.5), 'gain level (e-/ADU)') 1.5), 'gain level (e-/ADU)')
hdu_r.header['GAIN01'] = (np.float32( hdu_r.header['GAIN01'] = (np.float32(
...@@ -4552,7 +4552,7 @@ class MCIsimulator(): ...@@ -4552,7 +4552,7 @@ class MCIsimulator():
hdu_r.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)') hdu_r.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)')
###################################### ######################################
#### exposure and shutter information ##### # ### exposure and shutter information #####
hdu_r.header['EXPSTART'] = (np.float64( hdu_r.header['EXPSTART'] = (np.float64(
time2mjd(self.dt)), 'exposure start time (MJD)') time2mjd(self.dt)), 'exposure start time (MJD)')
hdu_r.header['EXPEND'] = (float(time2mjd(t2)), hdu_r.header['EXPEND'] = (float(time2mjd(t2)),
...@@ -4563,7 +4563,7 @@ class MCIsimulator(): ...@@ -4563,7 +4563,7 @@ class MCIsimulator():
self.information['exptime']), 'dark current time (s)') self.information['exptime']), 'dark current time (s)')
hdu_r.header['SHTSTAT'] = (bool(True), 'shutter status') hdu_r.header['SHTSTAT'] = (bool(True), 'shutter status')
###### satellite and its allitide information ############## # ##### satellite and its allitide information ##############
hdu_r.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') hdu_r.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version')
hdu_r.header['SATESWV'] = ( hdu_r.header['SATESWV'] = (
'softwave-1.0', 'satellite software version') 'softwave-1.0', 'satellite software version')
...@@ -4759,7 +4759,7 @@ class MCIsimulator(): ...@@ -4759,7 +4759,7 @@ class MCIsimulator():
ofd_i.header['FITSSWV'] = ( ofd_i.header['FITSSWV'] = (
'mci_sim_0.8.03', 'FITS creating software version') 'mci_sim_0.8.03', 'FITS creating software version')
############################################################################## ##############################################################################
######### Object information ######################################### # ######## Object information #########################################
ofd_i.header['OBJECT'] = ( ofd_i.header['OBJECT'] = (
self.information['name_obj'][:30], 'object name') self.information['name_obj'][:30], 'object name')
ofd_i.header['TARGET'] = ( ofd_i.header['TARGET'] = (
...@@ -4771,7 +4771,7 @@ class MCIsimulator(): ...@@ -4771,7 +4771,7 @@ class MCIsimulator():
ofd_i.header['DEC_OBJ'] = ( ofd_i.header['DEC_OBJ'] = (
float(self.information['dec_obj']), 'object Dec (deg)') float(self.information['dec_obj']), 'object Dec (deg)')
######## Telescope information ############### # ####### Telescope information ###############
ofd_i.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') ofd_i.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version')
ofd_i.header['DATE-OBS'] = (data_time[:21], ofd_i.header['DATE-OBS'] = (data_time[:21],
'observation date (yyyy-mm-ddThh:mm:ss.s)') 'observation date (yyyy-mm-ddThh:mm:ss.s)')
...@@ -4860,7 +4860,7 @@ class MCIsimulator(): ...@@ -4860,7 +4860,7 @@ class MCIsimulator():
ofd_i.header['CHECKSUM'] = (0, '0') ofd_i.header['CHECKSUM'] = (0, '0')
################################################ ################################################
########## finish header for 0 layer #################### # ######### finish header for 0 layer ####################
# finish header for 0 layer # finish header for 0 layer
...@@ -4871,7 +4871,7 @@ class MCIsimulator(): ...@@ -4871,7 +4871,7 @@ class MCIsimulator():
hdu_i = fits.ImageHDU(data=np.uint16(self.image_i)) hdu_i = fits.ImageHDU(data=np.uint16(self.image_i))
# need review thie header # need review thie header
### MCI header for image layer ### # ## MCI header for image layer ###
# need review thie header # need review thie header
hdu_i.header['XTENSION'] = ('IMAGE', 'extension type') hdu_i.header['XTENSION'] = ('IMAGE', 'extension type')
hdu_i.header['BITPIX'] = (np.int16(16), 'array data type') hdu_i.header['BITPIX'] = (np.int16(16), 'array data type')
...@@ -4906,7 +4906,7 @@ class MCIsimulator(): ...@@ -4906,7 +4906,7 @@ class MCIsimulator():
######################################## ########################################
#### instrument status and object information ###### # ### instrument status and object information ######
hdu_i.header['TELESCOP'] = ('CSST', 'telescope name') hdu_i.header['TELESCOP'] = ('CSST', 'telescope name')
hdu_i.header['INSTRUME'] = ('MCI', 'instrument name') hdu_i.header['INSTRUME'] = ('MCI', 'instrument name')
hdu_i.header['CHANNEL'] = ('C3', 'channel number') hdu_i.header['CHANNEL'] = ('C3', 'channel number')
...@@ -4930,9 +4930,9 @@ class MCIsimulator(): ...@@ -4930,9 +4930,9 @@ class MCIsimulator():
hdu_i.header['DEC_OBJ'] = ( hdu_i.header['DEC_OBJ'] = (
float(self.information['dec_obj']), 'object Dec (deg)') float(self.information['dec_obj']), 'object Dec (deg)')
##### detector and Filter information ##### # #### detector and Filter information #####
hdu_i.header['FILTER'] = (self.filter_i[:6], 'filter band') hdu_i.header['FILTER'] = (self.filter_i[:6], 'filter band')
##### Detector information #### # #### Detector information ####
hdu_i.header['DETSN'] = ( hdu_i.header['DETSN'] = (
'E2V-CCD-290-0000000', 'detector serial number') 'E2V-CCD-290-0000000', 'detector serial number')
hdu_i.header['DETNAME'] = ('opt-red', 'detector name') hdu_i.header['DETNAME'] = ('opt-red', 'detector name')
...@@ -4973,7 +4973,7 @@ class MCIsimulator(): ...@@ -4973,7 +4973,7 @@ class MCIsimulator():
hdu_i.header['BIN_X'] = (np.int16(1), 'bin number in X') hdu_i.header['BIN_X'] = (np.int16(1), 'bin number in X')
hdu_i.header['BIN_Y'] = (np.int16(1), 'bin number in Y') hdu_i.header['BIN_Y'] = (np.int16(1), 'bin number in Y')
#### World coordinate system information ##### # ### World coordinate system information #####
hdu_i.header['WCSAXES'] = (np.int16(2), 'number of WCS axes') hdu_i.header['WCSAXES'] = (np.int16(2), 'number of WCS axes')
hdu_i.header['CRPIX1'] = ( hdu_i.header['CRPIX1'] = (
...@@ -4998,7 +4998,7 @@ class MCIsimulator(): ...@@ -4998,7 +4998,7 @@ class MCIsimulator():
hdu_i.header['CD2_2'] = ( hdu_i.header['CD2_2'] = (
float(self.information['CD2_2']), 'transformation matrix element (deg/pix)') float(self.information['CD2_2']), 'transformation matrix element (deg/pix)')
###### Readout information ####### # ##### Readout information #######
hdu_i.header['GAINLVL'] = (np.float32( hdu_i.header['GAINLVL'] = (np.float32(
1.5), 'gain level (e-/ADU)') 1.5), 'gain level (e-/ADU)')
hdu_i.header['GAIN01'] = (np.float32( hdu_i.header['GAIN01'] = (np.float32(
...@@ -5103,7 +5103,7 @@ class MCIsimulator(): ...@@ -5103,7 +5103,7 @@ class MCIsimulator():
hdu_i.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)') hdu_i.header['ROSPEED'] = (np.float32(100), 'readout speed (MHz)')
###################################### ######################################
#### exposure and shutter information ##### # ### exposure and shutter information #####
hdu_i.header['EXPSTART'] = (np.float64( hdu_i.header['EXPSTART'] = (np.float64(
time2mjd(self.dt)), 'exposure start time (MJD)') time2mjd(self.dt)), 'exposure start time (MJD)')
hdu_i.header['EXPEND'] = (float(time2mjd(t2)), hdu_i.header['EXPEND'] = (float(time2mjd(t2)),
...@@ -5114,7 +5114,7 @@ class MCIsimulator(): ...@@ -5114,7 +5114,7 @@ class MCIsimulator():
self.information['exptime']), 'dark current time (s)') self.information['exptime']), 'dark current time (s)')
hdu_i.header['SHTSTAT'] = (bool(True), 'shutter status') hdu_i.header['SHTSTAT'] = (bool(True), 'shutter status')
###### satellite and its allitide information ############## # ##### satellite and its allitide information ##############
hdu_i.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version') hdu_i.header['REFFRAME'] = ('CSSTGSC-1.0', 'guiding catalog version')
hdu_i.header['SATESWV'] = ( hdu_i.header['SATESWV'] = (
'softwave-1.0', 'satellite software version') 'softwave-1.0', 'satellite software version')
...@@ -5373,7 +5373,7 @@ class MCIsimulator(): ...@@ -5373,7 +5373,7 @@ class MCIsimulator():
self.orbit_pars[self.orbit_exp_num, 6]) self.orbit_pars[self.orbit_exp_num, 6])
########################################## ##########################################
###self.dt=julian.from_jd(exptime_start_jd, fmt='jd') # ##self.dt=julian.from_jd(exptime_start_jd, fmt='jd')
self.dt = jd2time(exptime_start_jd) self.dt = jd2time(exptime_start_jd)
# str(self.dt.year)+'-'+str(self.dt.month)+'-'+str(self.dt.day) # str(self.dt.year)+'-'+str(self.dt.month)+'-'+str(self.dt.day)
...@@ -5387,7 +5387,7 @@ class MCIsimulator(): ...@@ -5387,7 +5387,7 @@ class MCIsimulator():
# self.earthshine_theta=30.0 # in degree # self.earthshine_theta=30.0 # in degree
################################### ###################################
if self.sky_shift_rot: if self.sky_shift_rot:
if self.information['simnumber'] > 1: if self.information['simnumber'] > 1:
...@@ -5448,14 +5448,14 @@ class MCIsimulator(): ...@@ -5448,14 +5448,14 @@ class MCIsimulator():
self.information['rotSkyPos'] self.information['rotSkyPos']
self.information['rotangle'] = theta.deg self.information['rotangle'] = theta.deg
#print('rotangle',theta.deg ) # print('rotangle',theta.deg )
self.log.info('dis_ra(in pixel)=%f, dis_dec(in pixel)=%f, sky_rot(in deg)=%f' % ( self.log.info('dis_ra(in pixel)=%f, dis_dec(in pixel)=%f, sky_rot(in deg)=%f' % (
dis_ra*3600/0.05, dis_dec*3600/0.05, theta.deg)) dis_ra*3600/0.05, dis_dec*3600/0.05, theta.deg))
####################################################################################################### ###################################################################
############################################################################# ##################################################################
######### simulate star images from star SED data ################## # ######## simulate star images from star SED data ###########
if self.sim_star: if self.sim_star:
########################### ###########################
...@@ -5659,8 +5659,8 @@ class MCIsimulator(): ...@@ -5659,8 +5659,8 @@ class MCIsimulator():
print('Souce is not correct and programe will return') print('Souce is not correct and programe will return')
sys.exit(1) sys.exit(1)
######################################################################################################################## ###########################################################################
######################################################################################################################## #########################################################################
# Apply flat-field large scale structure for one chip # Apply flat-field large scale structure for one chip
if self.flatfieldM: if self.flatfieldM:
...@@ -5727,7 +5727,7 @@ class MCIsimulator(): ...@@ -5727,7 +5727,7 @@ class MCIsimulator():
self.applyRadiationDamage() self.applyRadiationDamage()
print('applyRadiationDamage()') print('applyRadiationDamage()')
###################################################################### ######################################################################
###### apply readoutNoise ###### # ##### apply readoutNoise ######
if self.readoutNoise: if self.readoutNoise:
self.applyReadoutNoise() self.applyReadoutNoise()
......
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