1. add lensing magnification (53592d2e) · Commits · csst-sims / csst_msc_sim

.gitignore

+3 −1

Original line number	Diff line number	Diff line
		@@ -6,3 +6,5 @@ dist/*
		*disperse.c
		*interp.c
		!*libshao.so
		*.out
		pnodes
		No newline at end of file

ObservationSim/MockObject/CatalogBase.py

+62 −57

Original line number	Diff line number	Diff line
		@@ -87,27 +87,32 @@ class CatalogBase(metaclass=ABCMeta):
		return e1, e2, e_total

		@staticmethod
		def convert_sed(mag, sed, target_filt, norm_filt=None):
		def convert_sed(mag, sed, target_filt, norm_filt=None, mu=1.):
		bandpass = target_filt.bandpass_full

		if norm_filt is not None:
		norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001
		else:
		norm_filt = Table(
		np.array(np.array([bandpass.wave_list*10.0, bandpass.func(bandpass.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY'])
		np.array(np.array([bandpass.wave_list*10.0, bandpass.func(
		bandpass.wave_list)])).T, names=(['WAVELENGTH', 'SENSITIVITY'])
		)
		norm_thr_rang_ids = norm_filt['SENSITIVITY'] > 0.001

		sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=mag,
		spectrum=sed,
		norm_thr=norm_filt,
		sWave=np.floor(norm_filt[norm_thr_rang_ids][0][0]),
		sWave=np.floor(
		norm_filt[norm_thr_rang_ids][0][0]),
		eWave=np.ceil(norm_filt[norm_thr_rang_ids][-1][0]))
		sed_photon = copy.copy(sed)
		sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
		sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
		sed_photon = np.array(
		[sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
		sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
		sed_photon[:, 1] * mu), interpolant='nearest')
		# Get magnitude
		sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
		sed_photon = galsim.SED(sed_photon, wave_type='A',
		flux_type='1', fast=False)
		interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=bandpass)
		mag_csst = getABMAG(
		interFlux=interFlux,

ObservationSim/MockObject/Galaxy.py

+84 −54

Original line number	Diff line number	Diff line
		@@ -8,6 +8,7 @@ from ObservationSim.MockObject.MockObject import MockObject

		# import tracemalloc


		class Galaxy(MockObject):
		def __init__(self, param, logger=None):
		super().__init__(param, logger=logger)
		@@ -16,6 +17,11 @@ class Galaxy(MockObject):
		self.disk_sersic_idx = 1.
		if not hasattr(self, "bulge_sersic_idx"):
		self.bulge_sersic_idx = 4.
		if not hasattr(self, "mu"):
		if hasattr(self, "detA"):
		self.mu = 1./self.detA
		else:
		self.mu = 1.

		def unload_SED(self):
		"""(Test) free up SED memory
		@@ -24,14 +30,16 @@ class Galaxy(MockObject):

		def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150., fd_shear=None):
		if len(psf_list) != len(bandpass_list):
		raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.")
		raise ValueError(
		"!!!The number of PSF profiles and the number of bandpasses must be equal.")
		objs = []
		if nphotons_tot == None:
		nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
		# print("nphotons_tot = ", nphotons_tot)

		try:
		full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
		full = integrate_sed_bandpass(
		sed=self.sed, bandpass=filt.bandpass_full)
		except Exception as e:
		print(e)
		if self.logger:
		@@ -54,10 +62,12 @@ class Galaxy(MockObject):
		return -1

		psf = psf_list[i]
		disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0)
		disk = galsim.Sersic(n=self.disk_sersic_idx,
		half_light_radius=self.hlr_disk, flux=1.0)
		disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
		disk = disk.shear(disk_shape)
		bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0)
		bulge = galsim.Sersic(n=self.bulge_sersic_idx,
		half_light_radius=self.hlr_bulge, flux=1.0)
		bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
		bulge = bulge.shear(bulge_shape)

		@@ -67,13 +77,15 @@ class Galaxy(MockObject):
		gal = bulge
		else:
		gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk
		gal = gal.withFlux(nphotons)
		if fd_shear is not None:
		g1 += fd_shear.g1
		g2 += fd_shear.g2
		gal_shear = galsim.Shear(g1=g1, g2=g2)
		gal = gal.shear(gal_shear)
		# Magnification
		gal = gal.magnify(self.mu)
		gal = galsim.Convolve(psf, gal)
		gal = gal.withFlux(nphotons)

		objs.append(gal)
		final = galsim.Sum(objs)
		@@ -85,7 +97,8 @@ class Galaxy(MockObject):
		# print("nphotons_tot = ", nphotons_tot)

		try:
		full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
		full = integrate_sed_bandpass(
		sed=self.sed, bandpass=filt.bandpass_full)
		except Exception as e:
		print(e)
		if self.logger:
		@@ -121,10 +134,12 @@ class Galaxy(MockObject):
		is_updated = 0

		# Model the galaxy as disk + bulge
		disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
		disk = galsim.Sersic(
		n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
		disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
		disk = disk.shear(disk_shape)
		bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
		bulge = galsim.Sersic(
		n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
		bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
		bulge = bulge.shear(bulge_shape)

		@@ -155,7 +170,8 @@ class Galaxy(MockObject):
		# print("nphotons_sub-band_%d = %.2f"%(i, nphotons))

		# Get PSF model
		psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
		psf, pos_shear = psf_model.get_PSF(
		chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)

		if self.bfrac == 0:
		gal_temp = disk
		@@ -164,10 +180,13 @@ class Galaxy(MockObject):
		else:
		gal_temp = self.bfrac * bulge + (1.0 - self.bfrac) * disk
		gal_temp = gal_temp.shear(gal_shear)
		gal_temp = gal_temp.withFlux(nphotons)
		# Magnification
		gal_temp = gal_temp.magnify(self.mu)
		if not big_galaxy: # Not apply PSF for very big galaxy
		gal_temp = galsim.Convolve(psf, gal_temp)

		gal_temp = gal_temp.withFlux(nphotons)

		if i == 0:
		gal = gal_temp
		else:
		@@ -184,7 +203,8 @@ class Galaxy(MockObject):
		# ERROR happens
		return 2, pos_shear
		stamp.setCenter(x_nominal, y_nominal)
		bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)
		bounds = stamp.bounds & galsim.BoundsI(
		0, chip.npix_x - 1, 0, chip.npix_y - 1)
		if bounds.area() > 0:
		chip.img.setOrigin(0, 0)
		chip.img[bounds] += stamp[bounds]
		@@ -209,7 +229,8 @@ class Galaxy(MockObject):
		norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
		sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
		norm_thr=normFilter,
		sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
		sWave=np.floor(
		normFilter[norm_thr_rang_ids][0][0]),
		eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
		if sedNormFactor == 0:
		return 2, None
		@@ -230,7 +251,6 @@ class Galaxy(MockObject):

		chip_wcs_local = self.chip_wcs.local(self.real_pos)


		big_galaxy = False
		if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy
		big_galaxy = True
		@@ -244,7 +264,8 @@ class Galaxy(MockObject):

		flat_cube = chip.flat_cube

		xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
		xOrderSigPlus = {'A': 1.3909419820029296, 'B': 1.4760376591236062,
		'C': 4.035447379743442, 'D': 5.5684364343742825, 'E': 16.260021029735388}
		grating_split_pos_chip = 0 + grating_split_pos

		branges = np.zeros([len(bandpass_list), 2])
		@@ -267,10 +288,12 @@ class Galaxy(MockObject):
		brange = branges[i]

		# psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
		disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
		disk = galsim.Sersic(
		n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0, gsparams=gsp)
		disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
		disk = disk.shear(disk_shape)
		bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
		bulge = galsim.Sersic(
		n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0, gsparams=gsp)
		bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
		bulge = bulge.shear(bulge_shape)

		@@ -286,12 +309,13 @@ class Galaxy(MockObject):
		# kfrac = np.random.random()*(1.0 - self.bfrac)
		# gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots

		gal = gal.withFlux(tel.pupil_area * exptime)
		if fd_shear:
		g1 += fd_shear.g1
		g2 += fd_shear.g2
		gal_shear = galsim.Shear(g1=g1, g2=g2)
		gal = gal.shear(gal_shear)
		gal = gal.magnify(self.mu)
		gal = gal.withFlux(tel.pupil_area * exptime)
		# gal = galsim.Convolve(psf, gal)

		# if not big_galaxy: # Not apply PSF for very big galaxy
		@@ -299,17 +323,19 @@ class Galaxy(MockObject):
		# # if fd_shear is not None:
		# # gal = gal.shear(fd_shear)

		starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset,method = 'real_space')
		starImg = gal.drawImage(
		wcs=chip_wcs_local, offset=offset, method='real_space')

		origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
		x_nominal - (starImg.center.x - starImg.xmin)]
		starImg.setOrigin(0, 0)
		gal_origin = [origin_star[0], origin_star[1]]
		gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1]
		gal_end = [origin_star[0] + starImg.array.shape[0] -
		1, origin_star[1] + starImg.array.shape[1] - 1]

		if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
		subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
		## part img disperse
		# part img disperse

		subImg_p1 = starImg.array[:, 0:subSlitPos]
		star_p1 = galsim.Image(subImg_p1)
		@@ -332,7 +358,8 @@ class Galaxy(MockObject):
		grating_split_pos=grating_split_pos,
		local_wcs=chip_wcs_local, pos_img=pos_img)

		subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
		subImg_p2 = starImg.array[:,
		subSlitPos+1:starImg.array.shape[1]]
		star_p2 = galsim.Image(subImg_p2)
		star_p2.setOrigin(0, 0)
		origin_p2 = [origin_star[0], grating_split_pos_chip]
		@@ -391,11 +418,13 @@ class Galaxy(MockObject):
		def getGSObj(self, psf, g1=0, g2=0, flux=None, filt=None, tel=None, exptime=150.):
		if flux == None:
		flux = self.getElectronFluxFilt(filt, tel, exptime)
		disk = galsim.Sersic(n=self.disk_sersic_idx, half_light_radius=self.hlr_disk, flux=1.0)
		disk = galsim.Sersic(n=self.disk_sersic_idx,
		half_light_radius=self.hlr_disk, flux=1.0)
		disk_shape = galsim.Shear(g1=self.e1_disk, g2=self.e2_disk)
		disk = disk.shear(disk_shape)

		bulge = galsim.Sersic(n=self.bulge_sersic_idx, half_light_radius=self.hlr_bulge, flux=1.0)
		bulge = galsim.Sersic(n=self.bulge_sersic_idx,
		half_light_radius=self.hlr_bulge, flux=1.0)
		bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge)
		bulge = bulge.shear(bulge_shape)

		@@ -407,5 +436,6 @@ class Galaxy(MockObject):
		return final

		def getObservedEll(self, g1=0, g2=0):
		e1_obs, e2_obs, e_obs, theta = eObs(self.e1_total, self.e2_total, g1, g2)
		e1_obs, e2_obs, e_obs, theta = eObs(
		self.e1_total, self.e2_total, g1, g2)
		return self.e1_total, self.e2_total, g1, g2, e1_obs, e2_obs

ObservationSim/MockObject/Quasar.py

+40 −25

Original line number	Diff line number	Diff line
		import galsim
		import os, sys
		import os
		import sys
		import numpy as np
		import astropy.constants as cons
		from astropy.table import Table
		@@ -8,9 +9,15 @@ from scipy import interpolate
		from ObservationSim.MockObject.MockObject import MockObject
		from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG


		class Quasar(MockObject):
		def __init__(self, param, logger=None):
		super().__init__(param, logger=logger)
		if not hasattr(self, "mu"):
		if hasattr(self, "detA"):
		self.mu = 1./self.detA
		else:
		self.mu = 1.

		def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None):
		'''
		@@ -23,7 +30,8 @@ class Quasar(MockObject):
		if sed_templates is None:
		# Read SED data directly
		itype = objtypes[cosids == self.sed_type][0]
		sed_file = os.path.join(sed_path, itype + "_ID%s.sed"%(self.sed_type))
		sed_file = os.path.join(
		sed_path, itype + "_ID%s.sed" % (self.sed_type))
		if not os.path.exists(sed_file):
		raise ValueError("!!! No SED found.")
		sed_data = Table.read(sed_file, format="ascii")
		@@ -40,19 +48,25 @@ class Quasar(MockObject):
		wave, flux = sed_data[0], sed_data[1]

		flux_photon = flux * (wave / (cons.h.value * cons.c.value)) * 1e-13
		sed_photon = Table(np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
		sed_photon = Table(
		np.array([wave, flux_photon]).T, names=('WAVELENGTH', 'FLUX'))
		# Get scaling factor for SED
		sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'],
		spectrum=sed_photon,
		norm_thr=normFilter,
		sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
		sWave=np.floor(
		normFilter[norm_thr_rang_ids][0][0]),
		eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
		sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
		sed_photon = np.array(
		[sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T
		# Convert to galsim.SED object
		spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
		self.sed = galsim.SED(spec, wave_type='A', flux_type='1', fast=False)
		spec = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
		sed_photon[:, 1]), interpolant='nearest')
		self.sed = galsim.SED(spec, wave_type='A',
		flux_type='1', fast=False)
		# Get magnitude
		interFlux = integrate_sed_bandpass(sed=self.sed, bandpass=target_filt.bandpass_full)
		interFlux = integrate_sed_bandpass(
		sed=self.sed, bandpass=target_filt.bandpass_full)
		self.param['mag_%s' % target_filt.filter_type] = getABMAG(
		interFlux=interFlux,
		bandpass=target_filt.bandpass_full)
		@@ -61,7 +75,8 @@ class Quasar(MockObject):

		elif survey_type == "spectroscopic":
		if sed_templates is None:
		self.sedPhotons(sed_path=sed_path, cosids=cosids, objtypes=objtypes)
		self.sedPhotons(sed_path=sed_path,
		cosids=cosids, objtypes=objtypes)
		else:
		sed_data = sed_templates[self.sed_type]
		sed_data = getObservedSED(
		@@ -74,8 +89,8 @@ class Quasar(MockObject):
		y = speci(lamb)
		# erg/s/cm2/A --> photo/s/m2/A
		all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
		self.sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))

		self.sed = Table(
		np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))

		def unload_SED(self):
		"""(Test) free up SED memory

ObservationSim/MockObject/Star.py

+12 −6

Original line number	Diff line number	Diff line
		import galsim
		import os, sys
		import os
		import sys
		import numpy as np
		import astropy.constants as cons
		from astropy.table import Table
		@@ -8,9 +9,12 @@ from scipy import interpolate
		from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG, tag_sed
		from ObservationSim.MockObject.MockObject import MockObject


		class Star(MockObject):
		def __init__(self, param, logger=None):
		super().__init__(param, logger=logger)
		if not hasattr(self, "mu"):
		self.mu = 1.

		def unload_SED(self):
		"""(Test) free up SED memory
		@@ -28,13 +32,15 @@ class Star(MockObject):

		def getGSObj_multiband(self, tel, psf_list, bandpass_list, filt, nphotons_tot=None, g1=0, g2=0, exptime=150.):
		if len(psf_list) != len(bandpass_list):
		raise ValueError("!!!The number of PSF profiles and the number of bandpasses must be equal.")
		raise ValueError(
		"!!!The number of PSF profiles and the number of bandpasses must be equal.")
		objs = []
		if nphotons_tot == None:
		nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)

		try:
		full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
		full = integrate_sed_bandpass(
		sed=self.sed, bandpass=filt.bandpass_full)
		except Exception as e:
		print(e)
		self.logger.error(e)