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

update

parent 0b24318d
Pipeline #7137 failed with stage
in 0 seconds
Show whitespace changes
Inline Side-by-side
csst_ifs_sim/csst_ifs_sim.py
View file @ 5c4a3a79
...@@ -2407,8 +2407,8 @@ class IFSsimulator(): ...@@ -2407,8 +2407,8 @@ class IFSsimulator():
########################################################## ##########################################################
######################################################################### #########################################################################
#
#
def applyDarkCurrent(self): def applyDarkCurrent(self):
""" """
Returns Returns
...@@ -2507,13 +2507,12 @@ class IFSsimulator(): ...@@ -2507,13 +2507,12 @@ class IFSsimulator():
self.information['dark4_r'] self.information['dark4_r']
# #############################################################################
# ############################################################################# # #############################################################################
def applyPoissonNoise(self):
"""
def applyPoissonNoise(self):
"""
Returns Returns
------- -------
None. None.
...@@ -2542,8 +2541,8 @@ class IFSsimulator(): ...@@ -2542,8 +2541,8 @@ class IFSsimulator():
self.log.info('Added Poisson noise on channel red') self.log.info('Added Poisson noise on channel red')
self.image_r += residual self.image_r += residual
# ############################################################################
############################################################################## #
def applyCosmetics(self): def applyCosmetics(self):
""" """
...@@ -3096,7 +3095,8 @@ class IFSsimulator(): ...@@ -3096,7 +3095,8 @@ class IFSsimulator():
else: else:
self.log.error( self.log.error(
'Cannot include pre- and overscan because of an unknown quadrant!') 'Cannot include pre- and overscan because of an unknown quadrant!')
#
#
# ################################################################# # #################################################################
canvas = np.zeros((self.information['redsize'], canvas = np.zeros((self.information['redsize'],
...@@ -3324,14 +3324,11 @@ class IFSsimulator(): ...@@ -3324,14 +3324,11 @@ class IFSsimulator():
def CCDreadout(self): def CCDreadout(self):
""" """
Returns Returns
------- -------
None. None.
""" """
imgb = self.image_b.copy() imgb = self.image_b.copy()
prescan = int(self.information['prescan']) prescan = int(self.information['prescan'])
overscan = int(self.information['overscan']) overscan = int(self.information['overscan'])
...@@ -3385,7 +3382,9 @@ class IFSsimulator(): ...@@ -3385,7 +3382,9 @@ class IFSsimulator():
y1 = 2418+prescan y1 = 2418+prescan
y2 = y1+2048 y2 = y1+2048
temp[x1:x2, y1:y2] = imgb[0:1024, 0:2048] temp[x1:x2, y1:y2] = imgb[0:1024, 0:2048]
#
#
#
# ## third part, old OSE, down to up ################### # ## third part, old OSE, down to up ###################
x1 = 0 x1 = 0
x2 = x1+1024 x2 = x1+1024
...@@ -3957,7 +3956,7 @@ class IFSsimulator(): ...@@ -3957,7 +3956,7 @@ class IFSsimulator():
ofd_r.header['FITSSWV'] = ( ofd_r.header['FITSSWV'] = (
'ifs_sim_0.8.03', 'FITS creating software version') 'ifs_sim_0.8.03', 'FITS creating software version')
######### Object information ############# # ######## Object information #############
if self.source == 'SCI' or self.source == 'COMP': if self.source == 'SCI' or self.source == 'COMP':
ofd_r.header['OBJECT'] = ( ofd_r.header['OBJECT'] = (
self.information['name_obj'][:30], 'object name') self.information['name_obj'][:30], 'object name')
...@@ -4362,7 +4361,7 @@ class IFSsimulator(): ...@@ -4362,7 +4361,7 @@ class IFSsimulator():
'========================================================================', before='EXPT0_1') '========================================================================', before='EXPT0_1')
hdulist_b.writeto(self.file_b, output_verify='ignore', checksum=True) hdulist_b.writeto(self.file_b, output_verify='ignore', checksum=True)
######################################################################### # ##################################################################
hdulist_r = fits.HDUList([hdu2, hdu_r]) hdulist_r = fits.HDUList([hdu2, hdu_r])
hdu_r.header.add_comment( hdu_r.header.add_comment(
...@@ -4399,7 +4398,7 @@ class IFSsimulator(): ...@@ -4399,7 +4398,7 @@ class IFSsimulator():
""" """
# read solar template # 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(self.information['dir_path']+'IFS_inputdata/refs/solar_spec.dat', sep='\\s+',
header=None, comment='#') header=None, comment='#')
template_wave = solar_template[0].values template_wave = solar_template[0].values
template_flux = solar_template[1].values template_flux = solar_template[1].values
...@@ -4411,7 +4410,7 @@ class IFSsimulator(): ...@@ -4411,7 +4410,7 @@ class IFSsimulator():
surface_brightness = earthshine_curve[1].values surface_brightness = earthshine_curve[1].values
# read V-band throughtput # 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(self.information['dir_path']+'IFS_inputdata/refs/filter_Bessell_V.dat', sep='\\s+',
header=None, comment='#') header=None, comment='#')
filter_wave = cat_filter_V[0].values filter_wave = cat_filter_V[0].values
filter_response = cat_filter_V[1].values filter_response = cat_filter_V[1].values
...@@ -4510,7 +4509,7 @@ class IFSsimulator(): ...@@ -4510,7 +4509,7 @@ class IFSsimulator():
""" """
sizeout = 1000 # here image has the pixelscale =0.01 arcsec sizeout = 1000 # here image has the pixelscale =0.01 arcsec
################# create the slicer and slit file for high resolution 0.01arcsec simulation ################### # create the slicer and slit file for high resolution 0.01arcsec simulation ###################
maskSlice = dict() maskSlice = dict()
maskSlit = dict() maskSlit = dict()
...@@ -4526,15 +4525,16 @@ class IFSsimulator(): ...@@ -4526,15 +4525,16 @@ class IFSsimulator():
self.maskSliceHole = maskSlice self.maskSliceHole = maskSlice
self.maskSlitHole = maskSlit self.maskSlitHole = maskSlit
######################################################################### ################################################################
##### load slicer_Qe1K.fits #### # ########
# #### load slicer_Qe1K.fits ####
slicerfile = self.information['dir_path'] + \ slicerfile = self.information['dir_path'] + \
'IFS_inputdata/Flatfield/slicer_QE_1K_230625.fits' 'IFS_inputdata/Flatfield/slicer_QE_1K_230625.fits'
da = fits.open(slicerfile) da = fits.open(slicerfile)
self.log.info('hole Sim and slicer 1K QE file:%s' % (slicerfile)) self.log.info('hole Sim and slicer 1K QE file:%s' % (slicerfile))
slicer_Qe = da[0].data slicer_Qe = da[0].data
##### load hole mask ###### # #### load hole mask ######
# load hole mask matrix # load hole mask matrix
file = self.information['dir_path'] + \ file = self.information['dir_path'] + \
'IFS_inputdata/Hole/holemask_230612.fits' 'IFS_inputdata/Hole/holemask_230612.fits'
...@@ -4597,7 +4597,7 @@ class IFSsimulator(): ...@@ -4597,7 +4597,7 @@ class IFSsimulator():
Wave[ilam] = lam Wave[ilam] = lam
#print('i = ',ilam ) # print('i = ',ilam )
############################################### ###############################################
sumt = 0.0 sumt = 0.0
...@@ -4606,7 +4606,7 @@ class IFSsimulator(): ...@@ -4606,7 +4606,7 @@ class IFSsimulator():
t = getSpectrum(HgArst, lam+(k-5) * t = getSpectrum(HgArst, lam+(k-5) *
Width_lambda/11, self.HgArsigma) Width_lambda/11, self.HgArsigma)
sumt = sumt+t[0]/11 sumt = sumt+t[0]/11
#print('t=', t) # print('t=', t)
# the normalized HgAr light source in the range of 0 and 1.0 # the normalized HgAr light source in the range of 0 and 1.0
LightSt[ilam] = sumt LightSt[ilam] = sumt
...@@ -4688,13 +4688,13 @@ class IFSsimulator(): ...@@ -4688,13 +4688,13 @@ class IFSsimulator():
# consider the slice optical efficiency in different slicer channel # consider the slice optical efficiency in different slicer channel
img0 = img0*slicer_Qe img0 = img0*slicer_Qe
########## do the slice effect ################### # ######### do the slice effect ###################
for k in range(32): for k in range(32):
# do slice effect to get slice image # do slice effect to get slice image
img1 = img0*self.maskSliceHole[str(k)] img1 = img0*self.maskSliceHole[str(k)]
############ get opd2 and PSF2 ###################### # ########### get opd2 and PSF2 ######################
wavefront2 = self.opd2/lam wavefront2 = self.opd2/lam
psf2 = anySampledPSFnew(wavefront2, self.pupil, Q, 256) psf2 = anySampledPSFnew(wavefront2, self.pupil, Q, 256)
# do convolve # do convolve
...@@ -4704,11 +4704,11 @@ class IFSsimulator(): ...@@ -4704,11 +4704,11 @@ class IFSsimulator():
conv[conv < 0.0] = 0.0 conv[conv < 0.0] = 0.0
img2 = conv img2 = conv
##############do Slit Mask ########################### # #############do Slit Mask ###########################
img2 = img2*self.maskSlitHole[str(k)] img2 = img2*self.maskSlitHole[str(k)]
######### get opd3 and PSF3 ########################## # ######## get opd3 and PSF3 ##########################
wavefront3 = self.opd3/lam wavefront3 = self.opd3/lam
psf3 = anySampledPSFnew(wavefront3, self.pupil, Q, 256) psf3 = anySampledPSFnew(wavefront3, self.pupil, Q, 256)
...@@ -4721,14 +4721,14 @@ class IFSsimulator(): ...@@ -4721,14 +4721,14 @@ class IFSsimulator():
img3 = conv img3 = conv
######################## get subimag ##################### # ####################### get subimag #####################
subimage = galsim.Image(800, 800) subimage = galsim.Image(800, 800)
subimage.array[:, :] = img3[int( subimage.array[:, :] = img3[int(
sizeout/2)-400:int(sizeout/2)+400, int(sizeout/2)-400:int(sizeout/2)+400] sizeout/2)-400:int(sizeout/2)+400, int(sizeout/2)-400:int(sizeout/2)+400]
# fits.writeto('subimage.fits',subimage.array,overwrite=True) # fits.writeto('subimage.fits',subimage.array,overwrite=True)
####### get photons from sub-image ##################### # ###### get photons from sub-image #####################
subimage.scale = 0.01 # here pixelscale=0.01arcsec subimage.scale = 0.01 # here pixelscale=0.01arcsec
subimage.setOrigin(0, 0) subimage.setOrigin(0, 0)
...@@ -4777,7 +4777,7 @@ class IFSsimulator(): ...@@ -4777,7 +4777,7 @@ class IFSsimulator():
photons_red.y = photons_red.y / \ photons_red.y = photons_red.y / \
(self.pixelscale)+dy_red+self.slice_red['py'][k] (self.pixelscale)+dy_red+self.slice_red['py'][k]
#red_sensor.accumulate(photons_red, red_img) # red_sensor.accumulate(photons_red, red_img)
photons_red.addTo(red_img) photons_red.addTo(red_img)
...@@ -4801,7 +4801,7 @@ class IFSsimulator(): ...@@ -4801,7 +4801,7 @@ class IFSsimulator():
fits.writeto(bluefile, blue_img.array, overwrite=True) fits.writeto(bluefile, blue_img.array, overwrite=True)
fits.writeto(redfile, red_img.array, overwrite=True) fits.writeto(redfile, red_img.array, overwrite=True)
#### save flux #### # ### save flux ####
Flux = flux_1 Flux = flux_1
hdu1 = fits.PrimaryHDU(Flux) hdu1 = fits.PrimaryHDU(Flux)
...@@ -4822,7 +4822,7 @@ class IFSsimulator(): ...@@ -4822,7 +4822,7 @@ class IFSsimulator():
newd.writeto(filename, overwrite=True) newd.writeto(filename, overwrite=True)
return return
################################################################################## # ########################################################################
def sim_sky_img(self, skyfitsfilename, skyRa, skyDec, sky_rot, telRa, telDec, exposuretime, simnumber): def sim_sky_img(self, skyfitsfilename, skyRa, skyDec, sky_rot, telRa, telDec, exposuretime, simnumber):
""" """
...@@ -4885,7 +4885,7 @@ class IFSsimulator(): ...@@ -4885,7 +4885,7 @@ class IFSsimulator():
self.information['dec_pnt0'] = a[0].header['DEC'] + disDec self.information['dec_pnt0'] = a[0].header['DEC'] + disDec
############### calculate the earthshine and zodiacal noise ,new code 2023.11.1 ############ # calculate the earthshine and zodiacal noise ,new code 2023.11.1 ############
############### ###############
self.log.info('Real telescope pointing in Ra = %f, Dec = %f' % ( self.log.info('Real telescope pointing in Ra = %f, Dec = %f' % (
...@@ -4957,9 +4957,9 @@ class IFSsimulator(): ...@@ -4957,9 +4957,9 @@ class IFSsimulator():
# the main input data # the main input data
# the relatived wavelength which is converted from Unit A to nm # the relatived wavelength which is converted from Unit A to nm
Wave = 0.1*a[2].data Wave = 0.1*a[2].data
#print('Wave data header', hdr) # print('Wave data header', hdr)
###################################################################################### # #####################################################################
exptime = self.information['exptime'] # exposure time exptime = self.information['exptime'] # exposure time
...@@ -4989,7 +4989,7 @@ class IFSsimulator(): ...@@ -4989,7 +4989,7 @@ class IFSsimulator():
width_blue = 0 width_blue = 0
################################ ################################
############## doppler effect to photons.wavelength ############# # ############# doppler effect to photons.wavelength #############
if self.appbianpai: if self.appbianpai:
sn = self.simnumber-1 sn = self.simnumber-1
...@@ -5054,7 +5054,7 @@ class IFSsimulator(): ...@@ -5054,7 +5054,7 @@ class IFSsimulator():
# 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'])
### data read time is the exposure end time plus readouttime ### # ## data read time is the exposure end time plus readouttime ###
# data read end time is t3 # data read end time is t3
t2jd = time2jd(t2) t2jd = time2jd(t2)
...@@ -5185,7 +5185,7 @@ class IFSsimulator(): ...@@ -5185,7 +5185,7 @@ class IFSsimulator():
# if ilam==100: # if ilam==100:
# fits.writeto('/media/yan//IFSsim/IFSsim Data/original_Img.fits', Nimg, overwrite=True) # fits.writeto('/media/yan//IFSsim/IFSsim Data/original_Img.fits', Nimg, overwrite=True)
########################################################################## # #########################################################################
# shift image with photons position and rotate them round the image true center # shift image with photons position and rotate them round the image true center
img = galsim.Image(Nimg, copy=True) img = galsim.Image(Nimg, copy=True)
img.scale = self.pixelscale img.scale = self.pixelscale
...@@ -5206,7 +5206,7 @@ class IFSsimulator(): ...@@ -5206,7 +5206,7 @@ class IFSsimulator():
x2, y2 = rotation_yan(img.true_center.x, img.true_center.y, x2, y2 = rotation_yan(img.true_center.x, img.true_center.y,
photons.x/img.scale, photons.y/img.scale, sky_rot) photons.x/img.scale, photons.y/img.scale, sky_rot)
#print('rotation time=', t2-t1) # print('rotation time=', t2-t1)
photons2 = galsim.PhotonArray( photons2 = galsim.PhotonArray(
N=len(x2), x=x2*img.scale, y=y2*img.scale, flux=photons.flux) N=len(x2), x=x2*img.scale, y=y2*img.scale, flux=photons.flux)
...@@ -5217,7 +5217,7 @@ class IFSsimulator(): ...@@ -5217,7 +5217,7 @@ class IFSsimulator():
image0 = galsim.Image(sizeout, sizeout) image0 = galsim.Image(sizeout, sizeout)
image0.scale = self.pixelscale image0.scale = self.pixelscale
image0.setOrigin(0, 0) image0.setOrigin(0, 0)
#image.setCenter(int(np.mean(photons.x)), int(np.mean(photons.y))) # image.setCenter(int(np.mean(photons.x)), int(np.mean(photons.y)))
rotphotons = galsim.PhotonArray( rotphotons = galsim.PhotonArray(
N=len(photons.x), x=x2, y=y2, flux=photons2.flux) N=len(photons.x), x=x2, y=y2, flux=photons2.flux)
...@@ -5238,16 +5238,16 @@ class IFSsimulator(): ...@@ -5238,16 +5238,16 @@ class IFSsimulator():
hdu1.header['rotAngle'] = ( hdu1.header['rotAngle'] = (
sky_rot, 'sky rotation relative to telescope anticlockwise in degree') sky_rot, 'sky rotation relative to telescope anticlockwise in degree')
#PSFfilename='/media/yan//IFSsim/IFSsim Data/rot_shift_Img.fits' # 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 ) # 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_' + hdu1.writeto(self.result_path+'/shift_rot_sky/rot_shift_subImg_' +
str(simnumber)+'.fits', overwrite=True) str(simnumber)+'.fits', overwrite=True)
fits.writeto(self.result_path+'/shift_rot_sky/original_Img_' + fits.writeto(self.result_path+'/shift_rot_sky/original_Img_' +
str(simnumber)+'.fits', Nimg, overwrite=True) str(simnumber)+'.fits', Nimg, overwrite=True)
##################################################################### #####################################################################
### do convolve image0 with PSF0 from primay CSST ### # ## do convolve image0 with PSF0 from primay CSST ###
# calculate the PSF0 at this wavelength # calculate the PSF0 at this wavelength
Q = lam*1e-3*self.Fnum/self.pixelsize Q = lam*1e-3*self.Fnum/self.pixelsize
...@@ -5258,7 +5258,7 @@ class IFSsimulator(): ...@@ -5258,7 +5258,7 @@ class IFSsimulator():
psf0 = anySampledPSFnew(wavefront, self.pupil, Q, 64) psf0 = anySampledPSFnew(wavefront, self.pupil, Q, 64)
#psf01=anySampledPSFnew(wavefront, self.pupil , Q , 64) # psf01=anySampledPSFnew(wavefront, self.pupil , Q , 64)
conv = fftconvolve(image0.array, psf0, mode='same') conv = fftconvolve(image0.array, psf0, mode='same')
conv[conv < 0.0] = 0.0 conv[conv < 0.0] = 0.0
...@@ -5301,13 +5301,13 @@ class IFSsimulator(): ...@@ -5301,13 +5301,13 @@ class IFSsimulator():
lam2 = lam*(1+vv2/(3.0*1e5)) lam2 = lam*(1+vv2/(3.0*1e5))
img0 = img0*slicer_Qe img0 = img0*slicer_Qe
########## do the slice effect ################### # ######### do the slice effect ###################
for k in range(32): for k in range(32):
# do slice effect to get slice image # do slice effect to get slice image
img1 = img0*self.maskSlice[str(k)] img1 = img0*self.maskSlice[str(k)]
############ get opd2 and PSF2 ###################### # ########### get opd2 and PSF2 ######################
wavefront2 = self.opd2/lam wavefront2 = self.opd2/lam
psf2 = anySampledPSFnew(wavefront2, self.pupil, Q, 64) psf2 = anySampledPSFnew(wavefront2, self.pupil, Q, 64)
# here we choose reshape=False, the rotated image will # here we choose reshape=False, the rotated image will
...@@ -5319,11 +5319,11 @@ class IFSsimulator(): ...@@ -5319,11 +5319,11 @@ class IFSsimulator():
conv[conv < 0.0] = 0.0 conv[conv < 0.0] = 0.0
img2 = conv img2 = conv
##############do Slit Mask ########################### # #############do Slit Mask ###########################
img2 = img2*self.maskSlit[str(k)] img2 = img2*self.maskSlit[str(k)]
######### get opd3 and PSF3 ########################## # ######## get opd3 and PSF3 ##########################
wavefront3 = self.opd3/lam wavefront3 = self.opd3/lam
psf3 = anySampledPSFnew(wavefront3, self.pupil, Q, 64) psf3 = anySampledPSFnew(wavefront3, self.pupil, Q, 64)
...@@ -5338,9 +5338,9 @@ class IFSsimulator(): ...@@ -5338,9 +5338,9 @@ class IFSsimulator():
img3 = conv img3 = conv
#print('covtime =', t2-t1) # print('covtime =', t2-t1)
######################## get subimage ##################### # ### get subimage #####################
subimage = galsim.Image(80, 80) subimage = galsim.Image(80, 80)
subimage.array[:, :] = img3[int( subimage.array[:, :] = img3[int(
sizeout/2)-40:int(sizeout/2)+40, int(sizeout/2)-40:int(sizeout/2)+40] sizeout/2)-40:int(sizeout/2)+40, int(sizeout/2)-40:int(sizeout/2)+40]
...@@ -5355,17 +5355,17 @@ class IFSsimulator(): ...@@ -5355,17 +5355,17 @@ class IFSsimulator():
photons = galsim.PhotonArray.makeFromImage( photons = galsim.PhotonArray.makeFromImage(
subimage, max_flux=max(img3.max()/10000.0, 1)) subimage, max_flux=max(img3.max()/10000.0, 1))
############################################################################# ############################################################
# do something for each photons; # do something for each photons;
############################################################################## #############################################################
idx0 = np.where(photons.flux > 1e-3) idx0 = np.where(photons.flux > 1e-3)
# totla energy for slice image # totla energy for slice image
energy = energy+sum(photons.flux[idx0]) energy = energy+sum(photons.flux[idx0])
p_num = len(idx0[0]) p_num = len(idx0[0])
############################################################################################### ##############################################################
####################### code for test slice is right or not ################################### # ######### 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) # photons_2=galsim.PhotonArray(N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux = photons.flux[idx0], wavelength=lam)
# if k==0: # if k==0:
# sliceimage = galsim.Image(100,100) # sliceimage = galsim.Image(100,100)
...@@ -5425,8 +5425,8 @@ class IFSsimulator(): ...@@ -5425,8 +5425,8 @@ class IFSsimulator():
######################################################## ########################################################
########################### test code ################################ # ##################### test code ################################
###################################################################### # ################################################################
# fits.writeto('originalImg.fits',Nimg,overwrite=True) # fits.writeto('originalImg.fits',Nimg,overwrite=True)
# fits.writeto('PSFedImg.fits',temp,overwrite=True) # fits.writeto('PSFedImg.fits',temp,overwrite=True)
# fits.writeto('rotedphotonsImg.fits',image.array,overwrite=True) # fits.writeto('rotedphotonsImg.fits',image.array,overwrite=True)
...@@ -5434,10 +5434,10 @@ class IFSsimulator(): ...@@ -5434,10 +5434,10 @@ class IFSsimulator():
# fits.writeto('slicephotonsImg.fits',sliceimage.array,overwrite=True) # fits.writeto('slicephotonsImg.fits',sliceimage.array,overwrite=True)
# pause() # pause()
###################################################################### # #############################################################
###################### finish test code ########################### # ############## finish test code ###########################
#################################################################################### #####################################################################
# stray light will cover 2% of input total light; # stray light will cover 2% of input total light;
if self.sky_noise: if self.sky_noise:
...@@ -5480,7 +5480,7 @@ class IFSsimulator(): ...@@ -5480,7 +5480,7 @@ class IFSsimulator():
""" """
sizeout = 100 sizeout = 100
################# load the hole mask file ################### # ################ load the hole mask file ###################
# load HgAr data; # load HgAr data;
if self.source == 'LAMP': if self.source == 'LAMP':
...@@ -5548,7 +5548,7 @@ class IFSsimulator(): ...@@ -5548,7 +5548,7 @@ class IFSsimulator():
Wave[ilam] = lam Wave[ilam] = lam
# tab['wavelength'][ilam]=lam # tab['wavelength'][ilam]=lam
#print('i = ',ilam ) # print('i = ',ilam )
############################################### ###############################################
if self.source == 'LAMP': if self.source == 'LAMP':
...@@ -5558,7 +5558,7 @@ class IFSsimulator(): ...@@ -5558,7 +5558,7 @@ class IFSsimulator():
t = getSpectrum(HgArst, lam+(k-5) * t = getSpectrum(HgArst, lam+(k-5) *
Width_lambda/11, self.HgArsigma) Width_lambda/11, self.HgArsigma)
sumt = sumt+t[0]/11 sumt = sumt+t[0]/11
#print('t=', t) # print('t=', t)
# the normalized HgAr light source in the range of 0 and 1.0 # the normalized HgAr light source in the range of 0 and 1.0
LightSt[ilam] = sumt LightSt[ilam] = sumt
...@@ -5622,7 +5622,7 @@ class IFSsimulator(): ...@@ -5622,7 +5622,7 @@ class IFSsimulator():
psf0 = anySampledPSFnew(wavefront0, self.pupil, Q, 64) psf0 = anySampledPSFnew(wavefront0, self.pupil, Q, 64)
#psf01=anySampledPSFnew(wavefront, self.pupil , Q , 64) # psf01=anySampledPSFnew(wavefront, self.pupil , Q , 64)
conv = fftconvolve(slice_image, psf0, mode='same') conv = fftconvolve(slice_image, psf0, mode='same')
...@@ -5663,13 +5663,13 @@ class IFSsimulator(): ...@@ -5663,13 +5663,13 @@ class IFSsimulator():
self.information['dir_path']+'IFS_inputdata/Flatfield/slicer_QE230625.fits') self.information['dir_path']+'IFS_inputdata/Flatfield/slicer_QE230625.fits')
slicer_Qe = da[0].data slicer_Qe = da[0].data
img0 = img0*slicer_Qe img0 = img0*slicer_Qe
########## do the slice effect ################### # ######### do the slice effect ###################
for k in range(32): for k in range(32):
# do slice effect to get slice image # do slice effect to get slice image
img1 = img0*self.maskSlice[str(k)] img1 = img0*self.maskSlice[str(k)]
############ get opd2 and PSF2 ###################### # ########### get opd2 and PSF2 ######################
wavefront2 = self.opd2/lam wavefront2 = self.opd2/lam
psf2 = anySampledPSFnew(wavefront2, self.pupil, Q, 64) psf2 = anySampledPSFnew(wavefront2, self.pupil, Q, 64)
# do convolve # do convolve
...@@ -5679,11 +5679,11 @@ class IFSsimulator(): ...@@ -5679,11 +5679,11 @@ class IFSsimulator():
conv[conv < 0.0] = 0.0 conv[conv < 0.0] = 0.0
img2 = conv img2 = conv
##############do Slit Mask ########################### # #############do Slit Mask ###########################
img2 = img2*self.maskSlit[str(k)] img2 = img2*self.maskSlit[str(k)]
######### get opd3 and PSF3 ########################## # ######## get opd3 and PSF3 ##########################
wavefront3 = self.opd3/lam wavefront3 = self.opd3/lam
psf3 = anySampledPSFnew(wavefront3, self.pupil, Q, 64) psf3 = anySampledPSFnew(wavefront3, self.pupil, Q, 64)
...@@ -5696,14 +5696,14 @@ class IFSsimulator(): ...@@ -5696,14 +5696,14 @@ class IFSsimulator():
img3 = conv img3 = conv
######################## get subimage ##################### # # ## get subimage #####################
subimage = galsim.Image(80, 80) subimage = galsim.Image(80, 80)
subimage.array[:, :] = img3[int( subimage.array[:, :] = img3[int(
sizeout/2)-40:int(sizeout/2)+40, int(sizeout/2)-40:int(sizeout/2)+40] sizeout/2)-40:int(sizeout/2)+40, int(sizeout/2)-40:int(sizeout/2)+40]
# fits.writeto('subimage.fits',subimage.array,overwrite=True) # fits.writeto('subimage.fits',subimage.array,overwrite=True)
######################## get photons from sub-image ##################### # ##### get photons from sub-image #####################
subimage.scale = self.pixelscale subimage.scale = self.pixelscale
subimage.setOrigin(0, 0) subimage.setOrigin(0, 0)
...@@ -5711,16 +5711,16 @@ class IFSsimulator(): ...@@ -5711,16 +5711,16 @@ class IFSsimulator():
photons = galsim.PhotonArray.makeFromImage( photons = galsim.PhotonArray.makeFromImage(
subimage, max_flux=max(img3.max()/10000.0, 1)) subimage, max_flux=max(img3.max()/10000.0, 1))
############################################################################# # ###########################################################
# do something for each photons; # do something for each photons;
############################################################################## # ############################################################
idx0 = np.where(photons.flux > 1e-2) idx0 = np.where(photons.flux > 1e-2)
# energy=energy+sum(photons.flux[idx0]) ### totla energy for slice image # energy=energy+sum(photons.flux[idx0]) ### totla energy for slice image
p_num = len(idx0[0]) p_num = len(idx0[0])
############################################################################################### # #############################################################
####################### code for test slice is right or not ################################### # ##### 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) # photons_2=galsim.PhotonArray(N=p_num, x=photons.x[idx0], y=photons.y[idx0], flux = photons.flux[idx0], wavelength=lam)
# if k==0: # if k==0:
# sliceimage = galsim.Image(100,100) # sliceimage = galsim.Image(100,100)
...@@ -5751,7 +5751,7 @@ class IFSsimulator(): ...@@ -5751,7 +5751,7 @@ class IFSsimulator():
photons_blue.y = photons_blue.y/self.pixelscale + \ photons_blue.y = photons_blue.y/self.pixelscale + \
dy_blue+self.slice_blue['py'][k] dy_blue+self.slice_blue['py'][k]
#photons_blue.y =2048-photons_blue.y # photons_blue.y =2048-photons_blue.y
blue_sensor.accumulate(photons_blue, blue_img) blue_sensor.accumulate(photons_blue, blue_img)
...@@ -5774,7 +5774,7 @@ class IFSsimulator(): ...@@ -5774,7 +5774,7 @@ class IFSsimulator():
photons_red.y = photons_red.y/self.pixelscale + \ photons_red.y = photons_red.y/self.pixelscale + \
dy_red+self.slice_red['py'][k] dy_red+self.slice_red['py'][k]
#photons_red.y =3072-photons_red.y # photons_red.y =3072-photons_red.y
red_sensor.accumulate(photons_red, red_img) red_sensor.accumulate(photons_red, red_img)
...@@ -5798,7 +5798,7 @@ class IFSsimulator(): ...@@ -5798,7 +5798,7 @@ class IFSsimulator():
fits.writeto(bluefile, blue_img.array, overwrite=True) fits.writeto(bluefile, blue_img.array, overwrite=True)
fits.writeto(redfile, red_img.array, overwrite=True) fits.writeto(redfile, red_img.array, overwrite=True)
#### save flux #### # ### save flux ####
Flux = flux_1 Flux = flux_1
hdu1 = fits.PrimaryHDU(Flux) hdu1 = fits.PrimaryHDU(Flux)
...@@ -5829,7 +5829,7 @@ class IFSsimulator(): ...@@ -5829,7 +5829,7 @@ class IFSsimulator():
# fits.writeto(bluefile, self.image_b, overwrite=True) # fits.writeto(bluefile, self.image_b, overwrite=True)
# fits.writeto(redfile, self.image_r, overwrite=True) # fits.writeto(redfile, self.image_r, overwrite=True)
################################################################################################# # #############################################################################
def simulate(self, sourcein, simnumber): def simulate(self, sourcein, simnumber):
""" """
...@@ -5890,7 +5890,7 @@ class IFSsimulator(): ...@@ -5890,7 +5890,7 @@ class IFSsimulator():
################################################################### ###################################################################
if self.appbianpai: if self.appbianpai:
# load yunxingbianpai csv file # load yunxingbianpai csv file
############ load star data catlog ##################### # ########### load star data catlog #####################
starcat = self.information['bianpai_file'] starcat = self.information['bianpai_file']
# starcat='selection_20230517_concat.fits' # starcat='selection_20230517_concat.fits'
...@@ -5902,7 +5902,7 @@ class IFSsimulator(): ...@@ -5902,7 +5902,7 @@ class IFSsimulator():
df = pd.read_csv( df = pd.read_csv(
self.information['dir_path']+'IFS_inputdata/TianCe/'+starcat) self.information['dir_path']+'IFS_inputdata/TianCe/'+starcat)
################################################################### ##############################################################
sn = self.simnumber-1 sn = self.simnumber-1
arr = np.array(df['time'][sn*5:sn*5+5]) arr = np.array(df['time'][sn*5:sn*5+5])
self.exptime_start_jd, self.exptime_start_index = find_min(arr) self.exptime_start_jd, self.exptime_start_index = find_min(arr)
...@@ -5911,10 +5911,10 @@ class IFSsimulator(): ...@@ -5911,10 +5911,10 @@ class IFSsimulator():
# self.earthshine_theta=df['earth_angle'][sn*5+index] # in degree # self.earthshine_theta=df['earth_angle'][sn*5+index] # in degree
self.dt = jd2time(self.exptime_start_jd) self.dt = jd2time(self.exptime_start_jd)
self.bianpai_data = df self.bianpai_data = df
################################################################## # ############################################################
else: else:
#### load orbit parameters ##### # ### load orbit parameters #####
flag = 0 flag = 0
self.dt = datetime.utcnow() self.dt = datetime.utcnow()
self.dt_num = int( self.dt_num = int(
...@@ -5937,7 +5937,7 @@ class IFSsimulator(): ...@@ -5937,7 +5937,7 @@ class IFSsimulator():
if flag == 1: if flag == 1:
break break
################################################################## # ############################################################
if kk + self.dt_num < len(d[:, 0]): if kk + self.dt_num < len(d[:, 0]):
...@@ -6040,7 +6040,7 @@ class IFSsimulator(): ...@@ -6040,7 +6040,7 @@ class IFSsimulator():
else: else:
print('Souce is not correct and programe will return') print('Souce is not correct and programe will return')
raise ValueError('Souce is not correct and programe will return') raise ValueError('Souce is not correct and programe will return')
############################################################################### # ########################################################################
# save the original image firstly; # save the original image firstly;
self.image_b_ori = self.image_b+10.0 self.image_b_ori = self.image_b+10.0
self.image_r_ori = self.image_r+10.0 self.image_r_ori = self.image_r+10.0
...@@ -6066,7 +6066,7 @@ class IFSsimulator(): ...@@ -6066,7 +6066,7 @@ class IFSsimulator():
for idk in range(exptimes): for idk in range(exptimes):
############ add some effect to images ####################################### # ### add some effect to images #################################
self.image_b = self.image_b_ori-10.0 self.image_b = self.image_b_ori-10.0
self.image_r = self.image_r_ori-10.0 self.image_r = self.image_r_ori-10.0
...@@ -6192,7 +6192,7 @@ def runIFSsim(sourcein, configfile, dir_path, result_path, iLoop, debug, applyho ...@@ -6192,7 +6192,7 @@ def runIFSsim(sourcein, configfile, dir_path, result_path, iLoop, debug, applyho
else: else:
simulate[iLoop].information['holemask'] = 'no' simulate[iLoop].information['holemask'] = 'no'
### dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/') # ## dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
# simulate[iLoop].information['dir_path'] = dir_path # simulate[iLoop].information['dir_path'] = dir_path
# simulate[iLoop].information['debug'] = debug # simulate[iLoop].information['debug'] = debug
...@@ -6203,4 +6203,4 @@ def runIFSsim(sourcein, configfile, dir_path, result_path, iLoop, debug, applyho ...@@ -6203,4 +6203,4 @@ def runIFSsim(sourcein, configfile, dir_path, result_path, iLoop, debug, applyho
return 1 return 1
############################## end ########################################## # ######################### end ##########################################
tests/test_ifs_sim.py
View file @ 5c4a3a79
...@@ -38,7 +38,7 @@ class TestDemoFunction(unittest.TestCase): ...@@ -38,7 +38,7 @@ class TestDemoFunction(unittest.TestCase):
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' # ##configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
...@@ -77,7 +77,7 @@ class TestDemoFunction(unittest.TestCase): ...@@ -77,7 +77,7 @@ class TestDemoFunction(unittest.TestCase):
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' # ##configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
...@@ -115,7 +115,7 @@ class TestDemoFunction(unittest.TestCase): ...@@ -115,7 +115,7 @@ class TestDemoFunction(unittest.TestCase):
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' # ##configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
...@@ -153,7 +153,7 @@ class TestDemoFunction(unittest.TestCase): ...@@ -153,7 +153,7 @@ class TestDemoFunction(unittest.TestCase):
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' # ##configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
...@@ -192,7 +192,7 @@ class TestDemoFunction(unittest.TestCase): ...@@ -192,7 +192,7 @@ class TestDemoFunction(unittest.TestCase):
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' # ##configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
...@@ -229,7 +229,7 @@ class TestDemoFunction(unittest.TestCase): ...@@ -229,7 +229,7 @@ class TestDemoFunction(unittest.TestCase):
dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/') dir_path = os.path.join(os.environ['UNIT_TEST_DATA_ROOT'], 'ifs_sim/')
print(dir_path) print(dir_path)
print(sys.version) print(sys.version)
###configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config' # ##configfile = dir_path+'IFS_inputdata/configData/IFS_sim_C90.config'
configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config' configfile = './csst_ifs_sim/ifs_data/IFS_sim_C90.config'
......
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