Commit 1aa4f533 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

1. add option for galaxy milky way extinction in overall_config 

2. bug fix: planck map nested ordering
parent 5e8554bd

catalog/C9_Catalog.py

+14 −7
Original line number Diff line number Diff line
@@ -111,10 +111,14 @@ class Catalog(CatalogBase): 
        self.max_size = 0.



        # [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.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")

                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

        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

config/config_overall.yaml

+7 −3
Original line number Diff line number Diff line
@@ -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/mock_objects/ExtinctionMW.py

+3 −2
Original line number Diff line number Diff line
@@ -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):