Improved IGM absorption

Author: cschreib

CHANGELOG

@@ -1,3 +1,10 @@
+v1.1.0
+======
+
+New/updated feature:
+ - Improved the treatment of IGM absorption. Now the absorption is computed on a per-galaxy basis, based on the redshift. Before that, the IGM absorption was quite poorly implemented, as a single value of the absorption for the entire catalog with no redshift dependence. The old behavior of the code can be recovered by setting the "naive_igm" flag on the command line.
+
+
 v1.0.9
 ======
 

CMakeLists.txt

@@ -12,7 +12,7 @@ find_package(phypp)
 set(EGG_SHARE_DIR ${CMAKE_INSTALL_PREFIX}/share/egg)
 install(FILES
     share/ir_lib_ce01.fits share/ir_lib_cs17.fits
-    share/mass_func_candels.fits share/opt_lib_fast.fits
+    share/mass_func_candels.fits share/opt_lib_fast.fits share/opt_lib_fast_noigm.fits
     DESTINATION ${EGG_SHARE_DIR} COMPONENT data)
 install(DIRECTORY share/filter-db
     DESTINATION ${EGG_SHARE_DIR} COMPONENT data)

share/opt_lib_fast_noigm.fits

@@ -60,6 +60,11 @@ int phypp_main(int argc, char* argv[]) {
     bool no_dust = false;        // do not generate fluxes from dust origin
     bool no_random = false;      // disable all randomization of parameters
                                  // (just passive, M*, z, and position + pos angle)
+
+    // Approximations
+    bool naive_igm = false;      // use a fixed IGM absorption for all galaxies (faster)
+
+
     // Debug - testing
     bool no_passive_lir = false; // do not generate dust for passive galaxies
     bool magdis_tdust = false;   // use Magdis+12 Tdust evolution (lower)
@@ -75,7 +80,7 @@ int phypp_main(int argc, char* argv[]) {
 
     // SED libraries
     std::string ir_lib_file = egg_share_dir+"ir_lib_cs17.fits";
-    std::string opt_lib_file = egg_share_dir+"opt_lib_fast.fits";
+    std::string opt_lib_file;
 
     // Filter library
     std::string filter_db_file = filters_dir+"db.dat";
@@ -127,7 +132,7 @@ int phypp_main(int argc, char* argv[]) {
         verbose, name(tseed, "seed"), name(tcosmo, "cosmo"),
         name(input_cat_file, "input_cat"), selection_band, bands, rfbands, help, list_bands,
         clust_r0, clust_r1, clust_lambda, clust_eta, clust_fclust_mlim, clust_fclust_lom,
-        clust_fclust_him, clust_urnd_mlim, magdis_tdust
+        clust_fclust_him, clust_urnd_mlim, magdis_tdust, naive_igm
     ));
 
     if (help) {
@@ -249,6 +254,14 @@ int phypp_main(int argc, char* argv[]) {
         }
     }
 
+    if (opt_lib_file.empty()) {
+        if (naive_igm) {
+            opt_lib_file = egg_share_dir+"opt_lib_fast.fits";
+        } else {
+            opt_lib_file = egg_share_dir+"opt_lib_fast_noigm.fits";
+        }
+    }
+
     // Initialize random seed and cosmology
     // ------------------------------------
     auto seed = make_seed(tseed);
@@ -1582,6 +1595,47 @@ if (!no_flux) {
         bulge, disk
     };
 
+    auto apply_igm = [&](double z, const vec1f& lam, vec1f& sed) {
+        // http://adsabs.harvard.edu/abs/1995ApJ...441...18M
+        // TODO: check this implementation someday, I suspect this is wrong or
+        // very approximate (taken directly from FAST)
+
+        double da; {
+            double l0 = 1050.0*(1.0 + z);
+            double l1 = 1170.0*(1.0 + z);
+            uint_t nstep = 100;
+            vec1d tl = rgen(l0, l1, nstep);
+            vec1d ptau = exp(-3.6e-3*pow(tl/1216.0, 3.46));
+            da = total(ptau)*(l1-l0)/nstep/(120.0*(1.0 + z));
+        }
+
+        double db; {
+            double l0 = 920.0*(1.0 + z);
+            double l1 = 1015.0*(1.0 + z);
+            uint_t nstep = 100;
+            vec1d tl = rgen(l0, l1, nstep);
+            vec1d ptau = exp(-1.7e-3*pow(tl/1026.0, 3.46) - 1.2e-3*pow(tl/972.5, 3.46) -
+                9.3e-4*pow(tl/950.0, 3.46));
+            db = total(ptau)*(l1-l0)/nstep/(95.0*(1.0 + z));
+        }
+
+        uint_t l0 = lower_bound(lam, 0.0921);
+        uint_t l1 = lower_bound(lam, 0.1026);
+        uint_t l2 = lower_bound(lam, 0.1216);
+
+        for (uint_t l : range(l0)) {
+            sed.safe[l] = 0;
+        }
+        for (uint_t l : range(l0, l1)) {
+            sed.safe[l] *= db;
+        }
+        for (uint_t l : range(l1, l2)) {
+            sed.safe[l] *= da;
+        }
+
+        return lam;
+    };
+
     auto get_flux = [&](const vec1f& m, const vec1u& optsed, vec2f& flux, vec2f& rfmag,
         bool no_ir, component cmp) {
 
@@ -1606,6 +1660,11 @@ if (!no_flux) {
                 merge_add(orlam, irlam, orsed, irsed, rlam, rsed);
             }
 
+            // Apply IGM absorption
+            if (!naive_igm && !no_stellar) {
+                apply_igm(out.z[i], rlam, rsed);
+            }
+
             // Obtain rest-frame magnitudes
             if (!rffilters.empty()) {
                 // Redshift the SED @ 10pc
@@ -1755,7 +1814,7 @@ void print_help() {
         }
     };
 
-    print("egg-gencat v1.0.9");
+    print("egg-gencat v1.1.0");
     print("usage: egg-gencat [options]\n");
 
     print("List of generic options:");
@@ -1797,6 +1856,11 @@ void print_help() {
         "fully deterministic recipes");
     print("");
 
+    print("List of approximations:");
+    argdoc("naive_igm", "[flag]", "use a fixed IGM absorption for all galaxies instead of using "
+        "a redshift-dependent prescription.");
+    print("");
+
     print("List of sky position related options:");
     argdoc("ra0", "[double, degrees]", "sky position of the center of the field "
         "(right ascension, default: 53.558750)");