Merge branch 'milky_way_extinction' into 'develop'

Add Milky Way extinction and extend the range for SED integration

See merge request !27

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.zeros(len(ra_arr))
 
         for igals in range(ngals):
             # # (TEST)
@@ -544,14 +551,28 @@ 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
         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')

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"

observation_sim/instruments/data/filters/fgs_sub.list

-3000
+2000
 4500
 4750
 5000

observation_sim/instruments/data/filters/g_sub.list

-3800
+2000
 4217
 4432
 4631
@@ -6,4 +6,4 @@
 5002
 5179
 5354
-5799
+11000
\ No newline at end of file

observation_sim/instruments/data/filters/i_sub.list

-6600
+2000
 7061
 7255
 7448
@@ -6,4 +6,4 @@
 7833
 8027
 8226
-8999
+11000
\ No newline at end of file

observation_sim/instruments/data/filters/nuv_sub.list

-2513
+2000
 2621
 2716
 2805
@@ -6,4 +6,4 @@
 2969
 3050
 3132
-3499
+11000
\ No newline at end of file

observation_sim/instruments/data/filters/r_sub.list

-5100
+2000
 5642
 5821
 6001
@@ -6,4 +6,4 @@
 6363
 6547
 6735
-7199
+11000
\ No newline at end of file

observation_sim/instruments/data/filters/u_sub.list

-3000
+2000
 3277
 3380
 3485
@@ -6,4 +6,4 @@
 3703
 3813
 3918
-4499
+11000
\ No newline at end of file

observation_sim/instruments/data/filters/y_sub.list

-9000
+2000
 9322
 9405
 9489
@@ -6,4 +6,4 @@
 9695
 9832
 10024
-10590
+11000
\ No newline at end of file

observation_sim/instruments/data/filters/z_sub.list

-7800
+2000
 8494
 8638
 8790
@@ -6,4 +6,4 @@
 9141
 9363
 9663
-10590
+11000
\ No newline at end of file

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, 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
@@ -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, nest=True)
         return self.ebv_planck[pix] * self.Rv
 
     def apply_extinction(self, spec, Av=1.0):
@@ -90,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

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.)