Add HWO residual contrast and exozodi realism

Author: Ben Mazin

examples/hwo_space.toml

@@ -29,8 +29,8 @@
 #   * Detector cap is set far above the expected rates so the first pass
 #     is not saturation limited.
 #   * Space zodiacal backgrounds are included through `background.space_zodi`:
-#     local zodi defaults to an average V-band surface brightness, and the
-#     exozodi level comes from [scene.exozodi].
+#     local zodi defaults to an average V-band surface brightness, and exozodi
+#     uses a HOSTS-style level draw plus a simple inclined-disk multiplier.
 #   * Calibration floor and residual systematic are set to negligible
 #     positive values; HWO contrast is represented by the planet contrast,
 #     not by an imposed calibration floor.
@@ -111,9 +111,9 @@ separation_mas = 200.0
 position_angle_deg = 0.0
 
 [scene.exozodi]
-# Solar-System-like exozodi level. `background.space_zodi` treats this as a
-# linear multiplier on a 1-zodi visible-light surface brightness near the
-# Earth-equivalent insolation distance.
+# Fixed fallback exozodi level. The HWO zodi stage below opts into a
+# HOSTS-style sampled level for realism; this scene value is used only if
+# that stage is returned to exozodi_level_distribution="fixed".
 level_zodi = 1.0
 
 [observation]
@@ -167,9 +167,9 @@ pointing_jitter_rms_mas = 0.5
 # space-wavefront stage provides a near-unity Strehl spectrum. A real
 # multi-sub-aperture injection stage would replace this with independent
 # modes; the broadcast keeps raw null depth degenerate but lets
-# performance.standard run. Leave apply_strehl_coupling=false because the downstream
-# injection.single_mode_fiber stage applies the Ruilier-Cassaing SMF
-# coupling exactly once.
+# performance.standard run. Leave apply_strehl_coupling=false because the
+# downstream injection.single_mode_fiber stage applies the Ruilier-Cassaing
+# SMF coupling exactly once.
 id = "injection_seed"
 type = "injection.ideal"
 n_modes = 36
@@ -288,9 +288,10 @@ band = "I"
 [[pipeline.stages]]
 # Separation-dependent planet throughput + stellar-disk leakage for the
 # 36-mode monolithic-pupil surrogate above. With no atmosphere, AO, or
-# formation-flying residual stage, the realized-null MC is disabled. A flat
-# 1e-10 residual-starlight floor anchors the first-pass HWO design contrast;
-# later realism stages can replace it with a wavefront-derived contrast map.
+# formation-flying residual stage, the realized-null MC is disabled. The
+# residual-starlight floor uses a lightweight wavelength- and separation-
+# dependent HWO I-band contrast surface as a stand-in for a future
+# wavefront-derived contrast map.
 # The grid is HWO-specific: at 0.8 µm on an 8 m primary, λ/D = 20.6 mas.
 # The sub-mas samples keep the finite-stellar-disk integral honest for a
 # 1.86 mas diameter star; the plotted science range still extends through
@@ -317,23 +318,31 @@ n_azimuths = 1
 n_realizations_per_azimuth = 0
 n_mc_realizations = 0
 mc_seed = 0
-residual_starlight_contrast = [1.0e-10]
+residual_starlight_table_path = "data/hwo_reference/residual_contrast_i_band.csv"
 residual_starlight_combine_mode = "max"
 
 [[pipeline.stages]]
 # Space zodiacal backgrounds. The default local-zodi value is an average
 # V-band surface brightness converted to I-band with an approximate solar color;
 # this HWO example pins a low/ecliptic-pole value because deep direct-imaging
-# targets should avoid the high-zodi sky. The 1-zodi exozodi surface brightness
-# is scaled from the host star's I-band magnitude and then multiplied by
-# [scene.exozodi].level_zodi above. The background is routed through the same
-# off-axis dark-port throughput that sensitivity.detection_curve computed for
-# the planet separation. The solid angle defaults to 1 (lambda/D)^2.
+# targets should avoid the high-zodi sky. Exozodi uses a HOSTS-like lognormal
+# dust-level draw and an inclined smooth-disk surface-brightness multiplier,
+# then routes through the same off-axis dark-port throughput that
+# sensitivity.detection_curve computed for the planet separation. The solid
+# angle defaults to 1 (lambda/D)^2.
 id = "zodi"
 type = "background.space_zodi"
 band = "I"
 local_zodi_v_mag_arcsec2 = 23.3
 derive_exozodi_mag_from_scene_jmag = true
+exozodi_level_distribution = "hosts2020_lognormal"
+exozodi_hosts2020_median_zodi = 3.0
+exozodi_hosts2020_sigma_dex = 0.47712125472
+exozodi_level_max_zodi = 100.0
+exozodi_morphology = "smooth_inclined_disk"
+exozodi_inclination_distribution = "isotropic"
+exozodi_max_inclination_deg = 80.0
+exozodi_max_surface_brightness_multiplier = 5.0
 routing = "scene_planet_throughput"
 
 [[pipeline.stages]]

nullsim/data/hwo_reference/residual_contrast_i_band.csv

@@ -0,0 +1,40 @@
+separation_mas,wavelength_um,contrast
+0.05,0.76,1.0e-6
+0.05,0.80,8.0e-7
+0.05,0.84,1.0e-6
+0.10,0.76,5.0e-7
+0.10,0.80,4.0e-7
+0.10,0.84,5.0e-7
+0.20,0.76,1.5e-7
+0.20,0.80,1.2e-7
+0.20,0.84,1.5e-7
+0.50,0.76,3.0e-8
+0.50,0.80,2.5e-8
+0.50,0.84,3.0e-8
+1.00,0.76,1.0e-8
+1.00,0.80,8.0e-9
+1.00,0.84,1.0e-8
+2.00,0.76,3.0e-9
+2.00,0.80,2.5e-9
+2.00,0.84,3.0e-9
+5.00,0.76,1.0e-9
+5.00,0.80,8.0e-10
+5.00,0.84,1.0e-9
+10.00,0.76,4.0e-10
+10.00,0.80,3.0e-10
+10.00,0.84,4.0e-10
+20.60,0.76,1.6e-10
+20.60,0.80,1.2e-10
+20.60,0.84,1.6e-10
+50.00,0.76,1.2e-10
+50.00,0.80,1.0e-10
+50.00,0.84,1.2e-10
+100.00,0.76,1.1e-10
+100.00,0.80,1.0e-10
+100.00,0.84,1.1e-10
+200.00,0.76,1.0e-10
+200.00,0.80,1.0e-10
+200.00,0.84,1.0e-10
+500.00,0.76,1.4e-10
+500.00,0.80,1.2e-10
+500.00,0.84,1.4e-10

nullsim/stages/background/space_zodi.py

@@ -8,7 +8,12 @@
 
 import numpy as np
 from numpy.typing import NDArray
-from pydantic import BaseModel, ConfigDict, Field as PydanticField
+from pydantic import (
+    BaseModel,
+    ConfigDict,
+    Field as PydanticField,
+    model_validator,
+)
 
 from nullsim.physics.photometry import (
     photon_rate_per_m2_per_bin,
@@ -77,6 +82,29 @@ class Params(BaseModel):
         band: Literal["I", "Y", "J", "H", "K"] | None = None
         local_zodi_level: float = PydanticField(default=1.0, ge=0.0)
         exozodi_level_zodi: float | None = PydanticField(default=None, ge=0.0)
+        exozodi_level_distribution: Literal["fixed", "hosts2020_lognormal"] = "fixed"
+        exozodi_hosts2020_median_zodi: float = PydanticField(default=3.0, ge=0.0)
+        exozodi_hosts2020_sigma_dex: float = PydanticField(
+            default=math.log10(3.0),
+            ge=0.0,
+        )
+        exozodi_level_min_zodi: float = PydanticField(default=0.0, ge=0.0)
+        exozodi_level_max_zodi: float | None = PydanticField(default=None, ge=0.0)
+        exozodi_morphology: Literal[
+            "smooth_face_on",
+            "smooth_inclined_disk",
+        ] = "smooth_face_on"
+        exozodi_inclination_deg: float = PydanticField(default=0.0, ge=0.0, lt=90.0)
+        exozodi_inclination_distribution: Literal["fixed", "isotropic"] = "fixed"
+        exozodi_max_inclination_deg: float = PydanticField(
+            default=80.0,
+            gt=0.0,
+            lt=90.0,
+        )
+        exozodi_max_surface_brightness_multiplier: float = PydanticField(
+            default=5.0,
+            ge=1.0,
+        )
         local_zodi_v_mag_arcsec2: float = PydanticField(
             default=DEFAULT_LOCAL_ZODI_V_MAG_ARCSEC2,
             gt=0.0,
@@ -95,6 +123,17 @@ class Params(BaseModel):
         n_subapertures: int | None = PydanticField(default=None, ge=1)
         area_source: Literal["effective", "collecting_area_per_mode"] = "effective"
 
+        @model_validator(mode="after")
+        def _validate_distribution_bounds(self) -> "SpaceZodiBackground.Params":
+            if (
+                self.exozodi_level_max_zodi is not None
+                and self.exozodi_level_max_zodi < self.exozodi_level_min_zodi
+            ):
+                raise ValueError(
+                    "exozodi_level_max_zodi must be >= exozodi_level_min_zodi."
+                )
+            return self
+
     def apply(self, state: SimulationState) -> SimulationState:
         params = self.params_model
         wavelength_grid_um = _require_wavelength_grid_um(state)
@@ -142,10 +181,19 @@ def apply(self, state: SimulationState) -> SimulationState:
             routing_throughput=1.0,
             n_dark_ports=routing.n_dark_ports,
         )
-        exozodi_level = _resolve_exozodi_level(state, params.exozodi_level_zodi)
+        exozodi_level, exozodi_level_source = _resolve_exozodi_level(
+            state,
+            params,
+            self.context.rng,
+        )
+        exozodi_morphology = _resolve_exozodi_morphology(params, self.context.rng)
+        exozodi_effective_level = (
+            float(exozodi_level)
+            * float(exozodi_morphology["surface_brightness_multiplier"])
+        )
         exozodi_base_rate_per_bin_hz = _surface_brightness_rate_per_bin_hz(
             mag_arcsec2=exozodi_mag,
-            level=exozodi_level,
+            level=exozodi_effective_level,
             wavelength_grid_um=wavelength_grid_um,
             bandwidth_um=bandwidth_um,
             band=band,
@@ -176,6 +224,29 @@ def apply(self, state: SimulationState) -> SimulationState:
             "band": band,
             "local_zodi_level": float(params.local_zodi_level),
             "exozodi_level_zodi": float(exozodi_level),
+            "exozodi_effective_level_zodi": float(exozodi_effective_level),
+            "exozodi_level_source": exozodi_level_source,
+            "exozodi_level_distribution": params.exozodi_level_distribution,
+            "exozodi_hosts2020_median_zodi": float(
+                params.exozodi_hosts2020_median_zodi
+            ),
+            "exozodi_hosts2020_sigma_dex": float(
+                params.exozodi_hosts2020_sigma_dex
+            ),
+            "exozodi_level_min_zodi": float(params.exozodi_level_min_zodi),
+            "exozodi_level_max_zodi": (
+                None if params.exozodi_level_max_zodi is None
+                else float(params.exozodi_level_max_zodi)
+            ),
+            "exozodi_morphology": params.exozodi_morphology,
+            "exozodi_inclination_deg": float(exozodi_morphology["inclination_deg"]),
+            "exozodi_inclination_source": exozodi_morphology["inclination_source"],
+            "exozodi_surface_brightness_multiplier": float(
+                exozodi_morphology["surface_brightness_multiplier"]
+            ),
+            "exozodi_raw_surface_brightness_multiplier": float(
+                exozodi_morphology["raw_surface_brightness_multiplier"]
+            ),
             "local_zodi_v_mag_arcsec2": float(params.local_zodi_v_mag_arcsec2),
             "exozodi_v_mag_arcsec2": float(params.exozodi_v_mag_arcsec2),
             "local_zodi_mag_arcsec2": float(local_mag),
@@ -229,7 +300,9 @@ def apply(self, state: SimulationState) -> SimulationState:
                 "dark-port routing throughput. Local zodi uses an absolute "
                 "V-band surface brightness converted with approximate "
                 "solar colors; exozodi defaults to a 1-zodi surface "
-                "brightness scaled by the host-star band magnitude."
+                "brightness scaled by the host-star band magnitude. Optional "
+                "HOSTS-style zodi-level sampling and inclined smooth-disk "
+                "surface-brightness scaling are applied before routing."
             ),
         }
         payload.update(curve_background_payload)
@@ -501,13 +574,84 @@ def _exozodi_mag_arcsec2(
 
 def _resolve_exozodi_level(
     state: SimulationState,
-    configured: float | None,
-) -> float:
-    if configured is not None:
-        return float(configured)
-    if state.scene.exozodi_level_zodi is None:
-        return 0.0
-    return float(state.scene.exozodi_level_zodi)
+    params: SpaceZodiBackground.Params,
+    rng: np.random.Generator,
+) -> tuple[float, str]:
+    """Resolve or sample the physical exozodi dust level in zodis."""
+
+    if params.exozodi_level_distribution == "fixed":
+        if params.exozodi_level_zodi is not None:
+            return float(params.exozodi_level_zodi), "params.exozodi_level_zodi"
+        if state.scene.exozodi_level_zodi is None:
+            return 0.0, "scene.exozodi_level_zodi missing"
+        return float(state.scene.exozodi_level_zodi), "scene.exozodi_level_zodi"
+
+    if params.exozodi_level_distribution != "hosts2020_lognormal":
+        raise ValueError(
+            "background.space_zodi: unsupported exozodi_level_distribution="
+            f"{params.exozodi_level_distribution!r}."
+        )
+    median = float(params.exozodi_hosts2020_median_zodi)
+    if median <= 0.0:
+        sampled = 0.0
+    elif float(params.exozodi_hosts2020_sigma_dex) == 0.0:
+        sampled = median
+    else:
+        sampled = median * 10.0 ** float(
+            rng.normal(0.0, float(params.exozodi_hosts2020_sigma_dex))
+        )
+    clipped = max(float(params.exozodi_level_min_zodi), float(sampled))
+    if params.exozodi_level_max_zodi is not None:
+        clipped = min(float(params.exozodi_level_max_zodi), clipped)
+    return clipped, "hosts2020_lognormal"
+
+
+def _resolve_exozodi_morphology(
+    params: SpaceZodiBackground.Params,
+    rng: np.random.Generator,
+) -> dict[str, float | str]:
+    """Resolve simple exozodi morphology diagnostics and rate multiplier."""
+
+    if params.exozodi_morphology == "smooth_face_on":
+        return {
+            "inclination_deg": 0.0,
+            "inclination_source": "smooth_face_on",
+            "raw_surface_brightness_multiplier": 1.0,
+            "surface_brightness_multiplier": 1.0,
+        }
+    if params.exozodi_morphology != "smooth_inclined_disk":
+        raise ValueError(
+            "background.space_zodi: unsupported exozodi_morphology="
+            f"{params.exozodi_morphology!r}."
+        )
+
+    if params.exozodi_inclination_distribution == "isotropic":
+        max_inclination_rad = math.radians(float(params.exozodi_max_inclination_deg))
+        cos_i = float(rng.uniform(math.cos(max_inclination_rad), 1.0))
+        inclination_deg = math.degrees(math.acos(cos_i))
+        inclination_source = "isotropic"
+    elif params.exozodi_inclination_distribution == "fixed":
+        inclination_deg = float(params.exozodi_inclination_deg)
+        inclination_source = "params.exozodi_inclination_deg"
+    else:
+        raise ValueError(
+            "background.space_zodi: unsupported "
+            "exozodi_inclination_distribution="
+            f"{params.exozodi_inclination_distribution!r}."
+        )
+
+    cos_i = max(math.cos(math.radians(inclination_deg)), 1.0e-6)
+    raw_multiplier = 1.0 / cos_i
+    multiplier = min(
+        raw_multiplier,
+        float(params.exozodi_max_surface_brightness_multiplier),
+    )
+    return {
+        "inclination_deg": inclination_deg,
+        "inclination_source": inclination_source,
+        "raw_surface_brightness_multiplier": raw_multiplier,
+        "surface_brightness_multiplier": multiplier,
+    }
 
 
 def _surface_brightness_rate_per_bin_hz(