Commit b9723717 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

Add Milky Way extinction and extend the range for SED integration

parent 6e13d206
...@@ -111,10 +111,14 @@ class Catalog(CatalogBase): ...@@ -111,10 +111,14 @@ class Catalog(CatalogBase):
self.max_size = 0. self.max_size = 0.
# [TODO] Milky Way extinction # [TODO] Milky Way extinction
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 = ExtinctionMW()
self.mw_ext.init_ext_model(model_name="odonnell") self.mw_ext.init_ext_model(model_name="odonnell")
self.mw_ext.load_Planck_ext( self.mw_ext.load_Planck_ext(
file_path="/public/home/fangyuedong/project/ext_maps/planck/HFI_CompMap_ThermalDustModel_2048_R1.20.fits") 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"]: 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 # Get the cloest star catalog file
...@@ -255,7 +259,10 @@ class Catalog(CatalogBase): ...@@ -255,7 +259,10 @@ class Catalog(CatalogBase):
) )
# [TODO] get Milky Way extinction AVs # [TODO] get Milky Way extinction AVs
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) MW_Av_arr = self.mw_ext.Av_from_Planck(ra=ra_arr, dec=dec_arr)
else:
MW_Av_arr = np.zeros(len(ra_arr))
for igals in range(ngals): for igals in range(ngals):
# # (TEST) # # (TEST)
...@@ -544,14 +551,28 @@ class Catalog(CatalogBase): ...@@ -544,14 +551,28 @@ class Catalog(CatalogBase):
y = speci(lamb) y = speci(lamb)
# [TODO] Apply Milky Way extinction # [TODO] Apply Milky Way extinction
if obj.type != 'star': if obj.type != 'star' and ("enable_mw_ext_gal" in self.config["catalog_options"] and self.config["catalog_options"]["enable_mw_ext_gal"]):
self.mw_ext.apply_extinction(y, Av=obj.mw_Av) y = self.mw_ext.apply_extinction(y, Av=obj.mw_Av)
# erg/s/cm2/A --> photon/s/m2/A # erg/s/cm2/A --> photon/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
if obj.type == 'quasar': # 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 # integrate to get the magnitudes
if obj.type == 'quasar' or obj.type == 'galaxy':
sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array( sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(
sed_photon[:, 1]), interpolant='nearest') sed_photon[:, 1]), interpolant='nearest')
......
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
# can add some of the command-line arguments here as well; # can add some of the command-line arguments here as well;
# ok to pass either way or both, as long as they are consistent # ok to pass either way or both, as long as they are consistent
work_dir: "/public/home/fangyuedong/project/workplace/" 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 and run counter are used to name the outputs
project_cycle: 9 project_cycle: 9
...@@ -44,11 +44,15 @@ catalog_options: ...@@ -44,11 +44,15 @@ catalog_options:
star_only: NO star_only: NO
# Only simulate galaxies? # Only simulate galaxies?
galaxy_only: YES galaxy_only: NO
# rotate galaxy ellipticity # rotate galaxy ellipticity
rotateEll: 0. # [degree] 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 # Observation setting
############################################### ###############################################
...@@ -68,7 +72,7 @@ obs_setting: ...@@ -68,7 +72,7 @@ obs_setting:
run_pointings: [0, 1, 2, 3, 4] run_pointings: [0, 1, 2, 3, 4]
# Whether to enable astrometric modeling # Whether to enable astrometric modeling
enable_astrometric_model: True enable_astrometric_model: YES
# Cut by saturation magnitude in which band? # Cut by saturation magnitude in which band?
cut_in_band: "z" cut_in_band: "z"
......
3800 2000
4217 4217
4432 4432
4631 4631
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
5002 5002
5179 5179
5354 5354
5799 11000
\ No newline at end of file \ No newline at end of file
6600 2000
7061 7061
7255 7255
7448 7448
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
7833 7833
8027 8027
8226 8226
8999 11000
\ No newline at end of file \ No newline at end of file
2513 2000
2621 2621
2716 2716
2805 2805
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
2969 2969
3050 3050
3132 3132
3499 11000
\ No newline at end of file \ No newline at end of file
5100 2000
5642 5642
5821 5821
6001 6001
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
6363 6363
6547 6547
6735 6735
7199 11000
\ No newline at end of file \ No newline at end of file
3000 2000
3277 3277
3380 3380
3485 3485
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
3703 3703
3813 3813
3918 3918
4499 11000
\ No newline at end of file \ No newline at end of file
9000 2000
9322 9322
9405 9405
9489 9489
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
9695 9695
9832 9832
10024 10024
10590 11000
\ No newline at end of file \ No newline at end of file
7800 2000
8494 8494
8638 8638
8790 8790
...@@ -6,4 +6,4 @@ ...@@ -6,4 +6,4 @@
9141 9141
9363 9363
9663 9663
10590 11000
\ No newline at end of file \ No newline at end of file
...@@ -13,7 +13,7 @@ class ExtinctionMW(object): ...@@ -13,7 +13,7 @@ class ExtinctionMW(object):
@staticmethod @staticmethod
def radec2pix(ra, dec, NSIDE=2048): 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, nest=True, lonlat=True)
def init_ext_model(self, model_name="odonnell", Av=1.0, Rv=3.1, lamb=None): def init_ext_model(self, model_name="odonnell", Av=1.0, Rv=3.1, lamb=None):
self.model_name = model_name self.model_name = model_name
...@@ -79,7 +79,8 @@ class ExtinctionMW(object): ...@@ -79,7 +79,8 @@ class ExtinctionMW(object):
u.degree, frame='icrs').galactic u.degree, frame='icrs').galactic
l, b = c.l.radian, c.b.radian l, b = c.l.radian, c.b.radian
NSIDE = hp.pixelfunc.get_nside(self.ebv_planck) 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, nest=True)
return self.ebv_planck[pix] * self.Rv return self.ebv_planck[pix] * self.Rv
def apply_extinction(self, spec, Av=1.0): def apply_extinction(self, spec, Av=1.0):
...@@ -90,6 +91,6 @@ class ExtinctionMW(object): ...@@ -90,6 +91,6 @@ class ExtinctionMW(object):
raise ValueError( raise ValueError(
"Need to initialize the extinction model (init_ext_model) first") "Need to initialize the extinction model (init_ext_model) first")
scale = 10**(-.4 * self.ext * Av) scale = 10**(-.4 * self.ext * Av)
print("scale = ", scale) # print("scale = ", scale)
spec *= scale spec *= scale
return spec return spec
...@@ -27,7 +27,7 @@ def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param): ...@@ -27,7 +27,7 @@ def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param):
InputDark=None) InputDark=None)
else: else:
chip.img, _ = chip_utils.add_poisson(img=chip.img, chip.img, _ = chip_utils.add_poisson(img=chip.img,
chip=self, chip=chip,
exptime=exptime, exptime=exptime,
poisson_noise=chip.poisson_noise, poisson_noise=chip.poisson_noise,
dark_noise=0.) dark_noise=0.)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment