diff --git a/.gitignore b/.gitignore index 1845922d7c86c84658453190839554dd92f753b1..e6e54466e20bc7a1dfebf86feb2a28ed8209d68b 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,6 @@ dist/* *.so *disperse.c *interp.c -!*libshao.so \ No newline at end of file +!*libshao.so +*.out +pnodes \ No newline at end of file diff --git a/ObservationSim/MockObject/CatalogBase.py b/ObservationSim/MockObject/CatalogBase.py index a65943e7971b0c1ba8453b9783c4cc5df69578b8..fda284966584ca9a7946307fa8a298eb4a38dd3c 100644 --- a/ObservationSim/MockObject/CatalogBase.py +++ b/ObservationSim/MockObject/CatalogBase.py @@ -23,56 +23,56 @@ class CatalogBase(metaclass=ABCMeta): @abstractmethod def load_norm_filt(self, obj): return None - + @staticmethod def initialize_param(): param = { - "star":-1, - "id":-1, - "ra":0, - "dec":0., - "ra_orig":0., - "dec_orig":0., - "z":0., - "sed_type":-1, - "model_tag":"unknown", - "mag_use_normal":100., - "theta":0., - "kappa":0., - "g1":0., - "g2":0., - "bfrac":0., - "av":0., - "redden":0., - "hlr_bulge":0., - "hlr_disk":0., - "ell_bulge":0., - "ell_disk":0., - "ell_tot":0., - "e1_disk":0., - "e2_disk":0., - "e1_bulge":0., - "e2_bulge":0., - "teff":0., - "logg":0., - "feh":0., - "DM":0., - "stellarMass":1., + "star": -1, + "id": -1, + "ra": 0, + "dec": 0., + "ra_orig": 0., + "dec_orig": 0., + "z": 0., + "sed_type": -1, + "model_tag": "unknown", + "mag_use_normal": 100., + "theta": 0., + "kappa": 0., + "g1": 0., + "g2": 0., + "bfrac": 0., + "av": 0., + "redden": 0., + "hlr_bulge": 0., + "hlr_disk": 0., + "ell_bulge": 0., + "ell_disk": 0., + "ell_tot": 0., + "e1_disk": 0., + "e2_disk": 0., + "e1_bulge": 0., + "e2_bulge": 0., + "teff": 0., + "logg": 0., + "feh": 0., + "DM": 0., + "stellarMass": 1., # C6 galaxies parameters - "e1":0., - "e2":0., - "bulgemass":0., - "diskmass":0., - "size":0., - "detA":0., - "type":0, - "veldisp":0., + "e1": 0., + "e2": 0., + "bulgemass": 0., + "diskmass": 0., + "size": 0., + "detA": 0., + "type": 0, + "veldisp": 0., "coeff": np.zeros(20), # Astrometry related - "pmra":0., - "pmdec":0., - "rv":0., - "parallax":1e-9 + "pmra": 0., + "pmdec": 0., + "rv": 0., + "parallax": 1e-9 } return param @@ -85,29 +85,34 @@ class CatalogBase(metaclass=ABCMeta): e1 = e_total * np.cos(phi + 2*rotation) e2 = e_total * np.sin(phi + 2*rotation) 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]), - eWave=np.ceil(norm_filt[norm_thr_rang_ids][-1][0])) + spectrum=sed, + norm_thr=norm_filt, + 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, @@ -117,4 +122,4 @@ class CatalogBase(metaclass=ABCMeta): return sed_photon, mag_csst, interFlux elif target_filt.survey_type == "spectroscopic": del sed_photon - return sed, mag_csst, interFlux \ No newline at end of file + return sed, mag_csst, interFlux diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 84de317653670b75b52eb1e842f762534da6160b..27d56380d98fb0b0872012396a7e1eaa1e94a2c2 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -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: @@ -46,7 +54,7 @@ class Galaxy(MockObject): if self.logger: self.logger.error(e) return -1 - + ratio = sub/full if not (ratio == -1 or (ratio != ratio)): nphotons = ratio * nphotons_tot @@ -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,14 +77,16 @@ 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) return final @@ -85,19 +97,20 @@ 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: self.logger.error(e) return 2, None - + # # [C6 TEST] # print('hlr_disk = %.4f, hlr_bulge = %.4f'%(self.hlr_disk, self.hlr_bulge)) # tracemalloc.start() big_galaxy = False - if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy + if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy big_galaxy = True # Set Galsim Parameters @@ -121,13 +134,15 @@ 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) - + # Get shear and deal with shear induced by field distortion if fd_shear: g1 += fd_shear.g1 @@ -149,14 +164,15 @@ class Galaxy(MockObject): nphotons = ratio * nphotons_tot else: continue - + # nphotons_sum += nphotons # # [C6 TEST] # 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 elif self.bfrac == 1: @@ -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) - if not big_galaxy: # Not apply PSF for very big galaxy + # 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,21 +203,22 @@ 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] is_updated = 1 chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) del stamp - + if is_updated == 0: # Return code 0: object photons missed this detector - print("obj %s missed"%(self.id)) + print("obj %s missed" % (self.id)) if self.logger: - self.logger.info("obj %s missed"%(self.id)) + self.logger.info("obj %s missed" % (self.id)) return 0, pos_shear - + # # [C6 TEST] # print("nphotons_sum = ", nphotons_sum) return 1, pos_shear @@ -208,9 +228,10 @@ class Galaxy(MockObject): if normFilter is not None: 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]), - eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0])) + norm_thr=normFilter, + 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 else: @@ -230,9 +251,8 @@ 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 + if self.hlr_disk > 3.0 or self.hlr_bulge > 3.0: # Very big galaxy big_galaxy = True if self.getMagFilter(filt) <= 15 and (not big_galaxy): @@ -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,13 +288,15 @@ 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) - + if self.bfrac == 0: gal = disk elif self.bfrac == 1: @@ -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,40 +323,43 @@ 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) star_p1.setOrigin(0, 0) origin_p1 = origin_star - xcenter_p1 = min(x_nominal,grating_split_pos_chip-1) - 0 + xcenter_p1 = min(x_nominal, grating_split_pos_chip-1) - 0 ycenter_p1 = y_nominal-0 sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, - ycenter=ycenter_p1, origin=origin_p1, - tar_spec=normalSED, - band_start=brange[0], band_end=brange[1], - conf=chip.sls_conf[0], - isAlongY=0, - flat_cube=flat_cube) + ycenter=ycenter_p1, origin=origin_p1, + tar_spec=normalSED, + band_start=brange[0], band_end=brange[1], + conf=chip.sls_conf[0], + isAlongY=0, + flat_cube=flat_cube) # self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1], psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, - local_wcs=chip_wcs_local, pos_img = pos_img) + 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] @@ -351,11 +378,11 @@ class Galaxy(MockObject): pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2], psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, - local_wcs=chip_wcs_local, pos_img = pos_img) + local_wcs=chip_wcs_local, pos_img=pos_img) del sdp_p1 del sdp_p2 - elif grating_split_pos_chip<=gal_origin[1]: + elif grating_split_pos_chip <= gal_origin[1]: sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, @@ -367,9 +394,9 @@ class Galaxy(MockObject): pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, - local_wcs=chip_wcs_local, pos_img = pos_img) + local_wcs=chip_wcs_local, pos_img=pos_img) del sdp - elif grating_split_pos_chip>=gal_end[1]: + elif grating_split_pos_chip >= gal_end[1]: sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, ycenter=y_nominal - 0, origin=origin_star, tar_spec=normalSED, @@ -381,7 +408,7 @@ class Galaxy(MockObject): pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], psf_model=psf_model, bandNo=i + 1, grating_split_pos=grating_split_pos, - local_wcs=chip_wcs_local, pos_img = pos_img) + local_wcs=chip_wcs_local, pos_img=pos_img) del sdp # print(self.y_nominal, starImg.center.y, starImg.ymin) @@ -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 diff --git a/ObservationSim/MockObject/Quasar.py b/ObservationSim/MockObject/Quasar.py index b1fae7fb306ff62c0ed94660bbccd596743f4925..1cfd99d1402cea6c9a8798644a78ad96abc958de 100755 --- a/ObservationSim/MockObject/Quasar.py +++ b/ObservationSim/MockObject/Quasar.py @@ -1,5 +1,6 @@ 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): ''' @@ -22,8 +29,9 @@ class Quasar(MockObject): norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 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)) + itype = objtypes[cosids == self.sed_type][0] + 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") @@ -33,49 +41,56 @@ class Quasar(MockObject): sed_data = sed_templates[self.sed_type] # redshift, intrinsic extinction sed_data = getObservedSED( - sedCat=sed_data, - redshift=self.z, - av=self.param['av'], + sedCat=sed_data, + redshift=self.z, + av=self.param['av'], redden=self.param['redden']) 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]), - eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0])) - sed_photon = np.array([sed_photon['WAVELENGTH'], sed_photon['FLUX']*sedNormFactor]).T + spectrum=sed_photon, + norm_thr=normFilter, + 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 # 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) - self.param['mag_%s'%target_filt.filter_type] = getABMAG( - interFlux=interFlux, + 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) # print('mag_use_normal = ', self.param['mag_use_normal']) # print('mag_%s = '%target_filt.filter_type, self.param['mag_%s'%target_filt.filter_type]) 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( - sedCat=sed_data, - redshift=self.z, - av=self.param['av'], + sedCat=sed_data, + redshift=self.z, + av=self.param['av'], redden=self.param['redden']) speci = interpolate.interp1d(sed_data[0], sed_data[1]) lamb = np.arange(2500, 10001 + 0.5, 0.5) 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 @@ -88,4 +103,4 @@ class Quasar(MockObject): qso = galsim.Gaussian(sigma=1.e-8, flux=1.) qso = qso.withFlux(flux) final = galsim.Convolve(psf, qso) - return final \ No newline at end of file + return final diff --git a/ObservationSim/MockObject/Star.py b/ObservationSim/MockObject/Star.py index c6510eb77a49d2c5e82ea3796f8fa8b2adf65e79..04fc88ea7a3e7369b4d98f74fdb5531616e6d5ec 100755 --- a/ObservationSim/MockObject/Star.py +++ b/ObservationSim/MockObject/Star.py @@ -1,5 +1,6 @@ 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 @@ -25,16 +29,18 @@ class Star(MockObject): star = star.withFlux(flux) final = galsim.Convolve(psf, star) return final - + 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) @@ -49,7 +55,7 @@ class Star(MockObject): print(e) self.logger.error(e) return -1 - + ratio = sub/full if not (ratio == -1 or (ratio != ratio)): @@ -61,4 +67,4 @@ class Star(MockObject): star = galsim.Convolve(psf, star) objs.append(star) final = galsim.Sum(objs) - return final \ No newline at end of file + return final diff --git a/ObservationSim/sim_steps/add_objects.py b/ObservationSim/sim_steps/add_objects.py index 942a2a0cd7d5e8bcd55c7dc6ffd92706f692611c..7b21e5d562c9398b83b33da4f20ba74a8325db46 100644 --- a/ObservationSim/sim_steps/add_objects.py +++ b/ObservationSim/sim_steps/add_objects.py @@ -10,6 +10,7 @@ from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSL from astropy.time import Time from datetime import datetime, timezone + def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # Get exposure time @@ -20,9 +21,12 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # Load catalogues if catalog is None: - self.chip_output.Log_error("Catalog interface class must be specified for SCIE-OBS") - raise ValueError("Catalog interface class must be specified for SCIE-OBS") - cat = catalog(config=self.overall_config, chip=chip, pointing=pointing, chip_output=self.chip_output, filt=filt) + self.chip_output.Log_error( + "Catalog interface class must be specified for SCIE-OBS") + raise ValueError( + "Catalog interface class must be specified for SCIE-OBS") + cat = catalog(config=self.overall_config, chip=chip, + pointing=pointing, chip_output=self.chip_output, filt=filt) # Prepare output file(s) for this chip # [NOTE] Headers of output .cat file may be updated by Catalog instance @@ -31,15 +35,18 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # Prepare the PSF model if self.overall_config["psf_setting"]["psf_model"] == "Gauss": - psf_model = PSFGauss(chip=chip, psfRa=self.overall_config["psf_setting"]["psf_rcont"]) + psf_model = PSFGauss( + chip=chip, psfRa=self.overall_config["psf_setting"]["psf_rcont"]) elif self.overall_config["psf_setting"]["psf_model"] == "Interp": if chip.survey_type == "spectroscopic": - psf_model = PSFInterpSLS(chip, filt, PSF_data_prefix=self.overall_config["psf_setting"]["psf_sls_dir"]) + psf_model = PSFInterpSLS( + chip, filt, PSF_data_prefix=self.overall_config["psf_setting"]["psf_sls_dir"]) else: - psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.overall_config["psf_setting"]["psf_pho_dir"]) + psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, + PSF_data_file=self.overall_config["psf_setting"]["psf_pho_dir"]) else: self.chip_output.Log_error("unrecognized PSF model type!!", flush=True) - + # Apply field distortion model if obs_param["field_dist"] == True: fd_model = FieldDistortion(chip=chip, img_rot=pointing.img_pa.deg) @@ -50,21 +57,23 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # Get the filter which will be used for magnitude cut for ifilt in range(len(self.all_filters)): temp_filter = self.all_filters[ifilt] - temp_filter.update_limit_saturation_mags(exptime=pointing.get_full_depth_exptime(temp_filter.filter_type), chip=chip) + temp_filter.update_limit_saturation_mags( + exptime=pointing.get_full_depth_exptime(temp_filter.filter_type), chip=chip) if temp_filter.filter_type.lower() == self.overall_config["obs_setting"]["cut_in_band"].lower(): cut_filter = temp_filter # Read in shear values from configuration file if the constant shear type is used if self.overall_config["shear_setting"]["shear_type"] == "constant": g1_field, g2_field, _ = get_shear_field(config=self.overall_config) - self.chip_output.Log_info("Use constant shear: g1=%.5f, g2=%.5f"%(g1_field, g2_field)) + self.chip_output.Log_info( + "Use constant shear: g1=%.5f, g2=%.5f" % (g1_field, g2_field)) # Get chip WCS if not hasattr(self, 'h_ext'): _, _ = self.prepare_headers(chip=chip, pointing=pointing) - - chip_wcs = galsim.FitsWCS(header = self.h_ext) - + + chip_wcs = galsim.FitsWCS(header=self.h_ext) + # Loop over objects nobj = len(cat.objs) missed_obj = 0 @@ -80,17 +89,19 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): try: sed_data = cat.load_sed(obj) norm_filt = cat.load_norm_filt(obj) - obj.sed, obj.param["mag_%s"%filt.filter_type.lower()], obj.param["flux_%s"%filt.filter_type.lower()] = cat.convert_sed( + obj.sed, obj.param["mag_%s" % filt.filter_type.lower()], obj.param["flux_%s" % filt.filter_type.lower()] = cat.convert_sed( mag=obj.param["mag_use_normal"], sed=sed_data, - target_filt=filt, + target_filt=filt, norm_filt=norm_filt, + mu=obj.mu ) - _, obj.param["mag_%s"%cut_filter.filter_type.lower()], obj.param["flux_%s"%cut_filter.filter_type.lower()] = cat.convert_sed( + _, obj.param["mag_%s" % cut_filter.filter_type.lower()], obj.param["flux_%s" % cut_filter.filter_type.lower()] = cat.convert_sed( mag=obj.param["mag_use_normal"], sed=sed_data, - target_filt=cut_filter, + target_filt=cut_filter, norm_filt=norm_filt, + mu=obj.mu ) except Exception as e: traceback.print_exc() @@ -102,16 +113,19 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): # Exclude very bright/dim objects (for now) if cut_filter.is_too_bright( - mag=obj.param["mag_%s"%self.overall_config["obs_setting"]["cut_in_band"].lower()], - margin=self.overall_config["obs_setting"]["mag_sat_margin"]): - self.chip_output.Log_info("obj %s too birght!! mag_%s = %.3f"%(obj.id, cut_filter.filter_type, obj.param["mag_%s"%self.overall_config["obs_setting"]["cut_in_band"].lower()])) + mag=obj.param["mag_%s" % + self.overall_config["obs_setting"]["cut_in_band"].lower()], + margin=self.overall_config["obs_setting"]["mag_sat_margin"]): + self.chip_output.Log_info("obj %s too birght!! mag_%s = %.3f" % ( + obj.id, cut_filter.filter_type, obj.param["mag_%s" % self.overall_config["obs_setting"]["cut_in_band"].lower()])) bright_obj += 1 obj.unload_SED() continue if filt.is_too_dim( - mag=obj.getMagFilter(filt), - margin=self.overall_config["obs_setting"]["mag_lim_margin"]): - self.chip_output.Log_info("obj %s too dim!! mag_%s = %.3f"%(obj.id, filt.filter_type, obj.getMagFilter(filt))) + mag=obj.getMagFilter(filt), + margin=self.overall_config["obs_setting"]["mag_lim_margin"]): + self.chip_output.Log_info("obj %s too dim!! mag_%s = %.3f" % ( + obj.id, filt.filter_type, obj.getMagFilter(filt))) dim_obj += 1 obj.unload_SED() continue @@ -130,14 +144,16 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): raise ValueError("Unknown shear input") # Get position of object on the focal plane - pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext) + pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS( + img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext) # [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane # Otherwise they will be considered missed objects # if pos_img.x == -1 or pos_img.y == -1 or (not chip.isContainObj(x_image=pos_img.x, y_image=pos_img.y, margin=0.)): if pos_img.x == -1 or pos_img.y == -1: - self.chip_output.Log_info('obj_ra = %.6f, obj_dec = %.6f, obj_ra_orig = %.6f, obj_dec_orig = %.6f'%(obj.ra, obj.dec, obj.ra_orig, obj.dec_orig)) - self.chip_output.Log_error("Objected missed: %s"%(obj.id)) + self.chip_output.Log_info('obj_ra = %.6f, obj_dec = %.6f, obj_ra_orig = %.6f, obj_dec_orig = %.6f' % ( + obj.ra, obj.dec, obj.ra_orig, obj.dec_orig)) + self.chip_output.Log_error("Objected missed: %s" % (obj.id)) missed_obj += 1 obj.unload_SED() continue @@ -146,31 +162,32 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): try: if self.overall_config["run_option"]["out_cat_only"]: isUpdated = True - obj.real_pos = obj.getRealPos(chip.img, global_x=obj.posImg.x, global_y=obj.posImg.y, img_real_wcs=obj.chip_wcs) + obj.real_pos = obj.getRealPos( + chip.img, global_x=obj.posImg.x, global_y=obj.posImg.y, img_real_wcs=obj.chip_wcs) pos_shear = 0. elif chip.survey_type == "photometric" and not self.overall_config["run_option"]["out_cat_only"]: isUpdated, pos_shear = obj.drawObj_multiband( tel=tel, - pos_img=pos_img, - psf_model=psf_model, - bandpass_list=filt.bandpass_sub_list, - filt=filt, - chip=chip, - g1=obj.g1, - g2=obj.g2, + pos_img=pos_img, + psf_model=psf_model, + bandpass_list=filt.bandpass_sub_list, + filt=filt, + chip=chip, + g1=obj.g1, + g2=obj.g2, exptime=exptime, fd_shear=fd_shear) elif chip.survey_type == "spectroscopic" and not self.overall_config["run_option"]["out_cat_only"]: isUpdated, pos_shear = obj.drawObj_slitless( - tel=tel, - pos_img=pos_img, - psf_model=psf_model, - bandpass_list=filt.bandpass_sub_list, - filt=filt, - chip=chip, - g1=obj.g1, - g2=obj.g2, + tel=tel, + pos_img=pos_img, + psf_model=psf_model, + bandpass_list=filt.bandpass_sub_list, + filt=filt, + chip=chip, + g1=obj.g1, + g2=obj.g2, exptime=exptime, normFilter=norm_filt, fd_shear=fd_shear) @@ -181,9 +198,10 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): pass elif isUpdated == 0: missed_obj += 1 - self.chip_output.Log_error("Objected missed: %s"%(obj.id)) + self.chip_output.Log_error("Objected missed: %s" % (obj.id)) else: - self.chip_output.Log_error("Draw error, object omitted: %s"%(obj.id)) + self.chip_output.Log_error( + "Draw error, object omitted: %s" % (obj.id)) continue except Exception as e: traceback.print_exc() @@ -196,38 +214,47 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): del psf_model gc.collect() - self.chip_output.Log_info("Running checkpoint #1 (Object rendering finished): pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) )) + self.chip_output.Log_info("Running checkpoint #1 (Object rendering finished): pointing-%d chip-%d pid-%d memory-%6.2fGB" % + (pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024))) - self.chip_output.Log_info("# objects that are too bright %d out of %d"%(bright_obj, nobj)) - self.chip_output.Log_info("# objects that are too dim %d out of %d"%(dim_obj, nobj)) - self.chip_output.Log_info("# objects that are missed %d out of %d"%(missed_obj, nobj)) + self.chip_output.Log_info( + "# objects that are too bright %d out of %d" % (bright_obj, nobj)) + self.chip_output.Log_info( + "# objects that are too dim %d out of %d" % (dim_obj, nobj)) + self.chip_output.Log_info( + "# objects that are missed %d out of %d" % (missed_obj, nobj)) # Apply flat fielding (with shutter effects) flat_normal = np.ones_like(chip.img.array) if obs_param["flat_fielding"] == True: - flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array) + flat_normal = flat_normal * chip.flat_img.array / \ + np.mean(chip.flat_img.array) if obs_param["shutter_effect"] == True: flat_normal = flat_normal * chip.shutter_img flat_normal = np.array(flat_normal, dtype='float32') - self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True]) + self.updateHeaderInfo(header_flag='ext', keys=[ + 'SHTSTAT'], values=[True]) else: - self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN1','SHTCLOS0'], values = [True,self.h_ext['SHTCLOS1'],self.h_ext['SHTOPEN0']]) + self.updateHeaderInfo(header_flag='ext', keys=['SHTSTAT', 'SHTOPEN1', 'SHTCLOS0'], values=[ + True, self.h_ext['SHTCLOS1'], self.h_ext['SHTOPEN0']]) chip.img *= flat_normal del flat_normal - # renew header info datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) datetime_obs = datetime_obs.replace(tzinfo=timezone.utc) t_obs = Time(datetime_obs) - - ##ccd刷新2s,等待0.s，开始曝光 + + # ccd刷新2s,等待0.s，开始曝光 t_obs_renew = Time(t_obs.mjd - (2.+0.) / 86400., format="mjd") - t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) - self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]]) + t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp( + t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1)) + self.updateHeaderInfo(header_flag='prim', keys=[ + 'DATE-OBS'], values=[t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]]) - #dark time : 曝光时间+刷新后等带时间0.s+关快门后读出前等待0.s - self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.+0.+pointing.exp_time]) + # dark time : 曝光时间+刷新后等带时间0.s+关快门后读出前等待0.s + self.updateHeaderInfo(header_flag='ext', keys=[ + 'DARKTIME'], values=[0.+0.+pointing.exp_time]) - return chip, filt, tel, pointing \ No newline at end of file + return chip, filt, tel, pointing diff --git a/setup.py b/setup.py index 4aa75806caefb2a0527773a34e6a2a9ce47bfc78..356d247f73fafdb57583d8a45a07205c75c45de5 100644 --- a/setup.py +++ b/setup.py @@ -76,7 +76,7 @@ with open("requirements.txt", "r") as f: ] setup(name='csst_msc_sim', - version='3.0.0', + version='3.0.0rc', packages=find_packages(), # install_requires=[ # # 'numpy>=1.18.5',