From 641978da983f798fca5abedf60db2bfe86fee462 Mon Sep 17 00:00:00 2001 From: Jordan Mirocha Date: Tue, 27 Nov 2018 11:32:35 -0500 Subject: [PATCH] renamed litdata instances so year reflects publication, not submission. --- README.rst | 6 + ares/populations/GalaxyEnsemble.py | 30 +- ares/util/SetDefaultParameterValues.py | 1 + doc/example_litdata.rst | 16 +- doc/history.rst | 5 + doc/index.rst | 6 + input/litdata/mirocha2016.py | 240 ------------- input/litdata/mirocha2017.py | 474 +++++++++++-------------- input/litdata/mirocha2018.py | 376 +++++++++++++++----- input/litdata/mirocha2019.py | 86 +++++ 10 files changed, 629 insertions(+), 611 deletions(-) delete mode 100755 input/litdata/mirocha2016.py mode change 100644 => 100755 input/litdata/mirocha2017.py create mode 100644 input/litdata/mirocha2019.py diff --git a/README.rst b/README.rst index 82859ea74..4c9c27466 100755 --- a/README.rst +++ b/README.rst @@ -12,6 +12,12 @@ A few papers on how it works: - Parameter inference: `Mirocha, Harker, & Burns (2015) `_. - Galaxy luminosity functions: `Mirocha, Furlanetto, & Sun (2016) `_. +And some applications: + +- `Mirocha & Furlanetto (2019) `_ +- `Schneider (2018) `_ +- `Mirocha et al. (2018) `_ + Be warned: this code is still under active development -- use at your own risk! Correctness of results is not guaranteed. diff --git a/ares/populations/GalaxyEnsemble.py b/ares/populations/GalaxyEnsemble.py index 5ddefe263..9e600af06 100755 --- a/ares/populations/GalaxyEnsemble.py +++ b/ares/populations/GalaxyEnsemble.py @@ -290,8 +290,10 @@ def _gen_deterministic_histories(self): # SFR = (zform, time (but really redshift)) # So, halo identity is wrapped up in axis=0 # In Cohort, formation time defines initial mass and trajectory (in full) - histories = {'z': zall, 'w': nh, 'SFR': sfr, 'Mh': Mh, - 'MAR': mar, 'SFE': sfe, 'nh': nh} + zobs = np.array([zall] * nh.shape[0]) + histories = {'z': zall, 'zobs': zobs, 'w': nh, 'SFR': sfr, 'Mh': Mh, + 'MAR': mar, 'SFE': sfe, 'nh': nh, + 'Mg_c': Mh, 'Mg_h': Mh} # Add in formation redshifts to match shape (useful after thinning) histories['zform'] = self.tile(zall, thin) @@ -478,9 +480,9 @@ def _gen_indeterminate_histories(self): for i, zform in enumerate(all_zform): t, z = self.get_timestamps(zform) - print(zform, len(t)) + #print(zform, len(t)) - tdyn = self.halos.DynamicalTime(1e10, z) / s_per_yr + #tdyn = self.halos.DynamicalTime(1e10, z) / s_per_yr k = t.size # Grab the things that we can't modify (too much) @@ -498,6 +500,19 @@ def _gen_indeterminate_histories(self): jmax = t.size - 1 + # Sure would be nice to eliminate the following loop. + # Can only do that if there's no randomness that depends on + # previous timesteps. + + # This can go faster... + if self.pf['pop_bcycling'] == False: + # Really want an internal feedback determinism switch. + pass + + + + + Mh[i,0] = _Mh[0] Mg_h[i,0] = 0.0 Mg_c[i,0] = fbar * _Mh[0] @@ -722,8 +737,9 @@ def LuminosityFunction(self, z, x, mags=True, wave=1600., band=None): # For each unique object (or object bin), sum emission # over past episodes of star formation. - L[k] = np.trapz(L_per_msun * self.histories['SFR'][k,iz:], - dx=dt[0:-1]) + L[k] = np.sum(L_per_msun * self.histories['SFR'][k,iz:] * dt) + #L[k] = np.trapz(L_per_msun * self.histories['SFR'][k,iz:], + # dx=dt[0:-1]) else: # In this case, the time-stepping is different for each @@ -768,7 +784,7 @@ def LuminosityFunction(self, z, x, mags=True, wave=1600., band=None): if len(ages) == 0: print('ages has no elements!', k) continue - + # Need to be careful about interpolating within last # dynamical time. diff --git a/ares/util/SetDefaultParameterValues.py b/ares/util/SetDefaultParameterValues.py index 0fc8f3fda..c13302ccd 100644 --- a/ares/util/SetDefaultParameterValues.py +++ b/ares/util/SetDefaultParameterValues.py @@ -520,6 +520,7 @@ def PopulationParameters(): # Cluster-centric model "pop_poisson": False, "pop_bcycling": False, + "pop_internal_feedback": False, "pop_sf_via_inflow": True, "pop_sf_via_reservior": False, diff --git a/doc/example_litdata.rst b/doc/example_litdata.rst index 8b022f93a..adc0e5ae6 100755 --- a/doc/example_litdata.rst +++ b/doc/example_litdata.rst @@ -10,7 +10,7 @@ A very incomplete set of data from from the literature exist in ``$ARES/input/li If any of these papers ring a bell, you can check out the contents in the following way: :: - litdata = ares.util.read_lit('mirocha2016') # a self-serving example + litdata = ares.util.read_lit('mirocha2017') # a self-serving example or, look directly at the source code, which lives in ``$ARES/input/litdata``. Hopefully the contents of these files are fairly self-explanatory! @@ -57,11 +57,11 @@ Reproducing Models from *ares* Papers ------------------------------------- If you're interested in reproducing a model from a paper exactly, you can either (1) contact me directly for the model of interest, or preferably (someday) download it from my website, or (2) re-compute it yourself. In the latter case, you just need to make sure you supply the required parameters. To facilitate this, I store "parameter files" (just dictionaries) in the litdata framework as well. You can access them like any other dataset from the literature, e.g., :: - m16 = ares.util.read_lit('mirocha2016') + m17 = ares.util.read_lit('mirocha2017') A few of the models we focused on most get their own dictionary, for example our reference double power law model for the star-formation efficiency is stored in the ``dpl`` variable: :: - sim = ares.simulations.Global21cm(**m16.dpl) + sim = ares.simulations.Global21cm(**m17.dpl) sim.run() sim.GlobalSignature() # voila! @@ -69,11 +69,11 @@ Hopefully this results *exactly* in the solid black curve from Figure 2 of `Miro Alternatively, you can use the ``ParameterBundle`` framework, which also taps into our collection of data from the literature. To access the set of parameters for the "dpl" model, you simply do: :: - pars = ares.util.ParameterBundle('mirocha2016:dpl') + pars = ares.util.ParameterBundle('mirocha2017:dpl') -This tells *ares* to retrieve the ``dpl`` variable within the ``mirocha2016`` module. See :doc:`param_bundles` for more on these objects. +This tells *ares* to retrieve the ``dpl`` variable within the ``mirocha2017`` module. See :doc:`param_bundles` for more on these objects. -`Mirocha, Furlanetto, & Sun (2016) `_ (``mirocha2016``) +`Mirocha, Furlanetto, & Sun (2016) `_ (``mirocha2017``) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ This model has a few options: ``dpl``, and the extensions ``floor`` and ``steep``, as explored in the paper. @@ -95,11 +95,11 @@ To re-make the right-hand panel of Figure 1 from the paper, you could do somethi import ares - dpl = ares.util.ParameterBundle('mirocha2016:dpl') + dpl = ares.util.ParameterBundle('mirocha2017:dpl') ax = None for model in ['floor', 'dpl', 'steep']: - pars = dpl + ares.util.ParameterBundle('mirocha2016:%s' % model) + pars = dpl + ares.util.ParameterBundle('mirocha2017:%s' % model) pars.update() sim = ares.simulations.Global21cm(**pars) sim.run() diff --git a/doc/history.rst b/doc/history.rst index d69bc306a..27721475e 100755 --- a/doc/history.rst +++ b/doc/history.rst @@ -4,6 +4,11 @@ Development History Here's an attempt to keep track of major changes to the code over time, which will be tagged in the bitbucket repository with version numbers. I haven't followed conventions for version numbering so far. Instead, I've simply tagged commits with a version number when a paper is submitted using that version of the code (e.g., v0.1 and v0.2), or when a series of noteworthy improvements or bug fixes have been made (v0.3). +v0.5 +---- +- This is the version of the code used in `Mirocha & Furlanetto (2019) `_. +- Note that several ``litdata'' modules have been updated so that the year is reflective of the year the paper was *published*, not *submitted*! + v0.4 ---- - This is the version of the code used in `Mirocha et al. (2018) `_. The main addition is global Lyman-Werner feedback, which raises the minimum mass of star-forming halos self-consistently using an iterative technique. diff --git a/doc/index.rst b/doc/index.rst index 1aff1dcf3..717ca2f6a 100755 --- a/doc/index.rst +++ b/doc/index.rst @@ -17,6 +17,12 @@ A few papers on how it works: - Parameter inference: `Mirocha, Harker, & Burns (2015) `_. - Galaxy luminosity functions: `Mirocha, Furlanetto, & Sun (2017) `_. +And some applications: + +- `Mirocha & Furlanetto (2019) `_ +- `Schneider (2018) `_ +- `Mirocha et al. (2018) `_ + Be warned: this code is still under active development -- use at your own risk! Correctness of results is not guaranteed. This documentation is as much of a work in progress as the code itself, so if you encounter gaps or errors please do let me know. diff --git a/input/litdata/mirocha2016.py b/input/litdata/mirocha2016.py deleted file mode 100755 index 6a77afcd2..000000000 --- a/input/litdata/mirocha2016.py +++ /dev/null @@ -1,240 +0,0 @@ -""" -Mirocha, Furlanetto, and Sun (2017). - -Parameters defining the fiducial model (see Table 1). -""" - -from numpy import inf -from ares.physics.Constants import E_LyA - -# Calibration set! -dpl = \ -{ - - # Halos, MAR, etc. - 'pop_Tmin{0}': 1e4, - 'pop_Tmin{1}': 'pop_Tmin{0}', - 'pop_sfr_model{0}': 'sfe-func', - 'pop_sfr_model{1}': 'link:sfrd:0', - 'pop_MAR{0}': 'hmf', - 'pop_MAR_conserve_norm{0}': False, - - # Stellar pop + fesc - 'pop_sed{0}': 'eldridge2009', - 'pop_binaries{0}': False, - 'pop_Z{0}': 0.02, - 'pop_Emin{0}': 10.19, - 'pop_Emax{0}': 24.6, - 'pop_rad_yield{0}': 'from_sed', # EminNorm and EmaxNorm arbitrary now - # should make this automatic - - 'pop_fesc{0}': 0.2, - - # Solve LWB! - 'pop_solve_rte{0}': (10.2, 13.6), - - # SFE - 'pop_fstar{0}': 'pq[0]', - 'pq_func{0}[0]': 'dpl', - 'pq_func_var{0}[0]': 'Mh', - - ## - # IMPORTANT - ## - 'pq_func_par0{0}[0]': 0.05, # Table 1 in paper (last 4 rows) - 'pq_func_par1{0}[0]': 2.8e11, - 'pq_func_par2{0}[0]': 0.49, - 'pq_func_par3{0}[0]': -0.61, - 'pop_calib_L1600{0}': 1.0185e28, # Enforces Equation 13 in paper - ## - # - ## - - # Careful with X-ray heating - 'pop_sed{1}': 'mcd', - 'pop_Z{1}': 'pop_Z{0}', - 'pop_rad_yield{1}': 2.6e39, - 'pop_rad_yield_Z_index{1}': None, - 'pop_alpha{1}': -1.5, # not used unless fesc > 0 - 'pop_Emin{1}': 2e2, - 'pop_Emax{1}': 3e4, - 'pop_EminNorm{1}': 5e2, - 'pop_EmaxNorm{1}': 8e3, - 'pop_logN{1}': -inf, - - 'pop_solve_rte{1}': True, - 'tau_redshift_bins': 1000, - 'tau_approx': 'neutral', - - # Control parameters - 'include_He': True, - 'approx_He': True, - 'secondary_ionization': 3, - 'approx_Salpha': 3, - 'problem_type': 102, - 'photon_counting': True, - 'cgm_initial_temperature': 2e4, - 'cgm_recombination': 'B', - 'clumping_factor': 3., - #'smooth_derivative': 0.5, - 'final_redshift': 5, -} - -base = dpl - -_floor_specific = \ -{ - 'pq_val_floor{0}[0]': 0.005, -} - -floor = _floor_specific - -_steep_specific = \ -{ - 'pop_focc{0}': 'pq[5]', - 'pq_func{0}[5]': 'okamoto', - 'pq_func_var{0}[5]': 'Mh', - 'pq_func_par0{0}[5]': 1., - 'pq_func_par1{0}[5]': 1e9, -} - -steep = _steep_specific - -""" -Redshift-dependent options. -""" - -_flex = \ -{ - 'pq_func{0}[0]': 'dpl_arbnorm', - 'pq_func_var{0}[0]': 'Mh', - - 'pq_val_ceil{0}[0]': 1.0, - - # Standard dpl model at 10^8 Msun - 'pq_func_par0{0}[0]': 'pq[1]', - 'pq_func_par1{0}[0]': 'pq[2]', - 'pq_func_par2{0}[0]': 0.49, - 'pq_func_par3{0}[0]': -0.61, - 'pq_func_par4{0}[0]': 1e8, # Mass at which fstar,0 is defined - - # Evolving part - 'pq_func{0}[1]': 'pl', - 'pq_func_var{0}[1]': '1+z', - 'pq_func_par0{0}[1]': 0.00205, - 'pq_func_par1{0}[1]': 7., - 'pq_func_par2{0}[1]': 0., # power-law index! - - 'pq_func{0}[2]': 'pl', - 'pq_func_var{0}[2]': '1+z', - 'pq_func_par0{0}[2]': 2.8e11, - 'pq_func_par1{0}[2]': 7., - 'pq_func_par2{0}[2]': 0., # power-law index! - -} - -flex = _flex - -_flex2 = \ -{ - 'pq_func{0}[0]': 'dpl_arbnorm', - 'pq_func_var{0}[0]': 'Mh', - - 'pq_val_ceil{0}[0]': 1.0, - - # Standard dpl model at 10^8 Msun - 'pq_func_par0{0}[0]': 'pq[1]', - 'pq_func_par1{0}[0]': 'pq[2]', - 'pq_func_par2{0}[0]': 'pq[3]', - 'pq_func_par3{0}[0]': 'pq[4]', - 'pq_func_par4{0}[0]': 1e10, # Mass at which fstar,0 is defined - - # Evolving part - 'pq_func{0}[1]': 'pl', - 'pq_func_var{0}[1]': '1+z', - 'pq_func_par0{0}[1]': 0.019, # DPL model at Mh=1e10 - 'pq_func_par1{0}[1]': 7., - 'pq_func_par2{0}[1]': 0., # power-law index! - - 'pq_func{0}[2]': 'pl', - 'pq_func_var{0}[2]': '1+z', - 'pq_func_par0{0}[2]': 2.8e11, - 'pq_func_par1{0}[2]': 7., - 'pq_func_par2{0}[2]': 0., # power-law index! - - 'pq_func{0}[3]': 'pl', - 'pq_func_var{0}[3]': '1+z', - 'pq_func_par0{0}[3]': 0.49, - 'pq_func_par1{0}[3]': 7., - 'pq_func_par2{0}[3]': 0., # power-law index! - - 'pq_func{0}[4]': 'pl', - 'pq_func_var{0}[4]': '1+z', - 'pq_func_par0{0}[4]': -0.61, - 'pq_func_par1{0}[4]': 7., - 'pq_func_par2{0}[4]': 0., # power-law index! - - # Possibility of LF steepening. - 'pq_val_floor{0}[0]': 'pq[5]', - 'pq_func{0}[5]': 'pl', - 'pq_func_var{0}[5]': '1+z', - 'pq_func_par0{0}[5]': 0.0, # unused by default - 'pq_func_par1{0}[5]': 7., - 'pq_func_par2{0}[5]': 0., - - # Possibility of LF turn-over - 'pop_focc{0}': 'pq[6]', - 'pq_func{0}[6]': 'okamoto', - 'pq_func_var{0}[6]': 'Mh', - 'pq_func_par0{0}[6]': 'pq[7]', - 'pq_func_par1{0}[6]': 'pq[8]', - - 'pq_func{0}[7]': 'pl', - 'pq_func_var{0}[7]': '1+z', - 'pq_func_par0{0}[7]': 1., - 'pq_func_par1{0}[7]': 5., # effectively not in use - 'pq_func_par2{0}[7]': 0., # power-law index! - - 'pq_func{0}[8]': 'pl', - 'pq_func_var{0}[8]': '1+z', - 'pq_func_par0{0}[8]': 0., # Renders focc = 1 for all halos - 'pq_func_par1{0}[8]': 7., - 'pq_func_par2{0}[8]': 0., # power-law index! -} - -dflex = _flex2 - -fobsc = \ -{ - 'pop_fobsc{0}': 'pq[10]', - 'pop_fobsc_by_num{0}': False, # fraction of UV luminosity that gets out - 'pq_val_ceil{0}[10]': 1.0, - 'pq_val_floor{0}[10]': 0.0, - 'pq_func{0}[10]': 'log_tanh_abs', - 'pq_func_var{0}[10]': 'Mh', - 'pq_func_par0{0}[10]': 0.0, # minimal obscuration - 'pq_func_par1{0}[10]': 'pq[11]', # peak obscuration - 'pq_func_par2{0}[10]': 'pq[12]', # log transition mass - 'pq_func_par3{0}[10]': 'pq[13]', # dlogM - - 'pq_func{0}[11]': 'pl', - 'pq_func_var{0}[11]': '1+z', - 'pq_func_par0{0}[11]': 0.5, - 'pq_func_par1{0}[11]': 7., # effectively not in use - 'pq_func_par2{0}[11]': 0., # power-law index! - - 'pq_func{0}[12]': 'pl', - 'pq_func_var{0}[12]': '1+z', - 'pq_func_par0{0}[12]': 11., - 'pq_func_par1{0}[12]': 7., - 'pq_func_par2{0}[12]': 0., # power-law index! - - 'pq_func{0}[13]': 'pl', - 'pq_func_var{0}[13]': '1+z', - 'pq_func_par0{0}[13]': 1.0, - 'pq_func_par1{0}[13]': 7., - 'pq_func_par2{0}[13]': 0., # power-law index! -} - - - diff --git a/input/litdata/mirocha2017.py b/input/litdata/mirocha2017.py old mode 100644 new mode 100755 index a47f9aa55..6a77afcd2 --- a/input/litdata/mirocha2017.py +++ b/input/litdata/mirocha2017.py @@ -1,300 +1,240 @@ -import numpy as np -from mirocha2016 import dpl, flex - -_popII_models = {} -for model in ['fall', 'strong_fall', 'weak_fall', 'soft_fall']: - _popII_models[model] = flex.copy() - _popII_models[model]['pq_func_par2{0}[1]'] = 1 - -for model in ['rise', 'strong_rise', 'weak_rise', 'soft_rise']: - _popII_models[model] = flex.copy() - _popII_models[model]['pq_func_par2{0}[1]'] = -1 +""" +Mirocha, Furlanetto, and Sun (2017). -for model in ['dpl', 'strong', 'weak', 'strong_weak']: - _popII_models[model] = {} +Parameters defining the fiducial model (see Table 1). +""" -_popII_models['soft'] = {} -_popII_models['early'] = flex.copy() -_popII_models['strong_early'] = flex.copy() +from numpy import inf +from ares.physics.Constants import E_LyA -_sed_soft = \ +# Calibration set! +dpl = \ { - # Emits X-rays - "pop_lya_src{2}": False, - "pop_ion_src_cgm{2}": False, - "pop_ion_src_igm{2}": True, - "pop_heat_src_cgm{2}": False, - "pop_heat_src_igm{2}": True, - - "pop_sed{2}": 'pl', - "pop_alpha{2}": -1.5, - 'pop_logN{2}': -np.inf, - - "pop_Emin{2}": 2e2, - "pop_Emax{2}": 3e4, - "pop_EminNorm{2}": 5e2, - "pop_EmaxNorm{2}": 8e3, - - #"pop_Ex{2}": 500., - "pop_rad_yield{2}": 2.6e40, - "pop_rad_yield_units{2}": 'erg/s/SFR', -} - -_sed_soft['pop_sfr_model{2}'] = 'link:sfrd:0' -_sed_soft['pop_alpha{2}'] = -2.5 -_sed_soft['pop_solve_rte{2}'] = True -_sed_soft['tau_approx'] = 'neutral' - -_popII_updates = \ -{ - 'dpl': {}, - 'fall': {}, - 'rise': {}, - 'strong': {'pop_Z{0}': 1e-3, 'pop_rad_yield_Z_index{1}': -0.6}, - 'strong_rise': {'pop_Z{0}': 1e-3, 'pop_rad_yield_Z_index{1}': -0.6}, - 'strong_fall': {'pop_Z{0}': 1e-3, 'pop_rad_yield_Z_index{1}': -0.6}, - 'weak': {'pop_logN{1}': 22.5}, - 'weak_rise': {'pop_logN{1}': 22.5}, - 'weak_fall': {'pop_logN{1}': 22.5}, - 'strong_weak': {'pop_Z{0}': 1e-3, 'pop_rad_yield_Z_index{1}': -0.6, 'pop_logN{1}': 22.5}, - 'soft': _sed_soft, - 'soft_rise': _sed_soft, - 'soft_fall': _sed_soft, - 'early': {'pop_Tmin{0}': 500.}, - 'strong_early': {'pop_Tmin{0}': 500., 'pop_Z{0}': 1e-3, 'pop_rad_yield_Z_index{1}': -0.6}, + # Halos, MAR, etc. + 'pop_Tmin{0}': 1e4, + 'pop_Tmin{1}': 'pop_Tmin{0}', + 'pop_sfr_model{0}': 'sfe-func', + 'pop_sfr_model{1}': 'link:sfrd:0', + 'pop_MAR{0}': 'hmf', + 'pop_MAR_conserve_norm{0}': False, + + # Stellar pop + fesc + 'pop_sed{0}': 'eldridge2009', + 'pop_binaries{0}': False, + 'pop_Z{0}': 0.02, + 'pop_Emin{0}': 10.19, + 'pop_Emax{0}': 24.6, + 'pop_rad_yield{0}': 'from_sed', # EminNorm and EmaxNorm arbitrary now + # should make this automatic + + 'pop_fesc{0}': 0.2, + + # Solve LWB! + 'pop_solve_rte{0}': (10.2, 13.6), + + # SFE + 'pop_fstar{0}': 'pq[0]', + 'pq_func{0}[0]': 'dpl', + 'pq_func_var{0}[0]': 'Mh', + + ## + # IMPORTANT + ## + 'pq_func_par0{0}[0]': 0.05, # Table 1 in paper (last 4 rows) + 'pq_func_par1{0}[0]': 2.8e11, + 'pq_func_par2{0}[0]': 0.49, + 'pq_func_par3{0}[0]': -0.61, + 'pop_calib_L1600{0}': 1.0185e28, # Enforces Equation 13 in paper + ## + # + ## + + # Careful with X-ray heating + 'pop_sed{1}': 'mcd', + 'pop_Z{1}': 'pop_Z{0}', + 'pop_rad_yield{1}': 2.6e39, + 'pop_rad_yield_Z_index{1}': None, + 'pop_alpha{1}': -1.5, # not used unless fesc > 0 + 'pop_Emin{1}': 2e2, + 'pop_Emax{1}': 3e4, + 'pop_EminNorm{1}': 5e2, + 'pop_EmaxNorm{1}': 8e3, + 'pop_logN{1}': -inf, + + 'pop_solve_rte{1}': True, + 'tau_redshift_bins': 1000, + 'tau_approx': 'neutral', + + # Control parameters + 'include_He': True, + 'approx_He': True, + 'secondary_ionization': 3, + 'approx_Salpha': 3, + 'problem_type': 102, + 'photon_counting': True, + 'cgm_initial_temperature': 2e4, + 'cgm_recombination': 'B', + 'clumping_factor': 3., + #'smooth_derivative': 0.5, + 'final_redshift': 5, } -for _update in _popII_updates.keys(): - _popII_models[_update].update(_popII_updates[_update]) +base = dpl -popII_pars = _popII_models - -popII_markers = \ +_floor_specific = \ { -'dpl': 's', 'strong': '^', 'weak': 'v', -'fall': '<', 'rise': '>', 'strong_weak': 'D', -'strong_fall': '^', 'weak_rise': 'v', -'strong_rise': '^', 'weak_fall': 'v', 'soft': 's', -'soft_fall': '<', 'soft_rise': '>', -'early': '<', 'strong_early': '<', + 'pq_val_floor{0}[0]': 0.005, } -# Lotta trouble just to get 'dpl' first in the list... -_amp = ['weak', 'strong'] -_timing = ['rise', 'fall'] -_all = _amp + _timing -popII_models = ['dpl'] + _all + ['strong_weak'] + ['early', 'strong_early'] - -for e1 in _amp: - for e2 in _timing: - popII_models.append('{0!s}_{1!s}'.format(e1, e2)) +floor = _floor_specific -popII_models.extend(['soft', 'soft_fall', 'soft_rise']) - -# relative to mirocha2016:dpl -_generic_updates = \ +_steep_specific = \ { - 'initial_redshift': 60, - 'pop_zform{0}': 60, - 'pop_zform{1}': 60, - 'feedback_LW': True, - 'sam_dz': 0.1, - 'kill_redshift': 5.6, - - 'feedback_LW_Mmin_rtol': 0, - 'feedback_LW_sfrd_rtol': 5e-2, - 'feedback_LW_sfrd_popid': 2, - 'feedback_LW_maxiter': 50, - 'feedback_LW_mixup_freq': 5, - 'feedback_LW_mixup_delay': 10, - 'pop_time_limit_delay{2}': True, + 'pop_focc{0}': 'pq[5]', + 'pq_func{0}[5]': 'okamoto', + 'pq_func_var{0}[5]': 'Mh', + 'pq_func_par0{0}[5]': 1., + 'pq_func_par1{0}[5]': 1e9, } -# This can lead to pickling issues...argh -#halos = ares.physics.HaloMassFunction() -#barrier_A = lambda zz: halos.VirialMass(1e4, zz) +steep = _steep_specific """ -First model: constant SFR in minihalos. +Redshift-dependent options. """ -_low = \ -{ - - 'pop_zform{2}': 60, - 'pop_zform{3}': 60, - 'pop_sfr_model{2}': 'sfr-func', - 'pop_sfr{2}': 1e-5, - 'pop_sfr_cross_threshold{2}': False, - 'pop_sed{2}': 'eldridge2009', - 'pop_binaries{2}': False, - 'pop_Z{2}': 1e-3, - 'pop_Emin{2}': 10.19, - 'pop_Emax{2}': 24.6, +_flex = \ +{ + 'pq_func{0}[0]': 'dpl_arbnorm', + 'pq_func_var{0}[0]': 'Mh', + + 'pq_val_ceil{0}[0]': 1.0, + + # Standard dpl model at 10^8 Msun + 'pq_func_par0{0}[0]': 'pq[1]', + 'pq_func_par1{0}[0]': 'pq[2]', + 'pq_func_par2{0}[0]': 0.49, + 'pq_func_par3{0}[0]': -0.61, + 'pq_func_par4{0}[0]': 1e8, # Mass at which fstar,0 is defined + + # Evolving part + 'pq_func{0}[1]': 'pl', + 'pq_func_var{0}[1]': '1+z', + 'pq_func_par0{0}[1]': 0.00205, + 'pq_func_par1{0}[1]': 7., + 'pq_func_par2{0}[1]': 0., # power-law index! + + 'pq_func{0}[2]': 'pl', + 'pq_func_var{0}[2]': '1+z', + 'pq_func_par0{0}[2]': 2.8e11, + 'pq_func_par1{0}[2]': 7., + 'pq_func_par2{0}[2]': 0., # power-law index! - 'pop_heat_src_igm{2}': False, - 'pop_ion_src_igm{2}': False, +} - # Solve LWB! - 'pop_solve_rte{2}': (10.2, 13.6), +flex = _flex - 'pop_sed{2}': 'eldridge2009', - 'pop_rad_yield{2}': 'from_sed', - 'pop_Z{2}': 0.02, +_flex2 = \ +{ + 'pq_func{0}[0]': 'dpl_arbnorm', + 'pq_func_var{0}[0]': 'Mh', + + 'pq_val_ceil{0}[0]': 1.0, + + # Standard dpl model at 10^8 Msun + 'pq_func_par0{0}[0]': 'pq[1]', + 'pq_func_par1{0}[0]': 'pq[2]', + 'pq_func_par2{0}[0]': 'pq[3]', + 'pq_func_par3{0}[0]': 'pq[4]', + 'pq_func_par4{0}[0]': 1e10, # Mass at which fstar,0 is defined + + # Evolving part + 'pq_func{0}[1]': 'pl', + 'pq_func_var{0}[1]': '1+z', + 'pq_func_par0{0}[1]': 0.019, # DPL model at Mh=1e10 + 'pq_func_par1{0}[1]': 7., + 'pq_func_par2{0}[1]': 0., # power-law index! + + 'pq_func{0}[2]': 'pl', + 'pq_func_var{0}[2]': '1+z', + 'pq_func_par0{0}[2]': 2.8e11, + 'pq_func_par1{0}[2]': 7., + 'pq_func_par2{0}[2]': 0., # power-law index! + + 'pq_func{0}[3]': 'pl', + 'pq_func_var{0}[3]': '1+z', + 'pq_func_par0{0}[3]': 0.49, + 'pq_func_par1{0}[3]': 7., + 'pq_func_par2{0}[3]': 0., # power-law index! - # Radiative knobs - 'pop_fesc_LW{2}': 1., - 'pop_fesc{2}': 0.0, - 'pop_rad_yield{3}': 2.6e39, + 'pq_func{0}[4]': 'pl', + 'pq_func_var{0}[4]': '1+z', + 'pq_func_par0{0}[4]': -0.61, + 'pq_func_par1{0}[4]': 7., + 'pq_func_par2{0}[4]': 0., # power-law index! - # Other stuff needed for X-rays - 'pop_sfr_model{3}': 'link:sfrd:2', - 'pop_sed{3}': 'mcd', - 'pop_Z{3}': 1e-3, - 'pop_rad_yield_Z_index{3}': None, - 'pop_Emin{3}': 2e2, - 'pop_Emax{3}': 3e4, - 'pop_EminNorm{3}': 5e2, - 'pop_EmaxNorm{3}': 8e3, - 'pop_ion_src_cgm{3}': False, - 'pop_logN{3}': -np.inf, - - 'pop_solve_rte{3}': True, - - #'pop_Mmin{0}': 'link:Mmax_active:2', - - # Tmin here just an initial guess -- will get modified by feedback. - 'pop_Tmin{2}': 500., - 'pop_Tmin{0}': None, - 'pop_Mmin{0}': 'link:Mmax:2', - 'pop_Tmax_ceil{2}': 1e6, - 'pop_sfr_cross_threshold{2}': False, - - 'pop_time_limit{2}': 2.5, - 'pop_bind_limit{2}': 1e51, - 'pop_abun_limit{2}': None, + # Possibility of LF steepening. + 'pq_val_floor{0}[0]': 'pq[5]', + 'pq_func{0}[5]': 'pl', + 'pq_func_var{0}[5]': '1+z', + 'pq_func_par0{0}[5]': 0.0, # unused by default + 'pq_func_par1{0}[5]': 7., + 'pq_func_par2{0}[5]': 0., - # Acknowledging that our mean metallicity is kludgey - # Note that this renders the stellar mass meaningless (it'll be zero). - 'pop_mass_yield{2}': 1.0, - 'pop_metal_yield{2}': 1.0, + # Possibility of LF turn-over + 'pop_focc{0}': 'pq[6]', + 'pq_func{0}[6]': 'okamoto', + 'pq_func_var{0}[6]': 'Mh', + 'pq_func_par0{0}[6]': 'pq[7]', + 'pq_func_par1{0}[6]': 'pq[8]', + + 'pq_func{0}[7]': 'pl', + 'pq_func_var{0}[7]': '1+z', + 'pq_func_par0{0}[7]': 1., + 'pq_func_par1{0}[7]': 5., # effectively not in use + 'pq_func_par2{0}[7]': 0., # power-law index! + + 'pq_func{0}[8]': 'pl', + 'pq_func_var{0}[8]': '1+z', + 'pq_func_par0{0}[8]': 0., # Renders focc = 1 for all halos + 'pq_func_par1{0}[8]': 7., + 'pq_func_par2{0}[8]': 0., # power-law index! } -low = _generic_updates.copy() -low.update(_low) - -high = low.copy() -high['pop_sed{2}'] = 'schaerer2002' -high['pop_mass{2}'] = 120. - -med = low.copy() -med['pop_sed{2}'] = 'schaerer2002' -med['pop_mass{2}'] = 5. - -LWoff = low.copy() -LWoff['pop_fesc_LW{2}'] = 0 - -bb = low.copy() -bb['pop_sed{2}'] = 'bb' -bb['pop_rad_yield{2}'] = 1e5 -bb['pop_rad_yield_units{2}'] = 'photons/baryon' -bb['pop_EminNorm{2}'] = 11.2 -bb['pop_EmaxNorm{2}'] = 13.6 -bb['pop_Emin{2}'] = 1. -bb['pop_Emax{2}'] = 1e2 -bb['pop_temperature{2}'] = 1e5 - -popII_like = low.copy() -popII_like['pop_sed{2}'] = 'eldridge2009' -popII_like['pop_Z{2}'] = 1e-3 - -""" -Second model: constant SFE in minihalos. -""" -#_csfe_specific = \ -#{ -# 'pop_sfr_model{2}': 'sfe-func', -# 'pop_sfr{2}': None, -# 'pop_fstar{2}': 1e-4, -#} -# -#csfe = dpl.copy() -#csfe.update(_generic_updates) -#csfe.update(_csfr_specific) -#csfe.update(_csfe_specific) -# -#_csff_specific = \ -#{ -# 'pop_sfr{2}': 'pq[2]', -# 'pq_func{2}[2]': 'pl', -# 'pq_func_var{2}[2]': 'Mh', -# 'pq_func_par0{2}[2]': 1e-4, -# 'pq_func_par1{2}[2]': 1e8, -# 'pq_func_par2{2}[2]': 1., -#} -# -#csff = csfr.copy() -#csff.update(_csff_specific) +dflex = _flex2 -""" -Third model: extrapolated SFE in minihalos (i.e., same SFE as atomic halos). -""" -#_xsfe_specific = \ -#{ -# -# 'pop_fesc_LW{0}': 'pq[101]', -# 'pq_func{0}[101]': 'astep', -# 'pq_func_var{0}[101]': 'Mh', -# 'pq_func_par0{0}[101]': 1., -# 'pq_func_par1{0}[101]': 1., -# 'pq_func_par2{0}[101]': (barrier_A, 'z', 1), -# -# 'pop_fesc{0}': 'pq[102]', -# 'pq_func{0}[102]': 'astep', -# 'pq_func_var{0}[102]': 'Mh', -# 'pq_func_par0{0}[102]': 0., # No LyC from minihalos by default -# 'pq_func_par1{0}[102]': 0.1, -# 'pq_func_par2{0}[102]': (barrier_A, 'z', 1), -# -# 'pop_Tmin{0}': 500., -# 'pop_Mmin{1}': 'pop_Mmin{0}', -# -# 'pop_sfr_model{1}': 'link:sfe:0', -# -# # X-ray sources -# 'pop_rad_yield{1}': 'pq[103]', -# 'pq_func{1}[103]': 'astep', -# 'pq_func_var{1}[103]': 'Mh', -# 'pq_func_par0{1}[103]': 2.6e39, -# 'pq_func_par1{1}[103]': 2.6e39, -# 'pq_func_par2{1}[103]': (barrier_A, 'z', 1), -#} -# -#xsfe = dpl.copy() -#xsfe.update(_generic_updates) -#xsfe.update(_xsfe_specific) - -csfr_blobs = \ +fobsc = \ { - 'blob_names': ['popII_sfrd_tot', 'popIII_sfrd_tot', - 'popII_sfrd_bc', 'popIII_sfrd_bc', - 'popII_Mmin', 'popIII_Mmin', - 'popII_Mmax', 'popIII_Mmax', - 'popII_nh', 'popIII_nh'], - 'blob_ivars': ('z', np.arange(5, 60.1, 0.1)), - 'blob_funcs': ['pops[0].SFRD', 'pops[2].SFRD', 'pops[0].SFRD_at_threshold', - 'pops[2].SFRD_at_threshold', 'pops[0].Mmin', 'pops[2].Mmin', - 'pops[0].Mmax', 'pops[2].Mmax', 'pops[0].nactive', 'pops[2].nactive'], + 'pop_fobsc{0}': 'pq[10]', + 'pop_fobsc_by_num{0}': False, # fraction of UV luminosity that gets out + 'pq_val_ceil{0}[10]': 1.0, + 'pq_val_floor{0}[10]': 0.0, + 'pq_func{0}[10]': 'log_tanh_abs', + 'pq_func_var{0}[10]': 'Mh', + 'pq_func_par0{0}[10]': 0.0, # minimal obscuration + 'pq_func_par1{0}[10]': 'pq[11]', # peak obscuration + 'pq_func_par2{0}[10]': 'pq[12]', # log transition mass + 'pq_func_par3{0}[10]': 'pq[13]', # dlogM + + 'pq_func{0}[11]': 'pl', + 'pq_func_var{0}[11]': '1+z', + 'pq_func_par0{0}[11]': 0.5, + 'pq_func_par1{0}[11]': 7., # effectively not in use + 'pq_func_par2{0}[11]': 0., # power-law index! + + 'pq_func{0}[12]': 'pl', + 'pq_func_var{0}[12]': '1+z', + 'pq_func_par0{0}[12]': 11., + 'pq_func_par1{0}[12]': 7., + 'pq_func_par2{0}[12]': 0., # power-law index! + + 'pq_func{0}[13]': 'pl', + 'pq_func_var{0}[13]': '1+z', + 'pq_func_par0{0}[13]': 1.0, + 'pq_func_par1{0}[13]': 7., + 'pq_func_par2{0}[13]': 0., # power-law index! } -csfe_blobs = csfr_blobs -csff_blobs = csfr_blobs -dpl_blobs = \ -{ - 'blob_names': ['popII_sfrd_tot', 'popII_Mmin', 'popII_Mmax'], - 'blob_ivars': ('z', np.arange(5, 60.1, 0.1)), - 'blob_funcs': ['pops[0].SFRD', 'pops[0].Mmin', 'pops[0].Mmax'], -} diff --git a/input/litdata/mirocha2018.py b/input/litdata/mirocha2018.py index 55e90b62b..318d9b91d 100644 --- a/input/litdata/mirocha2018.py +++ b/input/litdata/mirocha2018.py @@ -1,102 +1,300 @@ -from mirocha2016 import dpl, dflex -from ares.util import ParameterBundle as PB -from ares.physics.Constants import nu_0_mhz, h_p, erg_per_ev - -## -# All need this! -## -base = PB(verbose=0, **dpl) \ - + PB(verbose=0, **dflex) \ - + PB(verbose=0, 'dust:var_beta') - -lf = \ +import numpy as np +from mirocha2017 import dpl, flex + +_popII_models = {} +for model in ['fall', 'strong_fall', 'weak_fall', 'soft_fall']: + _popII_models[model] = flex.copy() + _popII_models[model]['pq_func_par2{0}[1]'] = 1 + +for model in ['rise', 'strong_rise', 'weak_rise', 'soft_rise']: + _popII_models[model] = flex.copy() + _popII_models[model]['pq_func_par2{0}[1]'] = -1 + +for model in ['dpl', 'strong', 'weak', 'strong_weak']: + _popII_models[model] = {} + +_popII_models['soft'] = {} +_popII_models['early'] = flex.copy() +_popII_models['strong_early'] = flex.copy() + +_sed_soft = \ { - 'pq_func_par0{0}[12]': 11.6690040263, - 'pq_func_par0{0}[13]': 0.855092040361, - 'pq_func_par2{0}[0]': 1.03890453721, - 'pq_func_par2{0}[13]': -0.292169721106, - 'pq_func_par3{0}[0]': -0.372920485712, - 'pq_func_par2{0}[12]': -0.742521570574, - 'pq_func_par2{0}[11]': 0.0561650504878, - 'pq_func_par1{0}[0]': 1.8979710707e+11, - 'pq_func_par0{0}[11]': 0.249869796896, - 'pq_func_par0{0}[0]': 0.171411643892, + + # Emits X-rays + "pop_lya_src{2}": False, + "pop_ion_src_cgm{2}": False, + "pop_ion_src_igm{2}": True, + "pop_heat_src_cgm{2}": False, + "pop_heat_src_igm{2}": True, + + "pop_sed{2}": 'pl', + "pop_alpha{2}": -1.5, + 'pop_logN{2}': -np.inf, + + "pop_Emin{2}": 2e2, + "pop_Emax{2}": 3e4, + "pop_EminNorm{2}": 5e2, + "pop_EmaxNorm{2}": 8e3, + + #"pop_Ex{2}": 500., + "pop_rad_yield{2}": 2.6e40, + "pop_rad_yield_units{2}": 'erg/s/SFR', } -cold = \ +_sed_soft['pop_sfr_model{2}'] = 'link:sfrd:0' +_sed_soft['pop_alpha{2}'] = -2.5 +_sed_soft['pop_solve_rte{2}'] = True +_sed_soft['tau_approx'] = 'neutral' + +_popII_updates = \ { - 'approx_thermal_history': 'exp', - 'load_ics': 'parametric', - 'inits_Tk_p0': 194.002300947, - 'inits_Tk_p1': 1.20941098917, - 'inits_Tk_p2': -6.0088645858, - - 'pq_func_par2{0}[4]': 0.0661410442476, - 'pq_func_par0{0}[8]': 19701832.4424, - 'pq_func_par2{0}[5]': -0.624613431068, - 'pq_func_par2{0}[8]': -0.74856754179, - 'pq_func_par2{0}[7]': 1.1001036864, - 'inits_Tk_p1': 1.20941098917, - 'pq_func_par0{0}[7]': 1.27092398097, - 'inits_Tk_p2': -6.0088645858, - 'pq_func_par0{0}[4]': -0.459851109813, - 'pq_func_par2{0}[1]': -0.978610145341, - 'pop_rad_yield{1}': 1.46150981378e+40, - 'pq_func_par0{0}[5]': 0.0338201803874, - 'pq_func_par2{0}[2]': 0.451825099378, - 'inits_Tk_p0': 194.002300947, - 'pq_func_par0{0}[2]': 2.57313940031e+11, - 'pq_func_par2{0}[3]': 0.0305460186899, - 'pq_func_par0{0}[3]': 0.899806521888, - 'pop_Tmin{0}': 18875.1354673, - 'pq_func_par0{0}[1]': 0.0195652006263, + 'dpl': {}, + 'fall': {}, + 'rise': {}, + 'strong': {'pop_Z{0}': 1e-3, 'pop_rad_yield_Z_index{1}': -0.6}, + 'strong_rise': {'pop_Z{0}': 1e-3, 'pop_rad_yield_Z_index{1}': -0.6}, + 'strong_fall': {'pop_Z{0}': 1e-3, 'pop_rad_yield_Z_index{1}': -0.6}, + 'weak': {'pop_logN{1}': 22.5}, + 'weak_rise': {'pop_logN{1}': 22.5}, + 'weak_fall': {'pop_logN{1}': 22.5}, + 'strong_weak': {'pop_Z{0}': 1e-3, 'pop_rad_yield_Z_index{1}': -0.6, 'pop_logN{1}': 22.5}, + 'soft': _sed_soft, + 'soft_rise': _sed_soft, + 'soft_fall': _sed_soft, + 'early': {'pop_Tmin{0}': 500.}, + 'strong_early': {'pop_Tmin{0}': 500., 'pop_Z{0}': 1e-3, 'pop_rad_yield_Z_index{1}': -0.6}, +} + +for _update in _popII_updates.keys(): + _popII_models[_update].update(_popII_updates[_update]) +popII_pars = _popII_models + +popII_markers = \ +{ +'dpl': 's', 'strong': '^', 'weak': 'v', +'fall': '<', 'rise': '>', 'strong_weak': 'D', +'strong_fall': '^', 'weak_rise': 'v', +'strong_rise': '^', 'weak_fall': 'v', 'soft': 's', +'soft_fall': '<', 'soft_rise': '>', +'early': '<', 'strong_early': '<', } -E21 = nu_0_mhz * 1e6 * h_p / erg_per_ev -nu_min = 1e7 # in Hz -nu_max = 1e12 -Emin = nu_min * (h_p / erg_per_ev) # 1 GHz -> eV -Emax = nu_max * (h_p / erg_per_ev) # 3 GHz -> eV +# Lotta trouble just to get 'dpl' first in the list... +_amp = ['weak', 'strong'] +_timing = ['rise', 'fall'] +_all = _amp + _timing +popII_models = ['dpl'] + _all + ['strong_weak'] + ['early', 'strong_early'] + +for e1 in _amp: + for e2 in _timing: + popII_models.append('{0!s}_{1!s}'.format(e1, e2)) -radio = \ +popII_models.extend(['soft', 'soft_fall', 'soft_rise']) + +# relative to mirocha2016:dpl +_generic_updates = \ { - 'pop_sfr_model{2}': 'link:sfrd:0', - 'pop_sed{2}': 'pl', - 'pop_alpha{2}': -0.7, - 'pop_Emin{2}': Emin, - 'pop_Emax{2}': Emax, - 'pop_EminNorm{2}': None, - 'pop_EmaxNorm{2}': None, - 'pop_Enorm{2}': E21, # 1.4 GHz - 'pop_rad_yield_units{2}': 'erg/s/sfr/hz', + 'initial_redshift': 60, + 'pop_zform{0}': 60, + 'pop_zform{1}': 60, + 'feedback_LW': True, + 'sam_dz': 0.1, + 'kill_redshift': 5.6, - # Solution method - 'pop_solve_rte{2}': True, - 'pop_radio_src{2}': True, - 'pop_lw_src{2}': False, - 'pop_lya_src{2}': False, + 'feedback_LW_Mmin_rtol': 0, + 'feedback_LW_sfrd_rtol': 5e-2, + 'feedback_LW_sfrd_popid': 2, + 'feedback_LW_maxiter': 50, + 'feedback_LW_mixup_freq': 5, + 'feedback_LW_mixup_delay': 10, + 'pop_time_limit_delay{2}': True, +} + +# This can lead to pickling issues...argh +#halos = ares.physics.HaloMassFunction() +#barrier_A = lambda zz: halos.VirialMass(1e4, zz) + +""" +First model: constant SFR in minihalos. +""" +_low = \ +{ + + 'pop_zform{2}': 60, + 'pop_zform{3}': 60, + + 'pop_sfr_model{2}': 'sfr-func', + 'pop_sfr{2}': 1e-5, + 'pop_sfr_cross_threshold{2}': False, + 'pop_sed{2}': 'eldridge2009', + 'pop_binaries{2}': False, + 'pop_Z{2}': 1e-3, + 'pop_Emin{2}': 10.19, + 'pop_Emax{2}': 24.6, + 'pop_heat_src_igm{2}': False, 'pop_ion_src_igm{2}': False, - 'pop_ion_src_cgm{2}': False, - - # Best fitting parameters - 'pq_func_par2{0}[4]': -0.367507400514, - 'pq_func_par0{0}[8]': 88971.2720411, - 'pq_func_par2{0}[5]': 0.432165531607, - 'pq_func_par2{0}[8]': -0.459908108301, - 'pq_func_par2{0}[7]': 0.0804573261674, - 'pq_func_par0{0}[7]': 0.755153610936, - 'pop_Tmin{0}': 19345.5297538, - 'pq_func_par0{0}[4]': -0.332641612164, - 'pq_func_par2{0}[1]': -0.280720519413, - 'pop_rad_yield{1}': 4.20070368093e+40, - 'pq_func_par0{0}[5]': 0.0200865171317, - 'pq_func_par2{0}[2]': 0.100749514086, - 'pq_func_par0{0}[2]': 2.19664981719e+11, - 'pq_func_par2{0}[3]': -0.251920672275, - 'pq_func_par0{0}[3]': 0.772229741005, - 'pop_zdead{2}': 17.0252061155, - 'pop_rad_yield{2}': 1.37568148182e+32, - 'pq_func_par0{0}[1]': 0.0268026393973, + + # Solve LWB! + 'pop_solve_rte{2}': (10.2, 13.6), + + 'pop_sed{2}': 'eldridge2009', + 'pop_rad_yield{2}': 'from_sed', + 'pop_Z{2}': 0.02, + + # Radiative knobs + 'pop_fesc_LW{2}': 1., + 'pop_fesc{2}': 0.0, + 'pop_rad_yield{3}': 2.6e39, + + # Other stuff needed for X-rays + 'pop_sfr_model{3}': 'link:sfrd:2', + 'pop_sed{3}': 'mcd', + 'pop_Z{3}': 1e-3, + 'pop_rad_yield_Z_index{3}': None, + 'pop_Emin{3}': 2e2, + 'pop_Emax{3}': 3e4, + 'pop_EminNorm{3}': 5e2, + 'pop_EmaxNorm{3}': 8e3, + 'pop_ion_src_cgm{3}': False, + 'pop_logN{3}': -np.inf, + + 'pop_solve_rte{3}': True, + + #'pop_Mmin{0}': 'link:Mmax_active:2', + + # Tmin here just an initial guess -- will get modified by feedback. + 'pop_Tmin{2}': 500., + 'pop_Tmin{0}': None, + 'pop_Mmin{0}': 'link:Mmax:2', + 'pop_Tmax_ceil{2}': 1e6, + 'pop_sfr_cross_threshold{2}': False, + + 'pop_time_limit{2}': 2.5, + 'pop_bind_limit{2}': 1e51, + 'pop_abun_limit{2}': None, + + # Acknowledging that our mean metallicity is kludgey + # Note that this renders the stellar mass meaningless (it'll be zero). + 'pop_mass_yield{2}': 1.0, + 'pop_metal_yield{2}': 1.0, +} + +low = _generic_updates.copy() +low.update(_low) + +high = low.copy() +high['pop_sed{2}'] = 'schaerer2002' +high['pop_mass{2}'] = 120. + +med = low.copy() +med['pop_sed{2}'] = 'schaerer2002' +med['pop_mass{2}'] = 5. + +LWoff = low.copy() +LWoff['pop_fesc_LW{2}'] = 0 + +bb = low.copy() +bb['pop_sed{2}'] = 'bb' +bb['pop_rad_yield{2}'] = 1e5 +bb['pop_rad_yield_units{2}'] = 'photons/baryon' +bb['pop_EminNorm{2}'] = 11.2 +bb['pop_EmaxNorm{2}'] = 13.6 +bb['pop_Emin{2}'] = 1. +bb['pop_Emax{2}'] = 1e2 +bb['pop_temperature{2}'] = 1e5 + +popII_like = low.copy() +popII_like['pop_sed{2}'] = 'eldridge2009' +popII_like['pop_Z{2}'] = 1e-3 + +""" +Second model: constant SFE in minihalos. +""" +#_csfe_specific = \ +#{ +# 'pop_sfr_model{2}': 'sfe-func', +# 'pop_sfr{2}': None, +# 'pop_fstar{2}': 1e-4, +#} +# +#csfe = dpl.copy() +#csfe.update(_generic_updates) +#csfe.update(_csfr_specific) +#csfe.update(_csfe_specific) +# +#_csff_specific = \ +#{ +# 'pop_sfr{2}': 'pq[2]', +# 'pq_func{2}[2]': 'pl', +# 'pq_func_var{2}[2]': 'Mh', +# 'pq_func_par0{2}[2]': 1e-4, +# 'pq_func_par1{2}[2]': 1e8, +# 'pq_func_par2{2}[2]': 1., +#} +# +#csff = csfr.copy() +#csff.update(_csff_specific) + +""" +Third model: extrapolated SFE in minihalos (i.e., same SFE as atomic halos). +""" +#_xsfe_specific = \ +#{ +# +# 'pop_fesc_LW{0}': 'pq[101]', +# 'pq_func{0}[101]': 'astep', +# 'pq_func_var{0}[101]': 'Mh', +# 'pq_func_par0{0}[101]': 1., +# 'pq_func_par1{0}[101]': 1., +# 'pq_func_par2{0}[101]': (barrier_A, 'z', 1), +# +# 'pop_fesc{0}': 'pq[102]', +# 'pq_func{0}[102]': 'astep', +# 'pq_func_var{0}[102]': 'Mh', +# 'pq_func_par0{0}[102]': 0., # No LyC from minihalos by default +# 'pq_func_par1{0}[102]': 0.1, +# 'pq_func_par2{0}[102]': (barrier_A, 'z', 1), +# +# 'pop_Tmin{0}': 500., +# 'pop_Mmin{1}': 'pop_Mmin{0}', +# +# 'pop_sfr_model{1}': 'link:sfe:0', +# +# # X-ray sources +# 'pop_rad_yield{1}': 'pq[103]', +# 'pq_func{1}[103]': 'astep', +# 'pq_func_var{1}[103]': 'Mh', +# 'pq_func_par0{1}[103]': 2.6e39, +# 'pq_func_par1{1}[103]': 2.6e39, +# 'pq_func_par2{1}[103]': (barrier_A, 'z', 1), +#} +# +#xsfe = dpl.copy() +#xsfe.update(_generic_updates) +#xsfe.update(_xsfe_specific) + +csfr_blobs = \ +{ + 'blob_names': ['popII_sfrd_tot', 'popIII_sfrd_tot', + 'popII_sfrd_bc', 'popIII_sfrd_bc', + 'popII_Mmin', 'popIII_Mmin', + 'popII_Mmax', 'popIII_Mmax', + 'popII_nh', 'popIII_nh'], + 'blob_ivars': ('z', np.arange(5, 60.1, 0.1)), + 'blob_funcs': ['pops[0].SFRD', 'pops[2].SFRD', 'pops[0].SFRD_at_threshold', + 'pops[2].SFRD_at_threshold', 'pops[0].Mmin', 'pops[2].Mmin', + 'pops[0].Mmax', 'pops[2].Mmax', 'pops[0].nactive', 'pops[2].nactive'], +} + +csfe_blobs = csfr_blobs +csff_blobs = csfr_blobs + +dpl_blobs = \ +{ + 'blob_names': ['popII_sfrd_tot', 'popII_Mmin', 'popII_Mmax'], + 'blob_ivars': ('z', np.arange(5, 60.1, 0.1)), + 'blob_funcs': ['pops[0].SFRD', 'pops[0].Mmin', 'pops[0].Mmax'], } + diff --git a/input/litdata/mirocha2019.py b/input/litdata/mirocha2019.py new file mode 100644 index 000000000..6d420474a --- /dev/null +++ b/input/litdata/mirocha2019.py @@ -0,0 +1,86 @@ +from mirocha2017 import dpl, dflex +from ares.util import ParameterBundle as PB +from ares.physics.Constants import nu_0_mhz, h_p, erg_per_ev + +## +# All need this! +## +base = PB(verbose=0, **dpl) \ + + PB(verbose=0, **dflex) \ + + PB('dust:var_beta', verbose=0) + +cold = \ +{ + # New base_kwargs + 'approx_thermal_history': 'exp', + 'load_ics': 'parametric', + + # Copy-pasted over from best fits. + 'pq_func_par0{0}[1]': 0.018949141521, + 'pq_func_par0{0}[2]': 2.31016897023e+11, + 'pq_func_par0{0}[3]': 0.924929473207, + 'pq_func_par0{0}[4]': -0.345183026665, + 'pq_func_par2{0}[1]': -1.93710531526, + 'pq_func_par2{0}[2]': -0.0401589348891, + 'pq_func_par2{0}[3]': 0.496357034538, + 'pq_func_par2{0}[4]': -0.85185812704, + 'pq_func_par0{0}[5]': 0.0180819517762, + 'pq_func_par2{0}[5]': 0.633935667655, + 'pq_func_par0{0}[7]': 2.47340783415, + 'pq_func_par0{0}[8]': 22676116.9726, + 'pq_func_par2{0}[7]': 1.29978775053, + 'pq_func_par2{0}[8]': 2.51794223721, + 'inits_Tk_p0': 203.477407296, + 'inits_Tk_p1': 1.23612113669, + 'inits_Tk_p2': -6.89882925178, + 'pop_Tmin{0}': 18729.8083367, + 'pop_rad_yield{1}': 6.75648524393e+40, +} + +E21 = nu_0_mhz * 1e6 * h_p / erg_per_ev +nu_min = 1e7 # in Hz +nu_max = 1e12 +Emin = nu_min * (h_p / erg_per_ev) # 1 GHz -> eV +Emax = nu_max * (h_p / erg_per_ev) # 3 GHz -> eV + +radio = \ +{ + # New base_kwargs + 'pop_sfr_model{2}': 'link:sfrd:0', + 'pop_sed{2}': 'pl', + 'pop_alpha{2}': -0.7, + 'pop_Emin{2}': Emin, + 'pop_Emax{2}': Emax, + 'pop_EminNorm{2}': None, + 'pop_EmaxNorm{2}': None, + 'pop_Enorm{2}': E21, # 1.4 GHz + 'pop_rad_yield_units{2}': 'erg/s/sfr/hz', + + 'pop_solve_rte{2}': True, + 'pop_radio_src{2}': True, + 'pop_lw_src{2}': False, + 'pop_lya_src{2}': False, + 'pop_heat_src_igm{2}': False, + 'pop_ion_src_igm{2}': False, + 'pop_ion_src_cgm{2}': False, + + # Best fitting parameters + 'pq_func_par0{0}[1]': 0.0210546853809, + 'pq_func_par0{0}[2]': 3.3655857398e+11, + 'pq_func_par0{0}[3]': 0.83260420295, + 'pq_func_par0{0}[4]': -0.38893617798, + 'pq_func_par2{0}[1]': -1.65784835208, + 'pq_func_par2{0}[2]': 1.10105346147, + 'pq_func_par2{0}[3]': 0.213603033125, + 'pq_func_par2{0}[4]': -0.484706407852, + 'pq_func_par0{0}[5]': 0.0187421303099, + 'pq_func_par2{0}[5]': 0.390393567924, + 'pq_func_par0{0}[7]': 1.86134140738, + 'pq_func_par0{0}[8]': 632300.756957, + 'pq_func_par2{0}[7]': 0.404071491728, + 'pq_func_par2{0}[8]': -0.545718558631, + 'pop_rad_yield{2}': 1.23318228932e+33, + 'pop_zdead{2}': 16.3318857443, + 'pop_Tmin{0}': 23809.941444, + 'pop_rad_yield{1}': 9.74775993145e+41, +}