test (22e2bedd) · Commits · csst-sims / csst_mci_sim

csst_mci_sim/csst_mci_sim.py

+78 −77

Original line number	Diff line number	Diff line
		@@ -1738,8 +1738,6 @@ class MCIsimulator():
		# np.save('umag.npy', df3['umag'] )
		self.star=df3



		del da,df2,df3

		self.information['ra_obj'] = self.star['ra_gaia'].mean()
		@@ -1760,109 +1758,109 @@ class MCIsimulator():
		#################### cal_PSF_array #################################


		if self.save_starpsf==100:
		# if self.save_starpsf==100:

		self.log.info('calculate and save star PSF data......')
		# self.log.info('calculate and save star PSF data......')

		primary_g=fits.PrimaryHDU()
		PSF_g=fits.HDUList([primary_g])
		# primary_g=fits.PrimaryHDU()
		# PSF_g=fits.HDUList([primary_g])

		primary_r=fits.PrimaryHDU()
		PSF_r=fits.HDUList([primary_r])
		# primary_r=fits.PrimaryHDU()
		# PSF_r=fits.HDUList([primary_r])

		primary_i=fits.PrimaryHDU()
		PSF_i=fits.HDUList([primary_i])
		# primary_i=fits.PrimaryHDU()
		# PSF_i=fits.HDUList([primary_i])

		fov=0.05*min(self.information['xsize'],self.information['ysize'])
		# fov=0.05*min(self.information['xsize'],self.information['ysize'])

		dec_arr,ra_arr=make_c_coor(fov,25)
		# dec_arr,ra_arr=make_c_coor(fov,25)

		k=0;
		# k=0;

		#######################################################################
		for ii in range(len(ra_arr)):
		for jj in range(len(ra_arr)):
		# #######################################################################
		# for ii in range(len(ra_arr)):
		# for jj in range(len(ra_arr)):

		##################################################################
		# ##################################################################

		galRa = ra_arr[ii,jj] +center_ra*3600 # ra of PFS, arcsecond
		galDec = dec_arr[ii,jj]+center_dec*3600 # dec of PSF, arcsecond
		# galRa = ra_arr[ii,jj] +center_ra*3600 # ra of PFS, arcsecond
		# galDec = dec_arr[ii,jj]+center_dec*3600 # dec of PSF, arcsecond

		fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, galRa, galDec)
		# fsx,fsy=cal_pos(center_ra,center_dec, rotTelPos, rotSkyPos, galRa, galDec)

		############# do field distottion ##########
		# ############# do field distottion ##########

		if self.distortion:
		# if self.distortion:

		for i in range(3):
		ch=channel[i]
		fpx,fpy=distortField(fsx, fsy, ch)
		# for i in range(3):
		# ch=channel[i]
		# fpx,fpy=distortField(fsx, fsy, ch)

		else:
		# else:

		fpx=fsx
		fpy=fsy
		# fpx=fsx
		# fpy=fsy

		###################################################################
		# ###################################################################

		psf=dict()
		# psf=dict()

		for i in range(3):
		# for i in range(3):

		ch=channel[i]
		# ch=channel[i]

		psfmat=self.get_PSF(fpx, fpy, ch)
		# psfmat=self.get_PSF(fpx, fpy, ch)

		temp=0
		# temp=0

		for iwave in range(7):
		temp=temp+1/7.0*psfmat[:,:,iwave]
		# for iwave in range(7):
		# temp=temp+1/7.0*psfmat[:,:,iwave]

		temp=temp/temp.sum()
		# temp=temp/temp.sum()

		####rotate the PSF data
		if abs(theta.deg)>0:
		psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=False) # here we choose reshape=False, the rotated image will
		else:
		psf[ch]=temp
		###################################
		# ####rotate the PSF data
		# if abs(theta.deg)>0:
		# psf[ch]=ndimage.rotate(temp, theta.deg, order=1, reshape=False) # here we choose reshape=False, the rotated image will
		# else:
		# psf[ch]=temp
		# ###################################

		hdu_g=fits.ImageHDU(psf['g'])
		hdu_g.header['ra'] = ra_arr[ii,jj]
		hdu_g.header['dec'] = dec_arr[ii,jj]
		hdu_g.header['rot_deg']=-theta.deg
		PSF_g.append(hdu_g)
		###################################
		hdu_r=fits.ImageHDU(psf['r'])
		hdu_r.header['ra'] = ra_arr[ii,jj]
		hdu_r.header['dec'] = dec_arr[ii,jj]
		hdu_r.header['rot_deg']=-theta.deg
		PSF_r.append(hdu_r)
		###################################
		hdu_i=fits.ImageHDU(psf['i'])
		hdu_i.header['ra'] = ra_arr[ii,jj]
		hdu_i.header['dec']= dec_arr[ii,jj]
		hdu_i.header['rot_deg']=-theta.deg
		PSF_i.append(hdu_i)
		###################################
		del hdu_g
		del hdu_r
		del hdu_i
		k=k+1
		# hdu_g=fits.ImageHDU(psf['g'])
		# hdu_g.header['ra'] = ra_arr[ii,jj]
		# hdu_g.header['dec'] = dec_arr[ii,jj]
		# hdu_g.header['rot_deg']=-theta.deg
		# PSF_g.append(hdu_g)
		# ###################################
		# hdu_r=fits.ImageHDU(psf['r'])
		# hdu_r.header['ra'] = ra_arr[ii,jj]
		# hdu_r.header['dec'] = dec_arr[ii,jj]
		# hdu_r.header['rot_deg']=-theta.deg
		# PSF_r.append(hdu_r)
		# ###################################
		# hdu_i=fits.ImageHDU(psf['i'])
		# hdu_i.header['ra'] = ra_arr[ii,jj]
		# hdu_i.header['dec']= dec_arr[ii,jj]
		# hdu_i.header['rot_deg']=-theta.deg
		# PSF_i.append(hdu_i)
		# ###################################
		# del hdu_g
		# del hdu_r
		# del hdu_i
		# k=k+1

		############################################
		file_g=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C1.fits'
		PSF_g.writeto(file_g,overwrite=True)
		# ############################################
		# file_g=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C1.fits'
		# PSF_g.writeto(file_g,overwrite=True)

		file_r=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C2.fits'
		PSF_r.writeto(file_r,overwrite=True)
		# file_r=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C2.fits'
		# PSF_r.writeto(file_r,overwrite=True)

		file_i=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C3.fits'
		PSF_i.writeto(file_i,overwrite=True)
		# file_i=self.result_path+'/PSF_Data/'+'PSF_sim_No.'+str(self.information['simnumber'])+'_C3.fits'
		# PSF_i.writeto(file_i,overwrite=True)

		del PSF_g
		del PSF_r
		del PSF_i
		# del PSF_g
		# del PSF_r
		# del PSF_i

		########## finish save PSF fits ##########################################
		############################################################################
		@@ -1878,6 +1876,7 @@ class MCIsimulator():


		nsrcs=len(self.star['ra_gaia']) + 8

		self.log.info('load star catlog successfully')

		##################################################################
		@@ -1910,6 +1909,8 @@ class MCIsimulator():
		#######################################################################
		nlayccd = 0


		self.log.info('load star catlog successfully 0000000')
		############################################

		ra = self.information['ra_pnt0']
		@@ -1922,11 +1923,11 @@ class MCIsimulator():
		z_sat=float(self.orbit_pars[self.orbit_exp_num,3])


		self.log.info('load star catlog 111')
		self.log.info('load star catlog 111111111111111')

		wave0, zodi0 = self.zodiacal(ra, dec, self.TianCe_day) # erg/s/cm^2/A/arcsec^2

		self.log.info('load star catlog 2222')
		self.log.info('load star catlog 2222222222222222222')
		# EarthShine from straylight
		sl = StrayLight(self.information['dir_path'],jtime=time_jd, sat=np.array([x_sat, y_sat, z_sat]),
		radec=np.array([(rau.degree).value, (decu.degree).value]))