more case test (b4201fbc) · Commits · csst-sims / csst_ifs_sim

csst_ifs_sim/csst_ifs_sim.py

+230 −230

Original line number	Diff line number	Diff line
		@@ -117,50 +117,50 @@ def transRaDec2D(ra, dec):
		###############################################################################


		def flux2ill(wave, flux):
		"""
		# def flux2ill(wave, flux):
		# """


		Parameters
		----------
		wave : TYPE
		DESCRIPTION.
		flux : TYPE
		DESCRIPTION.
		# Parameters
		# ----------
		# wave : TYPE
		# DESCRIPTION.
		# flux : TYPE
		# DESCRIPTION.

		Returns
		-------
		E : TYPE
		DESCRIPTION.
		# Returns
		# -------
		# E : TYPE
		# DESCRIPTION.

		"""
		# """

		# erg/s/cm^2/A/arcsec^2 to W/m^2
		# # erg/s/cm^2/A/arcsec^2 to W/m^2

		# 1 W/m^2/sr/μm = 0.10 erg/cm^2/s/sr/A
		# 1 sr = 1 rad^2 = 4.25452e10 arcsec^2
		# 1 J/s = 1 W
		# 1 J = 10^7 erg
		# # 1 W/m^2/sr/μm = 0.10 erg/cm^2/s/sr/A
		# # 1 sr = 1 rad^2 = 4.25452e10 arcsec^2
		# # 1 J/s = 1 W
		# # 1 J = 10^7 erg

		# convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr
		flux1 = flux / (1/4.25452e10)
		# convert erg/s/cm^2/A/sr to W/m^2/sr/um
		flux2 = flux1 * 10
		# # convert erg/s/cm^2/A/arcsec^2 to erg/s/cm^2/A/sr
		# flux1 = flux / (1/4.25452e10)
		# # convert erg/s/cm^2/A/sr to W/m^2/sr/um
		# flux2 = flux1 * 10

		# 对接收面积积分，输出单位 W/m^2/nm
		D = 2 # meter
		f = 28 # meter
		flux3 = flux2 * np.pi * D2 / 4 / f2 / 10**3
		# # 对接收面积积分，输出单位 W/m^2/nm
		# D = 2 # meter
		# f = 28 # meter
		# flux3 = flux2 * np.pi * D2 / 4 / f2 / 10**3

		# 对波长积分
		f = interp1d(wave, flux3)
		wave_interp = np.arange(3800, 7800)
		flux3_interp = f(wave_interp)
		# 输出单位 W/m^2
		delta_lamba = 0.1 # nm
		E = np.sum(flux3_interp * delta_lamba)
		# # 对波长积分
		# f = interp1d(wave, flux3)
		# wave_interp = np.arange(3800, 7800)
		# flux3_interp = f(wave_interp)
		# # 输出单位 W/m^2
		# delta_lamba = 0.1 # nm
		# E = np.sum(flux3_interp * delta_lamba)

		return E
		# return E

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

		@@ -391,25 +391,25 @@ class StrayLight(object):
		###############################################################################


		def str2time(strTime):
		"""
		# def str2time(strTime):
		# """


		Parameters
		----------
		strTime : TYPE
		DESCRIPTION.
		# Parameters
		# ----------
		# strTime : TYPE
		# DESCRIPTION.

		Returns
		-------
		TYPE
		DESCRIPTION.
		# Returns
		# -------
		# TYPE
		# DESCRIPTION.

		"""
		if len(strTime) > 20: # 暂时未用到
		msec = int(float('0.'+strTime[20:])*1000000) # 微秒
		str2 = strTime[0:19]+' '+str(msec)
		return datetime.strptime(str2, '%Y %m %d %H %M %S %f')
		# """
		# if len(strTime) > 20: # 暂时未用到
		# msec = int(float('0.'+strTime[20:])*1000000) # 微秒
		# str2 = strTime[0:19]+' '+str(msec)
		# return datetime.strptime(str2, '%Y %m %d %H %M %S %f')

		# datetime类转mjd
		##########################################################################
		@@ -822,40 +822,40 @@ def centroid(data):
		return float(cx), float(cy)

		###############################################################################
		def centroidN(data):
		"""
		# def centroidN(data):
		# """


		Parameters
		----------
		data : TYPE
		DESCRIPTION.
		# Parameters
		# ----------
		# data : TYPE
		# DESCRIPTION.

		Returns
		-------
		cx : TYPE
		DESCRIPTION.
		cy : TYPE
		DESCRIPTION.
		# Returns
		# -------
		# cx : TYPE
		# DESCRIPTION.
		# cy : TYPE
		# DESCRIPTION.

		"""
		'''
		calculate the centroid of the input two-dimentional image data
		# """
		# '''
		# calculate the centroid of the input two-dimentional image data

		Parameters
		----------
		data : input image.
		# Parameters
		# ----------
		# data : input image.

		Returns
		-------
		cx: the centroid column number, in horizontal direction definet in python image show
		cy: the centroid row number , in vertical direction
		# Returns
		# -------
		# cx: the centroid column number, in horizontal direction definet in python image show
		# cy: the centroid row number , in vertical direction

		'''
		###
		from scipy import ndimage
		cy, cx = ndimage.center_of_mass(data)
		return cx, cy
		# '''
		# ###
		# from scipy import ndimage
		# cy, cx = ndimage.center_of_mass(data)
		# return cx, cy

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

		@@ -1494,42 +1494,42 @@ class IFSsimulator():
		return wave_A, spec_erg2
		##########################################################################

		def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)):
		"""
		# def smoothingWithChargeDiffusion(self, image, sigma=(0.32, 0.32)):
		# """


		Parameters
		----------
		image : TYPE
		DESCRIPTION.
		sigma : TYPE, optional
		DESCRIPTION. The default is (0.32, 0.32).
		# Parameters
		# ----------
		# image : TYPE
		# DESCRIPTION.
		# sigma : TYPE, optional
		# DESCRIPTION. The default is (0.32, 0.32).

		Returns
		-------
		TYPE
		DESCRIPTION.
		# Returns
		# -------
		# TYPE
		# DESCRIPTION.

		"""
		"""
		Smooths a given image with a gaussian kernel with widths given as sigmas.
		This smoothing can be used to mimic charge diffusion within the CCD.
		# """
		# """
		# Smooths a given image with a gaussian kernel with widths given as sigmas.
		# This smoothing can be used to mimic charge diffusion within the CCD.

		The default values are from Table 8-2 of CCD_273_Euclid_secification_1.0.130812.pdf converted
		to sigmas (FWHM / (2sqrt(2ln2)) and rounded up to the second decimal.
		# The default values are from Table 8-2 of CCD_273_Euclid_secification_1.0.130812.pdf converted
		# to sigmas (FWHM / (2sqrt(2ln2)) and rounded up to the second decimal.

		.. Note:: This method should not be called for the full image if the charge spreading
		has already been taken into account in the system PSF to avoid double counting.
		# .. Note:: This method should not be called for the full image if the charge spreading
		# has already been taken into account in the system PSF to avoid double counting.

		:param image: image array which is smoothed with the kernel
		:type image: ndarray
		:param sigma: widths of the gaussian kernel that approximates the charge diffusion [0.32, 0.32].
		:param sigma: tuple
		# :param image: image array which is smoothed with the kernel
		# :type image: ndarray
		# :param sigma: widths of the gaussian kernel that approximates the charge diffusion [0.32, 0.32].
		# :param sigma: tuple

		:return: smoothed image array
		:rtype: ndarray
		"""
		return ndimage.filters.gaussian_filter(image, sigma)
		# :return: smoothed image array
		# :rtype: ndarray
		# """
		# return ndimage.filters.gaussian_filter(image, sigma)
		###############################################################################

		def readCosmicRayInformation(self):
		@@ -1835,88 +1835,88 @@ class IFSsimulator():
		######################################################################

		##############################################################################
		def generateflat(self, ave=1.0, sigma=0.01):
		"""
		# def generateflat(self, ave=1.0, sigma=0.01):
		# """


		Parameters
		----------
		ave : TYPE, optional
		DESCRIPTION. The default is 1.0.
		sigma : TYPE, optional
		DESCRIPTION. The default is 0.01.
		# Parameters
		# ----------
		# ave : TYPE, optional
		# DESCRIPTION. The default is 1.0.
		# sigma : TYPE, optional
		# DESCRIPTION. The default is 0.01.

		Returns
		-------
		TYPE
		DESCRIPTION.
		TYPE
		DESCRIPTION.
		# Returns
		# -------
		# TYPE
		# DESCRIPTION.
		# TYPE
		# DESCRIPTION.

		"""
		"""
		Creates a flat field image with given properties.
		# """
		# """
		# Creates a flat field image with given properties.

		:return: flat field image
		:rtype: ndarray
		"""
		self.log.info('Generating a flat field...')
		self.log.info('The flat field has mean value of 1 and a given fluctuations, usually either 1 or 2 percent defined by sigma= %d...' % sigma)
		# :return: flat field image
		# :rtype: ndarray
		# """
		# self.log.info('Generating a flat field...')
		# self.log.info('The flat field has mean value of 1 and a given fluctuations, usually either 1 or 2 percent defined by sigma= %d...' % sigma)

		np.random.seed(5*self.simnumber)
		self.flat_b = np.random.normal(loc=ave, scale=sigma, size=(2048, 4096))
		# np.random.seed(5*self.simnumber)
		# self.flat_b = np.random.normal(loc=ave, scale=sigma, size=(2048, 4096))

		np.random.seed(55*self.simnumber)
		self.flat_r = np.random.normal(loc=ave, scale=sigma, size=(3072, 6144))
		# np.random.seed(55*self.simnumber)
		# self.flat_r = np.random.normal(loc=ave, scale=sigma, size=(3072, 6144))

		s1 = self.flat_b
		# s1 = self.flat_b

		hdu1 = fits.PrimaryHDU(s1)
		# hdu1 = fits.PrimaryHDU(s1)

		hdu1.header.set('sigma', sigma)
		# hdu1.header.set('sigma', sigma)

		dtime = datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S')
		hdu1.header.add_history(
		'flat image of blue channel is generated on :'+dtime)
		# dtime = datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S')
		# hdu1.header.add_history(
		# 'flat image of blue channel is generated on :'+dtime)

		f1 = '../flat_Blue_'+str(sigma)+'.fits'
		fits.writeto(f1, s1, header=hdu1.header, overwrite=True)
		# f1 = '../flat_Blue_'+str(sigma)+'.fits'
		# fits.writeto(f1, s1, header=hdu1.header, overwrite=True)

		s2 = self.flat_r
		# s2 = self.flat_r

		hdu1 = fits.PrimaryHDU(s2)
		# hdu1 = fits.PrimaryHDU(s2)

		hdu1.header.set('sigma', sigma)
		# hdu1.header.set('sigma', sigma)

		dtime = datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S')
		hdu1.header.add_history(
		'flat image of red channel is generated on :'+dtime)
		# dtime = datetime.utcnow().strftime('%Y -%m -%d %H: %M: %S')
		# hdu1.header.add_history(
		# 'flat image of red channel is generated on :'+dtime)

		f2 = '../flat_Red_'+str(sigma)+'.fits'
		fits.writeto(f2, s2, header=hdu1.header, overwrite=True)
		# f2 = '../flat_Red_'+str(sigma)+'.fits'
		# fits.writeto(f2, s2, header=hdu1.header, overwrite=True)

		return self.flat_b, self.flat_r
		# return self.flat_b, self.flat_r

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

		def addLampFlux(self):
		"""
		# def addLampFlux(self):
		# """


		Returns
		-------
		None.
		# Returns
		# -------
		# None.

		"""
		"""
		Include flux from the calibration source.
		"""
		# """
		# """
		# Include flux from the calibration source.
		# """

		self.image_b += fits.getdata(self.information['flatflux'])
		self.image_r += fits.getdata(self.information['flatflux'])
		# self.image_b += fits.getdata(self.information['flatflux'])
		# self.image_r += fits.getdata(self.information['flatflux'])

		self.log.info('Flux from the calibration unit included (%s)' %
		self.information['flatflux'])
		# self.log.info('Flux from the calibration unit included (%s)' %
		# self.information['flatflux'])

		#############################################################################
		def MakeFlatMatrix(self, img, seed):
		@@ -2111,69 +2111,69 @@ class IFSsimulator():
		##########################################################
		#########################################################################

		def addReadoutTrails(self):
		"""


		Returns
		-------
		None.

		"""
		"""
		Add readout trails resulting from reading out the shutter open.

		Quadrants assumed to be numbered:
		2 3
		0 1
		"""
		flux_ratio = self.information['readouttime'] / float(
		self.information['bluesize']) / self.information['exptime']

		# make a copy, this will be updated
		data = self.image_b.copy()
		# Amplifier at different positions depending on the quadrant number !
		# left side is 0, 2 and right side is 1, 3 starting from bottom i.e.
		# going clock wise from lower left we have 0, 2, 3, and 1 quadrants.
		if self.information['quadrant'] in (2, 3):
		data = data[::-1, :]

		data_shift = data.copy() * flux_ratio
		size1, size2 = data.shape

		for i in range(1, size2, 1):
		data_shift2 = np.roll(data_shift, i, axis=0)
		data_shift2[:i, :] = 0.0
		data += data_shift2

		if self.information['quadrant'] in (2, 3):
		self.image_b = data[::-1, :]
		else:
		self.image_b = data

		# def addReadoutTrails(self):
		# """

		flux_ratio = self.information['readouttime'] / float(
		self.information['redsize']) / self.information['exptime']
		# make a copy, this will be updated
		data = self.image_r.copy()
		# Amplifier at different positions depending on the quadrant number !
		# left side is 0, 2 and right side is 1, 3 starting from bottom i.e.
		# going clock wise from lower left we have 0, 2, 3, and 1 quadrants.
		if self.information['quadrant'] in (2, 3):
		data = data[::-1, :]

		data_shift = data.copy() * flux_ratio
		size1, size2 = data.shape
		# Returns
		# -------
		# None.

		for i in range(1, size2, 1):
		data_shift2 = np.roll(data_shift, i, axis=0)
		data_shift2[:i, :] = 0.0
		data += data_shift2
		# """
		# """
		# Add readout trails resulting from reading out the shutter open.

		if self.information['quadrant'] in (2, 3):
		self.image_r = data[::-1, :]
		else:
		self.image_r = data
		# Quadrants assumed to be numbered:
		# 2 3
		# 0 1
		# """
		# flux_ratio = self.information['readouttime'] / float(
		# self.information['bluesize']) / self.information['exptime']

		# # make a copy, this will be updated
		# data = self.image_b.copy()
		# # Amplifier at different positions depending on the quadrant number !
		# # left side is 0, 2 and right side is 1, 3 starting from bottom i.e.
		# # going clock wise from lower left we have 0, 2, 3, and 1 quadrants.
		# if self.information['quadrant'] in (2, 3):
		# data = data[::-1, :]

		# data_shift = data.copy() * flux_ratio
		# size1, size2 = data.shape

		# for i in range(1, size2, 1):
		# data_shift2 = np.roll(data_shift, i, axis=0)
		# data_shift2[:i, :] = 0.0
		# data += data_shift2

		# if self.information['quadrant'] in (2, 3):
		# self.image_b = data[::-1, :]
		# else:
		# self.image_b = data


		# flux_ratio = self.information['readouttime'] / float(
		# self.information['redsize']) / self.information['exptime']
		# # make a copy, this will be updated
		# data = self.image_r.copy()
		# # Amplifier at different positions depending on the quadrant number !
		# # left side is 0, 2 and right side is 1, 3 starting from bottom i.e.
		# # going clock wise from lower left we have 0, 2, 3, and 1 quadrants.
		# if self.information['quadrant'] in (2, 3):
		# data = data[::-1, :]

		# data_shift = data.copy() * flux_ratio
		# size1, size2 = data.shape

		# for i in range(1, size2, 1):
		# data_shift2 = np.roll(data_shift, i, axis=0)
		# data_shift2[:i, :] = 0.0
		# data += data_shift2

		# if self.information['quadrant'] in (2, 3):
		# self.image_r = data[::-1, :]
		# else:
		# self.image_r = data

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

		@@ -2259,24 +2259,24 @@ class IFSsimulator():

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

		def applyScatteredLight(self):
		"""
		# def applyScatteredLight(self):
		# """


		Returns
		-------
		None.
		# Returns
		# -------
		# None.

		"""
		"""
		Adds spatially uniform scattered light to the image.
		"""
		sl = self.information['exptime'] * self.information['scattered_light']
		# """
		# """
		# Adds spatially uniform scattered light to the image.
		# """
		# sl = self.information['exptime'] * self.information['scattered_light']

		self.image_b += sl
		self.image_r += sl
		# self.image_b += sl
		# self.image_r += sl

		self.log.info('Added scattered light = %f' % sl)
		# self.log.info('Added scattered light = %f' % sl)

		##############################################################################
		def applyPoissonNoise(self):