Radio excess from ACG, MCG and PBH

build_cffi.py

@@ -1,19 +1,28 @@
 """Build the C code with CFFI."""
-import os
+import os, platform
 from cffi import FFI
 
 ffi = FFI()
 LOCATION = os.path.dirname(os.path.abspath(__file__))
 CLOC = os.path.join(LOCATION, "src", "py21cmfast", "src")
 include_dirs = [CLOC]
 
-# =================================================================
-# Set compilation arguments dependent on environment... a bit buggy
-# =================================================================
+# ========================================================================
+# Set compilation arguments dependent on environment and OS... a bit buggy
+# ========================================================================
+OS_Name = platform.system()
+
 if "DEBUG" in os.environ:
-    extra_compile_args = ["-fopenmp", "-w", "-g", "-O0"]
+    if OS_Name == 'Darwin':
+        # fopenmp is known to fail on mac
+        extra_compile_args = ["-w", "-g", "-O0"]
+    else:
+        extra_compile_args = ["-fopenmp", "-w", "-g", "-O0"]
 else:
-    extra_compile_args = ["-fopenmp", "-Ofast", "-w"]
+    if OS_Name == 'Darwin':
+        extra_compile_args = ["-Ofast", "-w"]
+    else:
+        extra_compile_args = ["-fopenmp", "-Ofast", "-w"]
 
 # Set the C-code logging level.
 # If DEBUG is set, we default to the highest level, but if not,

src/py21cmfast/_utils.py

@@ -401,11 +401,15 @@ def __init__(self, *args, **kwargs):
                 setattr(self, "_" + k, v)
 
         if kwargs:
+            pass
+            '''
+            # JSC: cannot silence this using warnings module
             logger.warning(
                 "The following parameters to {thisclass} are not supported: {lst}".format(
                     thisclass=self.__class__.__name__, lst=list(kwargs.keys())
                 )
             )
+            '''
 
     def convert(self, key, val):
         """Make any conversions of values before saving to the instance."""

src/py21cmfast/inputs.py

@@ -484,13 +484,16 @@ class UserParams(StructWithDefaults):
     def USE_INTERPOLATION_TABLES(self):
         """Whether to use interpolation tables for integrals, speeding things up."""
         if self._USE_INTERPOLATION_TABLES is None:
+            '''
+            JSC: force-silence warnings
             warnings.warn(
                 "The USE_INTERPOLATION_TABLES setting has changed in v3.1.2 to be "
                 "default True. You can likely ignore this warning, but if you relied on"
                 "having USE_INTERPOLATION_TABLES=False by *default*, please set it "
                 "explicitly. To silence this warning, set it explicitly to True. This"
                 "warning will be removed in v4."
             )
+            '''
             self._USE_INTERPOLATION_TABLES = True
 
         return self._USE_INTERPOLATION_TABLES
@@ -518,6 +521,7 @@ def POWER_SPECTRUM(self):
         See :func:`power_spectrum_model` for a string representation.
         """
         if self.USE_RELATIVE_VELOCITIES:
+            ''' jsc: this is too noisy
             if self._POWER_SPECTRUM != 5 or (
                 isinstance(self._POWER_SPECTRUM, str)
                 and self._POWER_SPECTRUM.upper() != "CLASS"
@@ -526,6 +530,7 @@ def POWER_SPECTRUM(self):
                     "Automatically setting POWER_SPECTRUM to 5 (CLASS) as you are using "
                     "relative velocities"
                 )
+            '''
             return 5
         else:
             if isinstance(self._POWER_SPECTRUM, str):
@@ -626,6 +631,12 @@ class FlagOptions(StructWithDefaults):
         Determines whether to use a fixed vcb=VAVG (*regardless* of USE_RELATIVE_VELOCITIES). It includes the average effect of velocities but not its fluctuations. See Muñoz+21 (2110.13919).
     USE_VELS_AUX: bool, optional
         Auxiliary variable (not input) to check if minihaloes are being used without relative velocities and complain
+    USE_RADIO_ACG: bool, optional
+        Determines whether to use radio excess background from ACG, if True then AstroParams.fR is used
+    USE_RADIO_MCG: bool, optional
+        Determines whether to use radio excess background from MCG, if True then AstroParams.fR_mini is used
+    Calibrate_EoR_feedback: bool
+        Whether to calibrate EoR/photo-heating feedback for Pop II/III radio excess, used only when MINIHALO is activated
     """
 
     _ffi = ffi
@@ -640,6 +651,9 @@ class FlagOptions(StructWithDefaults):
         "M_MIN_in_Mass": False,
         "PHOTON_CONS": False,
         "FIX_VCB_AVG": False,
+        "USE_RADIO_ACG": False,
+        "USE_RADIO_MCG": False,
+        "Calibrate_EoR_feedback": False,
     }
 
     @property
@@ -796,6 +810,17 @@ class AstroParams(StructWithDefaults):
         Impact of the LW feedback on Mturn for minihaloes. Default is 22.8685 and 0.47 following Machacek+01, respectively. Latest simulations suggest 2.0 and 0.6. See Sec 2 of Muñoz+21 (2110.13919).
     A_VCB, BETA_VCB: float, optional
         Impact of the DM-baryon relative velocities on Mturn for minihaloes. Default is 1.0 and 1.8, and agrees between different sims. See Sec 2 of Muñoz+21 (2110.13919).
+    fR: float, optional
+        Radio efficiency for galaxies, normalised to 1 for modern day galaxies. Given in log10 units.
+    aR: float, optional
+        Power-law energy spectra index
+    fR_mini: float, optional
+        Radio efficiency for molecularly cooling galaxies (MCG) in mini-halos, normalised to 1 for modern day galaxies. Given in log10 units.
+    aR_mini: float, optional
+        Power-law energy spectra index for mini-halos
+        PBH X-ray emission efficiency
+    Radio_Zmin: float, optional
+        Turn off radio emmisivity below this redshift, a phenomenological param motivated by ARCADE2 upper limit
     """
 
     _ffi = ffi
@@ -824,6 +849,12 @@ class AstroParams(StructWithDefaults):
         "BETA_LW": 0.6,
         "A_VCB": 1.0,
         "BETA_VCB": 1.8,
+        # Radio excess and PBH params
+        "fR": -10.0,
+        "aR": 0.7,
+        "fR_mini": -10.0,
+        "aR_mini": 0.7,
+        "Radio_Zmin": 0.0,
     }
 
     def __init__(
@@ -846,6 +877,8 @@ def convert(self, key, val):
             "L_X",
             "L_X_MINI",
             "X_RAY_Tvir_MIN",
+            "fR",
+            "fR_mini",
         ]:
             return 10**val
         else:
@@ -953,3 +986,11 @@ def validate_all_inputs(
         logger.warning(
             "USE_MINI_HALOS needs USE_RELATIVE_VELOCITIES to get the right evolution!"
         )
+
+    if flag_options is not None:
+        # Some requirements for excess radio background
+        if flag_options.USE_RADIO_MCG and not flag_options.USE_MINI_HALOS:
+            raise ValueError("USE_RADIO_MCG requires USE_MINI_HALOS!")
+        if flag_options.Calibrate_EoR_feedback and not flag_options.USE_RADIO_MCG:
+            raise ValueError("Calibrate_EoR_feedback requires USE_RADIO_MCG!")
+

src/py21cmfast/outputs.py

@@ -410,6 +410,8 @@ def _get_box_structures(self) -> dict[str, dict | tuple[int]]:
             "x_e_box": shape,
             "Tk_box": shape,
             "J_21_LW_box": shape,
+            "Trad_box": shape,
+            "History_box": shape,
         }
 
     @cached_property
@@ -422,6 +424,14 @@ def global_Ts(self):
         else:
             return np.mean(self.Ts_box)
 
+    @cached_property
+    def global_Trad(self):
+        """Global Radio Temp."""
+        if "Trad_box" not in self._computed_arrays:
+            raise AttributeError("Something went wrong, maybe Trad not computed?")
+        else:
+            return np.mean(self.Trad_box)
+
     @cached_property
     def global_Tk(self):
         """Global (mean) Tk."""
@@ -457,6 +467,8 @@ def get_required_input_arrays(self, input_box: _BaseOutputStruct) -> list[str]:
             required += [
                 "Tk_box",
                 "x_e_box",
+                "Trad_box",
+                "History_box",
             ]
             if self.flag_options.USE_MINI_HALOS:
                 required += ["J_21_LW_box"]
@@ -567,7 +579,7 @@ def get_required_input_arrays(self, input_box: _BaseOutputStruct) -> list[str]:
         elif isinstance(input_box, PerturbedField):
             required += ["density"]
         elif isinstance(input_box, TsBox):
-            required += ["J_21_LW_box", "x_e_box", "Tk_box"]
+            required += ["J_21_LW_box", "x_e_box", "Tk_box", "Trad_box", "History_box"]
         elif isinstance(input_box, IonizedBox):
             required += ["z_re_box", "Gamma12_box"]
             if self.flag_options.INHOMO_RECO: