From 1aa4f5339045dae11b99fd09e67aedf9db731741 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Tue, 23 Jul 2024 08:01:37 +0800 Subject: [PATCH 1/5] 1. add option for galaxy milky way extinction in overall_config 2. bug fix: planck map nested ordering --- catalog/C9_Catalog.py | 21 +++++++++++++------- config/config_overall.yaml | 10 +++++++--- observation_sim/mock_objects/ExtinctionMW.py | 5 +++-- 3 files changed, 24 insertions(+), 12 deletions(-) diff --git a/catalog/C9_Catalog.py b/catalog/C9_Catalog.py index 978d74d..b028bdf 100644 --- a/catalog/C9_Catalog.py +++ b/catalog/C9_Catalog.py @@ -111,10 +111,14 @@ class Catalog(CatalogBase): self.max_size = 0. # [TODO] Milky Way extinction - self.mw_ext = ExtinctionMW() - self.mw_ext.init_ext_model(model_name="odonnell") - self.mw_ext.load_Planck_ext( - file_path="/public/home/fangyuedong/project/ext_maps/planck/HFI_CompMap_ThermalDustModel_2048_R1.20.fits") + if "enable_mw_ext_gal" in config["catalog_options"] and config["catalog_options"]["enable_mw_ext_gal"]: + if "planck_ebv_map" not in config["catalog_options"]: + raise ValueError( + "Planck dust map must be given to enable Milky Way extinction calculation for galaxies.") + self.mw_ext = ExtinctionMW() + self.mw_ext.init_ext_model(model_name="odonnell") + self.mw_ext.load_Planck_ext( + file_path=config["catalog_options"]["planck_ebv_map"]) if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: # Get the cloest star catalog file @@ -255,7 +259,10 @@ class Catalog(CatalogBase): ) # [TODO] get Milky Way extinction AVs - MW_Av_arr = self.mw_ext.Av_from_Planck(ra=ra_arr, dec=dec_arr) + if "enable_mw_ext_gal" in self.config["catalog_options"] and self.config["catalog_options"]["enable_mw_ext_gal"]: + MW_Av_arr = self.mw_ext.Av_from_Planck(ra=ra_arr, dec=dec_arr) + else: + MW_Av_arr = np.ones(len(ra_arr)) for igals in range(ngals): # # (TEST) @@ -544,8 +551,8 @@ class Catalog(CatalogBase): y = speci(lamb) # [TODO] Apply Milky Way extinction - if obj.type != 'star': - self.mw_ext.apply_extinction(y, Av=obj.mw_Av) + if obj.type != 'star' and ("enable_mw_ext_gal" in self.config["catalog_options"] and self.config["catalog_options"]["enable_mw_ext_gal"]): + y = self.mw_ext.apply_extinction(y, Av=obj.mw_Av) # erg/s/cm2/A --> photon/s/m2/A all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 diff --git a/config/config_overall.yaml b/config/config_overall.yaml index 5027e49..d22dc68 100644 --- a/config/config_overall.yaml +++ b/config/config_overall.yaml @@ -11,7 +11,7 @@ # can add some of the command-line arguments here as well; # ok to pass either way or both, as long as they are consistent work_dir: "/public/home/fangyuedong/project/workplace/" -run_name: "ext_test" +run_name: "ext_on" # Project cycle and run counter are used to name the outputs project_cycle: 9 @@ -44,11 +44,15 @@ catalog_options: star_only: NO # Only simulate galaxies? - galaxy_only: YES + galaxy_only: NO # rotate galaxy ellipticity rotateEll: 0. # [degree] + # Whether to apply milky way extinction to galaxies + enable_mw_ext_gal: YES + planck_ebv_map: "/public/home/fangyuedong/project/ext_maps/planck/HFI_CompMap_ThermalDustModel_2048_R1.20.fits" + ############################################### # Observation setting ############################################### @@ -68,7 +72,7 @@ obs_setting: run_pointings: [0, 1, 2, 3, 4] # Whether to enable astrometric modeling - enable_astrometric_model: True + enable_astrometric_model: YES # Cut by saturation magnitude in which band? cut_in_band: "z" diff --git a/observation_sim/mock_objects/ExtinctionMW.py b/observation_sim/mock_objects/ExtinctionMW.py index a98106a..a754deb 100644 --- a/observation_sim/mock_objects/ExtinctionMW.py +++ b/observation_sim/mock_objects/ExtinctionMW.py @@ -13,7 +13,7 @@ class ExtinctionMW(object): @staticmethod def radec2pix(ra, dec, NSIDE=2048): - return hp.pixelfunc.ang2pix(nside=NSIDE, theta=ra, phi=dec, lonlat=True) + return hp.pixelfunc.ang2pix(nside=NSIDE, theta=ra, phi=dec, nested=True, lonlat=True) def init_ext_model(self, model_name="odonnell", Av=1.0, Rv=3.1, lamb=None): self.model_name = model_name @@ -79,7 +79,8 @@ class ExtinctionMW(object): u.degree, frame='icrs').galactic l, b = c.l.radian, c.b.radian NSIDE = hp.pixelfunc.get_nside(self.ebv_planck) - pix = hp.pixelfunc.ang2pix(nside=NSIDE, theta=np.pi/2. - b, phi=l) + pix = hp.pixelfunc.ang2pix( + nside=NSIDE, theta=np.pi/2. - b, phi=l, nested=True) return self.ebv_planck[pix] * self.Rv def apply_extinction(self, spec, Av=1.0): -- GitLab From b8d1624727f2e38f5eee85db01fc80cb241ff0c9 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Wed, 24 Jul 2024 21:25:02 +0800 Subject: [PATCH 2/5] fix a typo in add_poisson_and_dark --- observation_sim/sim_steps/add_pattern_noise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/observation_sim/sim_steps/add_pattern_noise.py b/observation_sim/sim_steps/add_pattern_noise.py index dfe90d0..f92cf36 100644 --- a/observation_sim/sim_steps/add_pattern_noise.py +++ b/observation_sim/sim_steps/add_pattern_noise.py @@ -27,7 +27,7 @@ def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param): InputDark=None) else: chip.img, _ = chip_utils.add_poisson(img=chip.img, - chip=self, + chip=chip, exptime=exptime, poisson_noise=chip.poisson_noise, dark_noise=0.) -- GitLab From e66531a393106988a678adc7bb945599807e83e2 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Mon, 29 Jul 2024 06:54:42 +0800 Subject: [PATCH 3/5] bug fix C9_Catalog.py and ExtinctionMW.py --- catalog/C9_Catalog.py | 2 +- observation_sim/mock_objects/ExtinctionMW.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/catalog/C9_Catalog.py b/catalog/C9_Catalog.py index b028bdf..590dda6 100644 --- a/catalog/C9_Catalog.py +++ b/catalog/C9_Catalog.py @@ -262,7 +262,7 @@ class Catalog(CatalogBase): if "enable_mw_ext_gal" in self.config["catalog_options"] and self.config["catalog_options"]["enable_mw_ext_gal"]: MW_Av_arr = self.mw_ext.Av_from_Planck(ra=ra_arr, dec=dec_arr) else: - MW_Av_arr = np.ones(len(ra_arr)) + MW_Av_arr = np.zeros(len(ra_arr)) for igals in range(ngals): # # (TEST) diff --git a/observation_sim/mock_objects/ExtinctionMW.py b/observation_sim/mock_objects/ExtinctionMW.py index a754deb..85ae692 100644 --- a/observation_sim/mock_objects/ExtinctionMW.py +++ b/observation_sim/mock_objects/ExtinctionMW.py @@ -13,7 +13,7 @@ class ExtinctionMW(object): @staticmethod def radec2pix(ra, dec, NSIDE=2048): - return hp.pixelfunc.ang2pix(nside=NSIDE, theta=ra, phi=dec, nested=True, lonlat=True) + return hp.pixelfunc.ang2pix(nside=NSIDE, theta=ra, phi=dec, nest=True, lonlat=True) def init_ext_model(self, model_name="odonnell", Av=1.0, Rv=3.1, lamb=None): self.model_name = model_name @@ -80,7 +80,7 @@ class ExtinctionMW(object): l, b = c.l.radian, c.b.radian NSIDE = hp.pixelfunc.get_nside(self.ebv_planck) pix = hp.pixelfunc.ang2pix( - nside=NSIDE, theta=np.pi/2. - b, phi=l, nested=True) + nside=NSIDE, theta=np.pi/2. - b, phi=l, nest=True) return self.ebv_planck[pix] * self.Rv def apply_extinction(self, spec, Av=1.0): @@ -91,6 +91,6 @@ class ExtinctionMW(object): raise ValueError( "Need to initialize the extinction model (init_ext_model) first") scale = 10**(-.4 * self.ext * Av) - print("scale = ", scale) + # print("scale = ", scale) spec *= scale return spec -- GitLab From 36106f04b9468e914a43550f49dc7bc2f608ede6 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Tue, 30 Jul 2024 03:15:54 +0800 Subject: [PATCH 4/5] change bluest and reddest end points (to 2000 and 11000 respectively) for sub bandpasses. --- observation_sim/instruments/data/filters/fgs_sub.list | 2 +- observation_sim/instruments/data/filters/g_sub.list | 4 ++-- observation_sim/instruments/data/filters/i_sub.list | 4 ++-- observation_sim/instruments/data/filters/nuv_sub.list | 4 ++-- observation_sim/instruments/data/filters/r_sub.list | 4 ++-- observation_sim/instruments/data/filters/u_sub.list | 4 ++-- observation_sim/instruments/data/filters/y_sub.list | 4 ++-- observation_sim/instruments/data/filters/z_sub.list | 4 ++-- 8 files changed, 15 insertions(+), 15 deletions(-) diff --git a/observation_sim/instruments/data/filters/fgs_sub.list b/observation_sim/instruments/data/filters/fgs_sub.list index 5a826ee..c02a44f 100644 --- a/observation_sim/instruments/data/filters/fgs_sub.list +++ b/observation_sim/instruments/data/filters/fgs_sub.list @@ -1,4 +1,4 @@ -3000 +2000 4500 4750 5000 diff --git a/observation_sim/instruments/data/filters/g_sub.list b/observation_sim/instruments/data/filters/g_sub.list index f5dda57..1384ade 100644 --- a/observation_sim/instruments/data/filters/g_sub.list +++ b/observation_sim/instruments/data/filters/g_sub.list @@ -1,4 +1,4 @@ -3800 +2000 4217 4432 4631 @@ -6,4 +6,4 @@ 5002 5179 5354 -5799 \ No newline at end of file +11000 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/i_sub.list b/observation_sim/instruments/data/filters/i_sub.list index 715f74a..0449403 100644 --- a/observation_sim/instruments/data/filters/i_sub.list +++ b/observation_sim/instruments/data/filters/i_sub.list @@ -1,4 +1,4 @@ -6600 +2000 7061 7255 7448 @@ -6,4 +6,4 @@ 7833 8027 8226 -8999 \ No newline at end of file +11000 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/nuv_sub.list b/observation_sim/instruments/data/filters/nuv_sub.list index c89b5dc..b6c60c5 100644 --- a/observation_sim/instruments/data/filters/nuv_sub.list +++ b/observation_sim/instruments/data/filters/nuv_sub.list @@ -1,4 +1,4 @@ -2513 +2000 2621 2716 2805 @@ -6,4 +6,4 @@ 2969 3050 3132 -3499 \ No newline at end of file +11000 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/r_sub.list b/observation_sim/instruments/data/filters/r_sub.list index 0abf191..ec0852e 100644 --- a/observation_sim/instruments/data/filters/r_sub.list +++ b/observation_sim/instruments/data/filters/r_sub.list @@ -1,4 +1,4 @@ -5100 +2000 5642 5821 6001 @@ -6,4 +6,4 @@ 6363 6547 6735 -7199 \ No newline at end of file +11000 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/u_sub.list b/observation_sim/instruments/data/filters/u_sub.list index f3a2b72..3ec5a58 100644 --- a/observation_sim/instruments/data/filters/u_sub.list +++ b/observation_sim/instruments/data/filters/u_sub.list @@ -1,4 +1,4 @@ -3000 +2000 3277 3380 3485 @@ -6,4 +6,4 @@ 3703 3813 3918 -4499 \ No newline at end of file +11000 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/y_sub.list b/observation_sim/instruments/data/filters/y_sub.list index 81dc243..69a05ee 100644 --- a/observation_sim/instruments/data/filters/y_sub.list +++ b/observation_sim/instruments/data/filters/y_sub.list @@ -1,4 +1,4 @@ -9000 +2000 9322 9405 9489 @@ -6,4 +6,4 @@ 9695 9832 10024 -10590 \ No newline at end of file +11000 \ No newline at end of file diff --git a/observation_sim/instruments/data/filters/z_sub.list b/observation_sim/instruments/data/filters/z_sub.list index b0191a5..f105711 100644 --- a/observation_sim/instruments/data/filters/z_sub.list +++ b/observation_sim/instruments/data/filters/z_sub.list @@ -1,4 +1,4 @@ -7800 +2000 8494 8638 8790 @@ -6,4 +6,4 @@ 9141 9363 9663 -10590 \ No newline at end of file +11000 \ No newline at end of file -- GitLab From dfb71c23e1d3a99f280f5393b675daee4903d3ba Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Tue, 30 Jul 2024 20:06:33 +0800 Subject: [PATCH 5/5] magnitude ('mag_use_normal') calculations for galaxies and QSOs now is only based on their SEDs --- catalog/C9_Catalog.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/catalog/C9_Catalog.py b/catalog/C9_Catalog.py index 590dda6..146205d 100644 --- a/catalog/C9_Catalog.py +++ b/catalog/C9_Catalog.py @@ -557,8 +557,22 @@ class Catalog(CatalogBase): # erg/s/cm2/A --> photon/s/m2/A all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - if obj.type == 'quasar': - # integrate to get the magnitudes + # if obj.type == 'quasar': + # # integrate to get the magnitudes + # sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T + # sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array( + # sed_photon[:, 1]), interpolant='nearest') + # sed_photon = galsim.SED( + # sed_photon, wave_type='A', flux_type='1', fast=False) + # interFlux = integrate_sed_bandpass( + # sed=sed_photon, bandpass=self.filt.bandpass_full) + # obj.param['mag_use_normal'] = getABMAG( + # interFlux, self.filt.bandpass_full) + # # mag = getABMAG(interFlux, self.filt.bandpass_full) + # # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) + + # integrate to get the magnitudes + if obj.type == 'quasar' or obj.type == 'galaxy': sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array( sed_photon[:, 1]), interpolant='nearest') -- GitLab