Merge remote-tracking branch 'origin/new_sim_sls' into new_sim (716634ad) · Commits · Liu Dezi / csst_msc_sim

ObservationSim/Instrument/Chip/Chip.py

+4 −2

Original line number	Diff line number	Diff line
		@@ -270,7 +270,7 @@ class Chip(FocalPlane):
		noise = self.dark_noise * exptime + self.read_noise**2
		return noise

		def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='SCI', sky_map=None, tel=None, logger=None):
		def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='SCI', sky_map=None, post_flash_map=None, tel=None, logger=None):
		# Set random seeds
		SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"])
		SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"])
		@@ -309,6 +309,9 @@ class Chip(FocalPlane):
		if config["output_setting"]["flat_output"] == False:
		del flat_img

		if post_flash_map is not None:
		img = img + post_flash_map

		# Apply Shutter-effect for one chip
		if config["ins_effects"]["shutter_effect"] == True:
		chip_utils.log_info(msg=" Apply shutter effect", logger=self.logger)
		@@ -320,7 +323,6 @@ class Chip(FocalPlane):
		shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" % (chip_output.subdir, self.chipID))
		del shutt_gsimg
		del shuttimg

		# # Add Poisson noise to the resulting images
		# # (NOTE): this can only applied to the slitless image
		# # since it dose not use photon shooting to draw stamps

ObservationSim/MockObject/FlatLED.py

0 → 100755

+350 −0

Original line number	Diff line number	Diff line


		import galsim
		import os, sys
		import numpy as np
		from astropy.io import fits
		from scipy.interpolate import griddata
		import math
		import astropy.constants as cons
		from astropy.table import Table
		from ObservationSim.MockObject.SpecDisperser import SpecDisperser
		import time
		from scipy import interpolate

		from ObservationSim.MockObject.MockObject import MockObject
		# from ObservationSim.Straylight import calculateSkyMap_split_g

		# flatDir = '/Volumes/EAGET/LED_FLAT/'
		LED_name = ['LED1', 'LED2', 'LED3', 'LED4', 'LED5', 'LED6', 'LED7', 'LED8', 'LED9', 'LED10', 'LED11', 'LED12', 'LED13',
		'LED14']
		cwaves_name = {'LED1': '275', 'LED2': '310', 'LED3': '430', 'LED4': '505', 'LED5': '545', 'LED6': '590', 'LED7': '670',
		'LED8': '760', 'LED9': '940', 'LED10': '940', 'LED11': '1050', 'LED12': '1050',
		'LED13': '340', 'LED14': '365'}

		cwaves = {'LED1': 2750, 'LED2': 3100, 'LED3': 4300, 'LED4': 5050, 'LED5': 5250, 'LED6': 5900, 'LED7': 6700,
		'LED8': 7600, 'LED9': 8800, 'LED10': 9400, 'LED11': 10500, 'LED12': 15500, 'LED13': 3400, 'LED14': 3650}
		cwaves_fwhm = {'LED1': 110, 'LED2': 120, 'LED3': 200, 'LED4': 300, 'LED5': 300, 'LED6': 130, 'LED7': 210,
		'LED8': 260, 'LED9': 400, 'LED10': 370, 'LED11': 500, 'LED12': 1400, 'LED13': 90, 'LED14': 100}
		# LED_QE = {'LED1': 0.3, 'LED2': 0.4, 'LED13': 0.5, 'LED14': 0.5, 'LED10': 0.4}
		# e-/ms
		fluxLED = {'LED1': 0.16478729, 'LED2': 0.084220931, 'LED3': 2.263360617, 'LED4': 2.190623489, 'LED5': 0.703504768,
		'LED6': 0.446117963, 'LED7': 0.647122098, 'LED8': 0.922313442,
		'LED9': 0.987278143, 'LED10': 2.043989167, 'LED11': 0.612571429, 'LED12': 1.228915663, 'LED13': 0.17029384,
		'LED14': 0.27842925}

		mirro_eff = {'GU':0.61, 'GV':0.8, 'GI':0.8}

		class FlatLED(object):
		def __init__(self, chip,filt, flatDir = '/Users/zhangxin/Work/SlitlessSim/csst_sls_calibration/flat_field_cube/models/', logger=None):
		# self.led_type_list = led_type_list
		self.flatDir = flatDir
		self.filt = filt
		self.chip = chip
		self.logger = logger

		###
		### return LED flat, e/s
		###
		def getLEDImage(self, led_type='LED1'):
		# cwave = cwaves[led_type]
		flat = fits.open(self.flatDir + 'model_' + cwaves_name[led_type] + 'nm.fits')
		xlen = flat[0].header['NAXIS1']
		ylen = 601
		x = np.linspace(0, self.chip.npix_x * 6, xlen)
		y = np.linspace(0, self.chip.npix_y * 5, ylen)
		xx, yy = np.meshgrid(x, y)

		a1 = flat[0].data[0:ylen, 0:xlen]
		# z = np.sin((xx+yy+xx2+yy2))
		# fInterp = interp2d(xx, yy, z, kind='linear')

		X_ = np.hstack((xx.flatten()[:, None], yy.flatten()[:, None]))
		Z_ = a1.flatten()

		n_x = np.arange(0, self.chip.npix_x * 6, 1)
		n_y = np.arange(0, self.chip.npix_y * 5, 1)

		M, N = np.meshgrid(n_x, n_y)

		i = self.chip.rowID - 1
		j = self.chip.colID - 1
		U = griddata(X_, Z_, (
		M[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)],
		N[self.chip.npix_y * i:self.chip.npix_y * (i + 1), self.chip.npix_x * j:self.chip.npix_x * (j + 1)]),
		method='cubic')
		U = U/np.mean(U)
		flatImage = UfluxLED[led_type]1000
		return flatImage

		def drawObj_LEDFlat_img(self, led_type_list=['LED1'], exp_t_list=[0.1]):
		if len(led_type_list) > len(exp_t_list):
		return np.ones([self.chip.npix_y,self.chip.npix_x])

		ledFlat = np.zeros([self.chip.npix_y,self.chip.npix_x])
		for i in np.arange(len(led_type_list)):
		led_type = led_type_list[i]
		exp_t = exp_t_list[i]
		unitFlatImg = self.getLEDImage(led_type=led_type)
		led_wave = cwaves[led_type]
		led_fwhm = cwaves_fwhm[led_type]
		led_spec = self.gaussian1d_profile_led(led_wave, led_fwhm)
		speci = interpolate.interp1d(led_spec['WAVELENGTH'], led_spec['FLUX'])
		w_list = np.arange(self.filt.blue_limit, self.filt.red_limit, 0.5) #A

		f_spec = speci(w_list)
		ccd_bp = self.chip._getChipEffCurve(self.chip.filter_type)
		ccd_eff = ccd_bp.__call__(w_list / 10.)
		filt_bp = self.filt.filter_bandpass
		fil_eff = filt_bp.__call__(w_list / 10.)
		t_spec = np.trapz(f_specccd_efffil_eff, w_list)
		# print(i, np.mean(unitFlatImg), t_spec, exp_t)
		unitFlatImg = unitFlatImg * t_spec

		ledFlat = ledFlat+unitFlatImg*exp_t
		return ledFlat

		def drawObj_LEDFlat_slitless(self, led_type_list=['LED1'], exp_t_list=[0.1]):
		if len(led_type_list) != len(exp_t_list):
		return np.ones([self.chip.npix_y,self.chip.npix_x])

		ledFlat = np.zeros([self.chip.npix_y,self.chip.npix_x])
		for i in np.arange(len(led_type_list)):
		led_type = led_type_list[i]
		exp_t = exp_t_list[i]
		unitFlatImg = self.getLEDImage(led_type=led_type)
		ledFlat_ = unitFlatImg*exp_t
		ledFlat_ = ledFlat_ / mirro_eff[self.filt.filter_type]
		ledFlat_.astype(np.float32)
		led_wave = cwaves[led_type]
		led_fwhm = cwaves_fwhm[led_type]
		led_spec = self.gaussian1d_profile_led(led_wave, led_fwhm)
		ledspec_map = self.calculateLEDSpec(
		skyMap=ledFlat_,
		blueLimit=self.filt.blue_limit,
		redLimit=self.filt.red_limit,
		conf=self.chip.sls_conf,
		pixelSize=self.chip.pix_scale,
		isAlongY=0,
		flat_cube=self.chip.flat_cube, led_spec=led_spec)

		ledFlat = ledFlat + ledspec_map

		return ledFlat

		def drawObj_LEDFlat(self, led_type_list=['LED1'], exp_t_list=[0.1]):
		if self.chip.survey_type == "photometric":
		return self.drawObj_LEDFlat_img(led_type_list=led_type_list, exp_t_list=exp_t_list)
		elif self.chip.survey_type == "spectroscopic":
		return self.drawObj_LEDFlat_slitless(led_type_list=led_type_list, exp_t_list=exp_t_list)


		def gaussian1d_profile_led(self, xc=5050, fwhm=300):
		sigma = fwhm/2.355
		x_radii = int(5*sigma + 1)
		xlist = np.arange(xc-x_radii, xc+x_radii, 0.5)
		xlist_ = np.zeros(len(xlist) + 2)
		xlist_[1:-1] = xlist
		xlist_[0] = 2550
		xlist_[-1] = 10000
		data = np.exp((-(xlist-xc)(xlist-xc))/(2sigmasigma))/(np.sqrt(2math.pi)*sigma)
		data_ = np.zeros(len(xlist) + 2)
		data_[1:-1] = data
		return Table(np.array([xlist_.astype(np.float32), data_.astype(np.float32)]).T, names=('WAVELENGTH', 'FLUX'))

		def calculateLEDSpec(self, skyMap=None, blueLimit=4200, redLimit=6500,
		conf=[''], pixelSize=0.074, isAlongY=0,
		split_pos=3685, flat_cube=None, led_spec=None):

		conf1 = conf[0]
		conf2 = conf[0]
		if np.size(conf) == 2:
		conf2 = conf[1]

		skyImg = galsim.Image(skyMap, xmin=0, ymin=0)

		tbstart = blueLimit
		tbend = redLimit

		fimg = np.zeros_like(skyMap)

		fImg = galsim.Image(fimg)

		spec = led_spec
		if isAlongY == 0:
		directParm = 0
		if isAlongY == 1:
		directParm = 1

		if split_pos >= skyImg.array.shape[directParm]:
		skyImg1 = galsim.Image(skyImg.array)
		origin1 = [0, 0]
		# sdp = specDisperser.specDisperser(orig_img=skyImg1, xcenter=skyImg1.center.x, ycenter=skyImg1.center.y,
		# full_img=fimg, tar_spec=spec, band_start=tbstart, band_end=tbend,
		# origin=origin1,
		# conf=conf1)
		# sdp.compute_spec_orders()

		y_len = skyMap.shape[0]
		x_len = skyMap.shape[1]
		delt_x = 100
		delt_y = 100

		sub_y_start_arr = np.arange(0, y_len, delt_y)
		sub_y_end_arr = sub_y_start_arr + delt_y
		sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len)

		sub_x_start_arr = np.arange(0, x_len, delt_x)
		sub_x_end_arr = sub_x_start_arr + delt_x
		sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len)

		for i, k1 in enumerate(sub_y_start_arr):
		sub_y_s = k1
		sub_y_e = sub_y_end_arr[i]

		sub_y_center = (sub_y_s + sub_y_e) / 2.

		for j, k2 in enumerate(sub_x_start_arr):
		sub_x_s = k2
		sub_x_e = sub_x_end_arr[j]

		skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
		origin_sub = [sub_y_s, sub_x_s]
		sub_x_center = (sub_x_s + sub_x_e) / 2.

		sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center,
		origin=origin_sub,
		tar_spec=spec,
		band_start=tbstart, band_end=tbend,
		conf=conf2,
		flat_cube=flat_cube, ignoreBeam=['D', 'E'])

		spec_orders = sdp.compute_spec_orders()

		for k, v in spec_orders.items():
		img_s = v[0]
		origin_order_x = v[1]
		origin_order_y = v[2]
		ssImg = galsim.ImageF(img_s)
		ssImg.setOrigin(origin_order_x, origin_order_y)
		bounds = ssImg.bounds & fImg.bounds
		if bounds.area() == 0:
		continue
		fImg[bounds] = fImg[bounds] + ssImg[bounds]



		else:


		# sdp.compute_spec_orders()
		y_len = skyMap.shape[0]
		x_len = skyMap.shape[1]
		delt_x = 500
		delt_y = y_len

		sub_y_start_arr = np.arange(0, y_len, delt_y)
		sub_y_end_arr = sub_y_start_arr + delt_y
		sub_y_end_arr[-1] = min(sub_y_end_arr[-1], y_len)

		delt_x = split_pos - 0
		sub_x_start_arr = np.arange(0, split_pos, delt_x)
		sub_x_end_arr = sub_x_start_arr + delt_x
		sub_x_end_arr[-1] = min(sub_x_end_arr[-1], split_pos)

		for i, k1 in enumerate(sub_y_start_arr):
		sub_y_s = k1
		sub_y_e = sub_y_end_arr[i]

		sub_y_center = (sub_y_s + sub_y_e) / 2.

		for j, k2 in enumerate(sub_x_start_arr):
		sub_x_s = k2
		sub_x_e = sub_x_end_arr[j]
		# print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e)
		T1 = time.time()
		skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
		origin_sub = [sub_y_s, sub_x_s]
		sub_x_center = (sub_x_s + sub_x_e) / 2.

		sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center,
		origin=origin_sub,
		tar_spec=spec,
		band_start=tbstart, band_end=tbend,
		conf=conf1,
		flat_cube=flat_cube)

		spec_orders = sdp.compute_spec_orders()

		for k, v in spec_orders.items():
		img_s = v[0]
		origin_order_x = v[1]
		origin_order_y = v[2]
		ssImg = galsim.ImageF(img_s)
		ssImg.setOrigin(origin_order_x, origin_order_y)
		bounds = ssImg.bounds & fImg.bounds
		if bounds.area() == 0:
		continue
		fImg[bounds] = fImg[bounds] + ssImg[bounds]

		T2 = time.time()

		print('time: %s ms' % ((T2 - T1) * 1000))

		delt_x = x_len - split_pos
		sub_x_start_arr = np.arange(split_pos, x_len, delt_x)
		sub_x_end_arr = sub_x_start_arr + delt_x
		sub_x_end_arr[-1] = min(sub_x_end_arr[-1], x_len)

		for i, k1 in enumerate(sub_y_start_arr):
		sub_y_s = k1
		sub_y_e = sub_y_end_arr[i]

		sub_y_center = (sub_y_s + sub_y_e) / 2.

		for j, k2 in enumerate(sub_x_start_arr):
		sub_x_s = k2
		sub_x_e = sub_x_end_arr[j]
		# print(i,j,sub_y_s, sub_y_e,sub_x_s,sub_x_e)

		T1 = time.time()

		skyImg_sub = galsim.Image(skyImg.array[sub_y_s:sub_y_e, sub_x_s:sub_x_e])
		origin_sub = [sub_y_s, sub_x_s]
		sub_x_center = (sub_x_s + sub_x_e) / 2.

		sdp = SpecDisperser(orig_img=skyImg_sub, xcenter=sub_x_center, ycenter=sub_y_center,
		origin=origin_sub,
		tar_spec=spec,
		band_start=tbstart, band_end=tbend,
		conf=conf2,
		flat_cube=flat_cube)

		spec_orders = sdp.compute_spec_orders()

		for k, v in spec_orders.items():
		img_s = v[0]
		origin_order_x = v[1]
		origin_order_y = v[2]
		ssImg = galsim.ImageF(img_s)
		ssImg.setOrigin(origin_order_x, origin_order_y)
		bounds = ssImg.bounds & fImg.bounds
		if bounds.area() == 0:
		continue
		fImg[bounds] = fImg[bounds] + ssImg[bounds]
		T2 = time.time()

		print('time: %s ms' % ((T2 - T1) * 1000))

		if isAlongY == 1:
		fimg, tmx, tmy = rotate90(array_orig=fImg.array, xc=0, yc=0, isClockwise=0)
		else:
		fimg = fImg.array

		fimg = fimg * pixelSize * pixelSize

		return fimg

ObservationSim/MockObject/SpecDisperser/SpecDisperser.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -239,7 +239,7 @@ class SpecDisperser(object):
		# else:
		# beam_flat[k] = self.flat_cube[:, originOut_y + i, originOut_x + j]

		status = disperse.disperse_grism_object(self.thumb_img,
		status = disperse.disperse_grism_object(self.thumb_img.astype(np.float32),
		flat_index[nonz], yfrac_beam[nonz],
		sensitivity_beam[nonz],
		modelf, x0,

ObservationSim/MockObject/init.py

+1 −0

Original line number	Diff line number	Diff line
		@@ -4,5 +4,6 @@ from .CatalogBase import CatalogBase
		from .Quasar import Quasar
		from .Star import Star
		from .Stamp import Stamp
		from .FlatLED import FlatLED
		# from .SkybackgroundMap import *
		# from .CosmicRay import CosmicRay

ObservationSim/ObservationSim.py

+150 −7

File changed.

Preview size limit exceeded, changes collapsed.