Skip to content

API Reference

Generated directly from the docstrings in the source — if a signature here disagrees with the prose guides, the source (and therefore this page) wins.

sed_model.grid

sed_model.grid

sed_model.grid

Loads a SED_Tools-prepared stellar atmosphere grid (flux_cube.bin + lookup_table.csv) into memory and exposes the axes and flux array for downstream interpolation.

Binary format (flux_cube.bin)

Header : 4 × int32 -> (n_teff, n_logg, n_meta, n_lambda) Axes : n_teff + n_logg + n_meta + n_lambda float64 values (written consecutively, no padding) Payload : n_teff × n_logg × n_meta × n_lambda float64 values stored in (W, M, L, T) order on disk, loaded into (T, L, M, W).

This matches the layout written by SED_Tools and read by the MESA colors module (colors_utils.f90 :: load_flux_cube).

AtmosphereGrid dataclass 

AtmosphereGrid(teff_grid, logg_grid, meta_grid, wavelengths, flux, model_dir)

Immutable container for a loaded stellar atmosphere grid.

Attributes:

Name Type Description
teff_grid (ndarray, shape(n_T))

Effective temperature grid points in Kelvin.
logg_grid (ndarray, shape(n_L))

log10(g / cm s^-2) grid points.
meta_grid (ndarray, shape(n_M))

Metallicity [M/H] grid points.
wavelengths (ndarray, shape(n_W))

Wavelength grid in Angstroms.
flux (ndarray, shape(n_T, n_L, n_M, n_W))

Surface flux in erg/s/cm^2/Å at each (Teff, logg, [M/H]) node.
model_dir Path

Directory from which the grid was loaded.

in_bounds 

in_bounds(teff, logg, meta)

Return True if (teff, logg, meta) is within the grid axes.

Source code in sed_model/grid.py 
76
77
78
79
80
81
82
def in_bounds(self, teff: float, logg: float, meta: float) -> bool:
    """Return True if (teff, logg, meta) is within the grid axes."""
    return (
        self.teff_bounds[0] <= teff <= self.teff_bounds[1]
        and self.logg_bounds[0] <= logg <= self.logg_bounds[1]
        and self.meta_bounds[0] <= meta <= self.meta_bounds[1]
    )

clamp 

clamp(teff, logg, meta)

Clamp (teff, logg, meta) to the grid boundary.

Source code in sed_model/grid.py 
84
85
86
87
88
89
90
91
def clamp(
    self, teff: float, logg: float, meta: float
) -> tuple[float, float, float]:
    """Clamp (teff, logg, meta) to the grid boundary."""
    teff = float(np.clip(teff, *self.teff_bounds))
    logg = float(np.clip(logg, *self.logg_bounds))
    meta = float(np.clip(meta, *self.meta_bounds))
    return teff, logg, meta

interp_radius 

interp_radius(teff, logg, meta)

Euclidean distance in normalised parameter space from the nearest grid point. Mirrors the Interp_rad diagnostic produced by the MESA colors module.

Source code in sed_model/grid.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def interp_radius(
    self, teff: float, logg: float, meta: float
) -> float:
    """Euclidean distance in normalised parameter space from the
    nearest grid point.  Mirrors the ``Interp_rad`` diagnostic
    produced by the MESA colors module."""
    t_norm = (self.teff_grid[-1] - self.teff_grid[0]) or 1.0
    l_norm = (self.logg_grid[-1] - self.logg_grid[0]) or 1.0
    m_norm = (self.meta_grid[-1] - self.meta_grid[0]) or 1.0

    i_t = int(np.argmin(np.abs(self.teff_grid - teff)))
    i_l = int(np.argmin(np.abs(self.logg_grid - logg)))
    i_m = int(np.argmin(np.abs(self.meta_grid - meta)))

    dt = (teff - self.teff_grid[i_t]) / t_norm
    dl = (logg - self.logg_grid[i_l]) / l_norm
    dm = (meta - self.meta_grid[i_m]) / m_norm

    return float(np.sqrt(dt**2 + dl**2 + dm**2))

load_grid 

load_grid(model_dir)

Load a SED_Tools atmosphere grid from model_dir.

The directory must contain: - flux_cube.bin — binary flux cube - lookup_table.csv — parameter index (used for validation only)

Parameters:

Name Type Description Default
model_dir str | Path

Path to the directory produced by sed-tools rebuild.

 required

Returns:

Type Description
AtmosphereGrid

Raises:

Type Description
FileNotFoundError

If flux_cube.bin is missing.
ValueError

If the binary header is invalid or the file is truncated.
Source code in sed_model/grid.py 
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
def load_grid(model_dir: str | Path) -> AtmosphereGrid:
    """Load a SED_Tools atmosphere grid from *model_dir*.

    The directory must contain:
      - ``flux_cube.bin``   — binary flux cube
      - ``lookup_table.csv`` — parameter index (used for validation only)

    Parameters
    ----------
    model_dir:
        Path to the directory produced by ``sed-tools rebuild``.

    Returns
    -------
    AtmosphereGrid

    Raises
    ------
    FileNotFoundError
        If ``flux_cube.bin`` is missing.
    ValueError
        If the binary header is invalid or the file is truncated.
    """
    model_dir = Path(model_dir).resolve()
    cube_path = model_dir / "flux_cube.bin"

    if not cube_path.exists():
        raise FileNotFoundError(
            f"flux_cube.bin not found in {model_dir}. "
            "Run 'sed-tools rebuild' to generate it."
        )

    teff_grid, logg_grid, meta_grid, wavelengths, flux = _read_flux_cube(cube_path)

    return AtmosphereGrid(
        teff_grid=teff_grid,
        logg_grid=logg_grid,
        meta_grid=meta_grid,
        wavelengths=wavelengths,
        flux=flux,
        model_dir=model_dir,
    )

validate_lookup_table 

validate_lookup_table(grid)

Read and return the lookup_table.csv from grid.model_dir.

Useful for sanity-checking that the grid axes in the binary cube are consistent with the CSV metadata. Not called at load time to avoid the pandas overhead on every initialisation.

Source code in sed_model/grid.py 
236
237
238
239
240
241
242
243
244
245
246
247
248
249
def validate_lookup_table(grid: AtmosphereGrid) -> pd.DataFrame:
    """Read and return the ``lookup_table.csv`` from *grid.model_dir*.

    Useful for sanity-checking that the grid axes in the binary cube
    are consistent with the CSV metadata.  Not called at load time to
    avoid the pandas overhead on every initialisation.
    """
    lookup_path = grid.model_dir / "lookup_table.csv"
    if not lookup_path.exists():
        raise FileNotFoundError(f"lookup_table.csv not found in {grid.model_dir}")

    df = pd.read_csv(lookup_path, dtype=str)
    df.columns = [c.strip().lstrip("#").strip() for c in df.columns]
    return df

sed_model.filters

sed_model.filters

sed_model.filters

Loads photometric filter transmission curves from SED_Tools-prepared .dat files and precomputes Vega, AB, and ST zero-points using the same photon-counting integrals as the MESA colors module (colors/private/synthetic.f90).

Zero-point definitions

Vega : F_zp = ∫ F_vega(λ) T(λ) λ dλ / ∫ T(λ) λ dλ AB : F_zp = ∫ F_AB(λ) T(λ) λ dλ / ∫ T(λ) λ dλ where F_AB(λ) = 3.631e-20 × c / λ² [erg/s/cm²/Å] ST : F_zp = 3.63e-9 (flat F_lambda, constant)

All integrals use the trapezoid rule for wavelength grids that may be non-uniform. The Fortran runtime uses adaptive Simpson; the difference is negligible for the dense wavelength grids produced by SED_Tools.

Filter dataclass 

Filter(name, path, wavelengths, transmission, vega_zero_point=-1.0, ab_zero_point=-1.0, st_zero_point=-1.0)

A single photometric filter with precomputed zero-points.

Attributes:

Name Type Description
name str

Short filter identifier (filename stem, e.g. 'B', 'Gbp').
path Path

Absolute path to the source .dat file.
wavelengths (ndarray, shape(n))

Filter wavelength grid in Angstroms.
transmission (ndarray, shape(n))

Dimensionless transmission in [0, 1].
vega_zero_point float

Photon-counting flux zero-point for the Vega magnitude system. -1.0 if no Vega SED was supplied at load time.
ab_zero_point float

Photon-counting flux zero-point for the AB magnitude system.
st_zero_point float

Photon-counting flux zero-point for the ST magnitude system.

zero_point 

zero_point(system)

Return the zero-point for system ('Vega', 'AB', or 'ST').

Source code in sed_model/filters.py 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
def zero_point(self, system: str) -> float:
    """Return the zero-point for *system* ('Vega', 'AB', or 'ST')."""
    s = system.upper()
    if s == "VEGA":
        if self.vega_zero_point < 0:
            raise ValueError(
                f"Filter '{self.name}': Vega zero-point not available. "
                "Supply a Vega SED path to load_filters()."
            )
        return self.vega_zero_point
    if s == "AB":
        return self.ab_zero_point
    if s == "ST":
        return self.st_zero_point
    raise ValueError(f"Unknown magnitude system '{system}'. Choose Vega, AB, or ST.")

load_filters 

load_filters(filter_paths, vega_sed_path=None)

Load a list of filter .dat files and return a list of :class:Filter objects with precomputed zero-points.

Parameters:

Name Type Description Default
filter_paths list[str | Path]

Paths to individual *.dat filter transmission files in the SED_Tools two-column format (wavelength Å, transmission 0–1). Comments starting with # are skipped.

 required
vega_sed_path str | Path | None

Optional path to a Vega reference SED CSV file (wavelength,flux header, wavelength in Å, flux in erg/s/cm²/Å). Required for Vega zero-points.

 None

Returns:

Type Description
list of Filter, in the same order as *filter_paths*.
Source code in sed_model/filters.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
def load_filters(
    filter_paths: list[str | Path],
    vega_sed_path: str | Path | None = None,
) -> list[Filter]:
    """Load a list of filter ``.dat`` files and return a list of
    :class:`Filter` objects with precomputed zero-points.

    Parameters
    ----------
    filter_paths:
        Paths to individual ``*.dat`` filter transmission files in the
        SED_Tools two-column format (wavelength Å, transmission 0–1).
        Comments starting with ``#`` are skipped.
    vega_sed_path:
        Optional path to a Vega reference SED CSV file
        (``wavelength,flux`` header, wavelength in Å,
        flux in erg/s/cm²/Å).  Required for Vega zero-points.

    Returns
    -------
    list of Filter, in the same order as *filter_paths*.
    """
    vega_wave: np.ndarray | None = None
    vega_flux: np.ndarray | None = None
    if vega_sed_path is not None:
        vega_wave, vega_flux = _load_vega_sed(Path(vega_sed_path))

    filters: list[Filter] = []
    for p in filter_paths:
        p = Path(p).resolve()
        wave, trans = _load_filter_dat(p)
        filt = Filter(
            name=p.stem,
            path=p,
            wavelengths=wave,
            transmission=trans,
        )
        filt.ab_zero_point = _compute_ab_zero_point(wave, trans)
        filt.st_zero_point = _compute_st_zero_point(wave, trans)
        if vega_wave is not None and vega_flux is not None:
            filt.vega_zero_point = _compute_vega_zero_point(
                vega_wave, vega_flux, wave, trans
            )
        filters.append(filt)

    return filters

load_filters_from_instrument_dir 

load_filters_from_instrument_dir(instrument_dir, vega_sed_path=None)

Load all filters listed in a SED_Tools instrument directory.

The directory must contain an index file whose name matches the last component of instrument_dir (e.g. Johnson/Johnson) that lists one .dat filename per line. This mirrors the structure expected by the MESA colors module.

Parameters:

Name Type Description Default
instrument_dir str | Path

Path to a data/filters/<Facility>/<Instrument>/ directory.

 required
vega_sed_path str | Path | None

Optional Vega reference SED for Vega zero-point computation.

 None

Returns:

Type Description
list of Filter
Source code in sed_model/filters.py 
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
def load_filters_from_instrument_dir(
    instrument_dir: str | Path,
    vega_sed_path: str | Path | None = None,
) -> list[Filter]:
    """Load all filters listed in a SED_Tools instrument directory.

    The directory must contain an index file whose name matches the
    last component of *instrument_dir* (e.g. ``Johnson/Johnson``) that
    lists one ``.dat`` filename per line.  This mirrors the structure
    expected by the MESA colors module.

    Parameters
    ----------
    instrument_dir:
        Path to a ``data/filters/<Facility>/<Instrument>/`` directory.
    vega_sed_path:
        Optional Vega reference SED for Vega zero-point computation.

    Returns
    -------
    list of Filter
    """
    instrument_dir = Path(instrument_dir).resolve()
    index_name = instrument_dir.name
    index_file = instrument_dir / index_name

    if not index_file.exists():
        # Fallback: load every .dat in the directory
        dat_files = sorted(instrument_dir.glob("*.dat"))
        if not dat_files:
            raise FileNotFoundError(
                f"No index file '{index_name}' and no .dat files in {instrument_dir}"
            )
    else:
        dat_files = []
        for line in index_file.read_text().splitlines():
            line = line.strip()
            if not line or line.startswith("#"):
                continue
            dat_files.append(instrument_dir / line)

    return load_filters(dat_files, vega_sed_path=vega_sed_path)

sed_model.forward

sed_model.forward

sed_model.forward

Forward model: stellar parameters → SED + synthetic photometry.

Pipeline

(Teff, logg, [M/H]) → SED interpolation (Fortran, Hermite or linear) → distance dilution (Fortran, (R/d)²) → extinction (Python, optional, sed_extinction) → bolometric (Fortran) → filter convolution (Fortran, per filter) → ForwardResult

All five physical parameters — Teff, logg, [M/H], Av, distance — can be fixed or free, described by a FitParams object (sed_model.params). This shared vocabulary is what makes the module bidirectional: the inverse model unpacks a theta vector with FitParams.unpack and passes the result directly to run_forward.

Backward-compatible call

The original positional signature still works::

run_forward(teff, logg, meta, R, d, grid, filters)

FitParams call (used by the inverse model)::

run_forward(fit_params=params, theta=theta_vec,
            R=R_sun, grid=grid, filters=filters)

ForwardResult dataclass 

ForwardResult(wavelengths, surface_flux, observed_flux, magnitudes, band_fluxes, bol_flux, bol_mag, interp_radius, clamped, teff, logg, meta, R, d, a_v=0.0, mag_system='AB', extinction_applied=False)

Output of one forward-model evaluation.

The five physical parameters that produced this result are all stored so the result is self-describing. The inverse model uses ForwardResult.magnitudes to compute the likelihood; the forward model reads parameters from FitParams — both sides work with the same data structure.

run_forward 

run_forward(teff=None, logg=None, meta=None, R=None, d=None, grid=None, filters=None, mag_system='AB', interp_method='hermite', extinction=None, fit_params=None, theta=None)

Evaluate the forward model for one set of stellar parameters.

Two calling conventions are supported.

Classic (backward-compatible)::

run_forward(teff, logg, meta, R, d, grid, filters)

FitParams (used by the inverse model)::

run_forward(fit_params=params, theta=theta_vec,
            R=R_sun, grid=grid, filters=filters)

In the FitParams convention theta contains only the free parameters in canonical order (teff, logg, meta, a_v, d — skipping fixed ones). fit_params.unpack(theta) fills in the fixed values. If Av or distance are free parameters they are taken from theta, not from any keyword argument.

Parameters:

Name Type Description Default
teff float

Atmospheric parameters. Required in classic mode.

 None
logg float

Atmospheric parameters. Required in classic mode.

 None
meta float

Atmospheric parameters. Required in classic mode.

 None
R float

Stellar radius in cm. Always required.

 None
d float

Distance in cm. Required in classic mode; ignored if distance is a free parameter in fit_params.

 None
grid AtmosphereGrid
None
filters list of Filter
None
mag_system ('AB', 'Vega', 'ST')
'AB'
interp_method ('hermite', 'linear')
'hermite'
extinction ExtinctionModel or None

Applied after dilution and before filter convolution. When Av is free in fit_params, the model's stored a_v is overridden by the value from theta at each call.

 None
fit_params FitParams or None
None
theta array - like or None

Required when fit_params is provided.

 None
Source code in sed_model/forward.py 
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
def run_forward(
    # classic positional args
    teff:    Optional[float] = None,
    logg:    Optional[float] = None,
    meta:    Optional[float] = None,
    R:       Optional[float] = None,
    d:       Optional[float] = None,
    grid:    Optional[AtmosphereGrid] = None,
    filters: Optional[list] = None,
    # options
    mag_system:    str = "AB",
    interp_method: str = "hermite",
    extinction:    Optional["ExtinctionModel"] = None,
    # FitParams interface
    fit_params: Optional[FitParams] = None,
    theta:      Optional[np.ndarray] = None,
) -> ForwardResult:
    """Evaluate the forward model for one set of stellar parameters.

    Two calling conventions are supported.

    **Classic** (backward-compatible)::

        run_forward(teff, logg, meta, R, d, grid, filters)

    **FitParams** (used by the inverse model)::

        run_forward(fit_params=params, theta=theta_vec,
                    R=R_sun, grid=grid, filters=filters)

    In the FitParams convention ``theta`` contains only the *free* parameters
    in canonical order (teff, logg, meta, a_v, d — skipping fixed ones).
    ``fit_params.unpack(theta)`` fills in the fixed values.  If Av or
    distance are free parameters they are taken from theta, not from any
    keyword argument.

    Parameters
    ----------
    teff, logg, meta : float
        Atmospheric parameters.  Required in classic mode.
    R : float
        Stellar radius in cm.  Always required.
    d : float
        Distance in cm.  Required in classic mode; ignored if distance is
        a free parameter in fit_params.
    grid : AtmosphereGrid
    filters : list of Filter
    mag_system : {'AB', 'Vega', 'ST'}
    interp_method : {'hermite', 'linear'}
    extinction : ExtinctionModel or None
        Applied after dilution and before filter convolution.
        When Av is free in fit_params, the model's stored a_v is
        overridden by the value from theta at each call.
    fit_params : FitParams or None
    theta : array-like or None
        Required when fit_params is provided.
    """
    # ------------------------------------------------------------------
    # Resolve parameters from whichever calling convention is used
    # ------------------------------------------------------------------
    if fit_params is not None:
        if theta is None:
            raise ValueError("theta must be supplied when fit_params is used")
        p     = fit_params.unpack(np.asarray(theta, dtype=np.float64))
        teff  = float(p['teff'])
        logg  = float(p['logg'])
        meta  = float(p['meta'])
        d_use = float(p['d'])
        av    = float(p['a_v'])
    else:
        if any(x is None for x in (teff, logg, meta, R, d, grid, filters)):
            raise ValueError(
                "teff, logg, meta, R, d, grid, and filters must all be "
                "supplied when fit_params is not used."
            )
        d_use = float(d)
        av    = 0.0

    cc = _get_cc_api()

    # ------------------------------------------------------------------
    # Compute interp_radius and clamped flag via grid helpers
    # (same logic as the original forward.py before the extinction rewrite)
    # ------------------------------------------------------------------
    clamped    = not grid.in_bounds(teff, logg, meta)
    interp_rad = grid.interp_radius(teff, logg, meta)
    teff_q, logg_q, meta_q = grid.clamp(teff, logg, meta)

    # ------------------------------------------------------------------
    # 1. SED interpolation (Fortran)
    #
    # cc.interp_sed_hermite / interp_sed_linear return (result_flux, ierr).
    # The flux array is C-contiguous float64; Fortran expects the cube as
    # (nt, nl, nm, nw) which is grid.flux's native shape.
    # ------------------------------------------------------------------
    flux_cube  = np.asfortranarray(grid.flux, dtype=np.float64)
    teff_grid  = np.ascontiguousarray(grid.teff_grid,   dtype=np.float64)
    logg_grid  = np.ascontiguousarray(grid.logg_grid,   dtype=np.float64)
    meta_grid  = np.ascontiguousarray(grid.meta_grid,   dtype=np.float64)
    wavelengths = np.ascontiguousarray(grid.wavelengths, dtype=np.float64)

    if interp_method == "hermite":
        surface_flux, ierr = cc.interp_sed_hermite(
            teff_q, logg_q, meta_q,
            teff_grid, logg_grid, meta_grid,
            flux_cube,
        )
    elif interp_method == "linear":
        surface_flux, ierr = cc.interp_sed_linear(
            teff_q, logg_q, meta_q,
            teff_grid, logg_grid, meta_grid,
            flux_cube,
        )
    else:
        raise ValueError(
            f"Unknown interp_method '{interp_method}'. "
            "Choose 'hermite' or 'linear'."
        )

    if ierr != 0:
        clamped = True

    # ------------------------------------------------------------------
    # 2. Distance dilution: F_obs = F_surface × (R/d)²  (Fortran)
    # ------------------------------------------------------------------
    observed_flux = cc.dilute_flux(surface_flux, float(R), d_use)

    # ------------------------------------------------------------------
    # 3. Extinction (Python, optional)
    #
    # Applied AFTER dilution, BEFORE filter convolution.
    # Physics: star → distance → dust column → telescope → filter.
    #
    # When Av is a free parameter, rebuild the ExtinctionModel with the
    # Av value from theta so every likelihood call uses the correct value.
    # ------------------------------------------------------------------
    extinction_applied = False
    if extinction is not None and getattr(extinction.config, 'enabled', False):
        if fit_params is not None:
            from dataclasses import replace as _replace
            from .sed_extinction import ExtinctionModel as _EM
            # In FitParams mode, FitParams is the source of truth for Av,
            # whether Av is fixed or free.  This prevents the fixed Av in
            # FitParams and the Av stored in ExtinctionModel from drifting apart.
            new_cfg = _replace(extinction.config, a_v=av)
            extinction = _EM(new_cfg)

        observed_flux      = extinction.apply(wavelengths, observed_flux)
        extinction_applied = True
        av                 = extinction.config.a_v

    # ------------------------------------------------------------------
    # 4. Bolometric quantities (Fortran)
    # ------------------------------------------------------------------
    bol_flux, bol_mag, _ = cc.bolometric(wavelengths, observed_flux)

    # ------------------------------------------------------------------
    # 5. Synthetic photometry per filter (Fortran)
    # ------------------------------------------------------------------
    magnitudes:  dict = {}
    band_fluxes: dict = {}

    for filt in filters:
        zp         = filt.zero_point(mag_system)
        filt_wave  = np.ascontiguousarray(filt.wavelengths,  dtype=np.float64)
        filt_trans = np.ascontiguousarray(filt.transmission, dtype=np.float64)
        mag, band_flux, _ = cc.synthetic_magnitude(
            wavelengths, observed_flux,
            filt_wave, filt_trans,
            zp,
        )
        magnitudes[filt.name]  = float(mag)
        band_fluxes[filt.name] = float(band_flux)

    return ForwardResult(
        wavelengths=wavelengths,
        surface_flux=surface_flux,
        observed_flux=observed_flux,
        magnitudes=magnitudes,
        band_fluxes=band_fluxes,
        bol_flux=float(bol_flux),
        bol_mag=float(bol_mag),
        interp_radius=float(interp_rad),
        clamped=clamped,
        teff=float(teff),
        logg=float(logg),
        meta=float(meta),
        R=float(R),
        d=d_use,
        a_v=av,
        mag_system=mag_system,
        extinction_applied=extinction_applied,
    )

run_forward_batch 

run_forward_batch(params, R, d, grid, filters, mag_system='AB', interp_method='hermite', extinction=None)

Run run_forward over an array of (teff, logg, meta) rows.

Parameters:

Name Type Description Default
params (ndarray, shape(N, 3))
required

Returns:

Type Description
list of ForwardResult, length N.
Source code in sed_model/forward.py 
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
def run_forward_batch(
    params: np.ndarray,
    R: float,
    d: float,
    grid: AtmosphereGrid,
    filters: list,
    mag_system: str = "AB",
    interp_method: str = "hermite",
    extinction: Optional["ExtinctionModel"] = None,
) -> list:
    """Run run_forward over an array of (teff, logg, meta) rows.

    Parameters
    ----------
    params : ndarray, shape (N, 3)
    Returns
    -------
    list of ForwardResult, length N.
    """
    params = np.atleast_2d(params)
    if params.shape[1] != 3:
        raise ValueError("params must have shape (N, 3): teff, logg, meta")
    return [
        run_forward(
            teff=float(r[0]), logg=float(r[1]), meta=float(r[2]),
            R=R, d=d, grid=grid, filters=filters,
            mag_system=mag_system, interp_method=interp_method,
            extinction=extinction,
        )
        for r in params
    ]

sed_model.inverse

sed_model.inverse

sed_model.inverse

Inverse model: observed magnitudes → posterior on stellar parameters.

The inverse model is the mirror of the forward model. Both share the same FitParams object (from sed_model.params) which declares each physical parameter as fixed or free. The MCMC samples only the free parameters; fixed ones are threaded to the forward model unchanged.

This bidirectionality is concrete, not cosmetic: - The forward model maps FitParams.unpack(theta) → magnitudes. - The inverse model maps magnitudes → FitParams + posterior on theta. - The same FitParams instance flows in both directions.

Parameter modes

Every physical quantity — Teff, logg, [M/H], Av, distance — is described by a ParamSpec inside a FitParams. Three behaviours:

fixed(value) Not sampled. Passed unchanged to the forward model at every step.

free(lo, hi) Sampled by MCMC with a flat prior over [lo, hi].

bounded(lo, hi) Alias for free — use it to be explicit about hard limits.

The simplest use::

from sed_model import run_inverse
from sed_model.params import fit_params_from_grid, PC_TO_CM

params = fit_params_from_grid(grid)                # Teff/logg/meta free
result = run_inverse(obs_mags, obs_errs, filter_names,
                     R=R_sun, fit_params=params,
                     grid=grid, filters=filters)

To fit with extinction free::

params = fit_params_from_grid(grid, a_v=(0.0, 3.0))   # Av free 0–3 mag
ext    = ExtinctionModel(enabled=True, law='fitzpatrick99', a_v=0.0)
result = run_inverse(..., fit_params=params, extinction=ext)

To fit with distance free (expect Av–d degeneracy in the posterior)::

params = fit_params_from_grid(
    grid,
    a_v=(0.0, 2.0),
    d_cm=(100*PC_TO_CM, 2000*PC_TO_CM),
)

run_inverse 

run_inverse(obs_magnitudes, obs_uncertainties, filter_names, R, grid, filters, fit_params=None, d=None, extinction_law=None, extinction=None, mag_system='AB', interp_method='hermite', n_walkers=32, n_steps=2000, n_burn=500, n_thin=1, p0_centre=None, p0_teff=None, p0_logg=None, p0_meta=None, p0_scatter=0.02, seed=None, progress=True)

Infer stellar parameters from observed broadband photometry.

Parameters:

Name Type Description Default
obs_magnitudes (array - like, shape(n_filters))

Observed magnitudes in the same order as filter_names.

 required
obs_uncertainties (array - like, shape(n_filters))

1-sigma uncertainties. Must all be > 0.

 required
filter_names list of str

Filter names matching elements of filters.

 required
R float

Stellar radius in cm (always required; not currently sampled).

 required
grid AtmosphereGrid
required
filters list of Filter
required
fit_params FitParams or None

Full parameter specification. If None, one is built from grid using the d convenience argument (default: 1 pc, fixed).

 None
d float or None

Distance in cm. Used only when fit_params is None to build a default FitParams with fixed distance. If fit_params is supplied this argument is ignored — set distance via fit_params.

 None
extinction ExtinctionModel or None

Dust extinction model. Required if Av is free in fit_params. Ignored if Av is fixed at 0.

 None
mag_system ('AB', 'Vega', 'ST')
'AB'
interp_method ('hermite', 'linear')
'hermite'
n_walkers int

emcee settings.

 32
n_steps int

emcee settings.

 32
n_burn int

emcee settings.

 32
n_thin int

emcee settings.

 32
p0_centre dict or None

Optional starting point keyed by parameter name, e.g. {'teff': 5800, 'logg': 4.4, 'meta': 0.0}.

 None
p0_scatter float

Width of initial ball as fraction of each parameter's range.

 0.02
seed int or None
None
progress bool
True

Returns:

Type Description
InverseResult
Source code in sed_model/inverse.py 
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
def run_inverse(
    obs_magnitudes:    list,
    obs_uncertainties: list,
    filter_names:      list,
    R:                 float,
    grid:              AtmosphereGrid,
    filters:           list,
    # FitParams — the shared parameter spec
    fit_params:        Optional[FitParams] = None,
    # convenience: if fit_params not supplied, these are used to build one
    d:                 Optional[float] = None,
    extinction_law:    Optional[str]   = None,
    # extinction model (required if a_v is free or if you want reddening)
    extinction:        Optional["ExtinctionModel"] = None,
    # sampler settings
    mag_system:        str   = "AB",
    interp_method:     str   = "hermite",
    n_walkers:         int   = 32,
    n_steps:           int   = 2000,
    n_burn:            int   = 500,
    n_thin:            int   = 1,
    p0_centre:         Optional[dict] = None,
    p0_teff:           Optional[float] = None,
    p0_logg:           Optional[float] = None,
    p0_meta:           Optional[float] = None,
    p0_scatter:        float = 0.02,
    seed:              Optional[int]  = None,
    progress:          bool  = True,
) -> InverseResult:
    """Infer stellar parameters from observed broadband photometry.

    Parameters
    ----------
    obs_magnitudes : array-like, shape (n_filters,)
        Observed magnitudes in the same order as ``filter_names``.
    obs_uncertainties : array-like, shape (n_filters,)
        1-sigma uncertainties.  Must all be > 0.
    filter_names : list of str
        Filter names matching elements of ``filters``.
    R : float
        Stellar radius in cm (always required; not currently sampled).
    grid : AtmosphereGrid
    filters : list of Filter
    fit_params : FitParams or None
        Full parameter specification.  If None, one is built from ``grid``
        using the ``d`` convenience argument (default: 1 pc, fixed).
    d : float or None
        Distance in cm.  Used only when ``fit_params`` is None to build a
        default ``FitParams`` with fixed distance.  If ``fit_params`` is
        supplied this argument is ignored — set distance via fit_params.
    extinction : ExtinctionModel or None
        Dust extinction model.  Required if Av is free in fit_params.
        Ignored if Av is fixed at 0.
    mag_system : {'AB', 'Vega', 'ST'}
    interp_method : {'hermite', 'linear'}
    n_walkers, n_steps, n_burn, n_thin : int
        emcee settings.
    p0_centre : dict or None
        Optional starting point keyed by parameter name, e.g.
        ``{'teff': 5800, 'logg': 4.4, 'meta': 0.0}``.
    p0_scatter : float
        Width of initial ball as fraction of each parameter's range.
    seed : int or None
    progress : bool

    Returns
    -------
    InverseResult
    """
    try:
        import emcee
    except ImportError as exc:
        raise ImportError(
            "emcee is required for the inverse model. "
            "Install with: pip install emcee"
        ) from exc

    # ------------------------------------------------------------------
    # Input validation
    # ------------------------------------------------------------------
    obs_mag = np.asarray(obs_magnitudes,   dtype=np.float64)
    obs_err = np.asarray(obs_uncertainties, dtype=np.float64)

    if obs_mag.shape != obs_err.shape:
        raise ValueError(
            "obs_magnitudes and obs_uncertainties must have the same length."
        )
    if len(filter_names) != len(obs_mag):
        raise ValueError(
            "filter_names must have the same length as obs_magnitudes."
        )
    if not np.all(obs_err > 0):
        raise ValueError(
            "All obs_uncertainties must be positive (> 0). "
            "Zero or negative uncertainties make the likelihood undefined."
        )
    if n_walkers < 6:
        raise ValueError("n_walkers must be >= 6 (2 × n_dim).")
    if n_burn >= n_steps:
        raise ValueError("n_burn must be < n_steps.")

    available = {f.name for f in filters}
    missing   = set(filter_names) - available
    if missing:
        raise ValueError(
            f"Filter names not in the supplied filters: {missing}. "
            f"Available: {available}"
        )

    # ------------------------------------------------------------------
    # Build FitParams if not supplied
    # ------------------------------------------------------------------
    if fit_params is None:
        d_cm = d if d is not None else PC_TO_CM
        fit_params = fit_params_from_grid(grid, d_cm=d_cm)

    n_dim = fit_params.n_free
    if n_dim == 0:
        raise ValueError(
            "All parameters are fixed — nothing to sample. "
            "Make at least one parameter free in fit_params."
        )
    if n_walkers < 2 * n_dim:
        raise ValueError(
            f"n_walkers ({n_walkers}) must be >= 2 × n_free ({2*n_dim}). "
            f"Free parameters: {fit_params.free_names}"
        )

    # ------------------------------------------------------------------
    # Build a default extinction model if Av is active but the caller did
    # not provide one.
    # ------------------------------------------------------------------
    av_active = fit_params.a_v.is_free or (fit_params.a_v.is_fixed and fit_params.a_v.value > 0.0)
    if extinction is None and av_active:
        from .sed_extinction import make_extinction_model
        extinction = make_extinction_model(
            enabled=True,
            law=extinction_law or 'fitzpatrick99',
            a_v=0.0 if fit_params.a_v.is_free else fit_params.a_v.value,
        )

    # ------------------------------------------------------------------
    # Log what is being sampled vs fixed
    # ------------------------------------------------------------------
    print(fit_params.summary())
    if extinction is not None and getattr(extinction.config, 'enabled', False):
        cfg = extinction.config
        av_mode = "free (sampled)" if fit_params.a_v.is_free else f"fixed={fit_params.a_v.value:.3f}"
        print(
            f"[inverse] Extinction: law={cfg.law}, Rv={cfg.r_v:.2f}, "
            f"Av {av_mode}"
        )
    else:
        print("[inverse] Extinction: disabled")

    # ------------------------------------------------------------------
    # Initial walker positions
    # Merge individual p0_teff/logg/meta kwargs into p0_centre dict.
    # p0_centre takes precedence; individual kwargs are a convenience alias.
    # ------------------------------------------------------------------
    centre: dict = {}
    if p0_teff is not None: centre['teff'] = p0_teff
    if p0_logg is not None: centre['logg'] = p0_logg
    if p0_meta is not None: centre['meta'] = p0_meta
    if p0_centre:
        centre.update(p0_centre)   # p0_centre wins on conflict

    rng = np.random.default_rng(seed)
    p0  = fit_params.initial_ball(
        n_walkers, centre=centre or None, scatter=p0_scatter, rng=rng
    )

    # ------------------------------------------------------------------
    # Run MCMC
    # ------------------------------------------------------------------
    sampler = emcee.EnsembleSampler(
        n_walkers, n_dim, _log_posterior,
        args=(
            obs_mag, obs_err, list(filter_names),
            R, fit_params, grid, filters,
            mag_system, interp_method, extinction,
        ),
    )
    sampler.run_mcmc(p0, n_steps, progress=progress)

    # Acceptance fraction diagnostic
    acc = sampler.acceptance_fraction
    if np.any(acc < 0.1) or np.any(acc > 0.9):
        warnings.warn(
            f"Some walkers have extreme acceptance fractions "
            f"(min={acc.min():.2f}, max={acc.max():.2f}). "
            f"Consider adjusting n_walkers or p0_scatter.",
            UserWarning, stacklevel=2,
        )

    try:
        tau = sampler.get_autocorr_time(quiet=True)
    except Exception:
        tau = None

    flat_samples  = sampler.get_chain(discard=n_burn, thin=n_thin, flat=True)
    flat_log_prob = sampler.get_log_prob(discard=n_burn, thin=n_thin, flat=True)

    fixed_params = {name: fit_params._spec(name).value for name in fit_params.fixed_names}

    return InverseResult(
        samples=flat_samples,
        log_prob=flat_log_prob,
        filter_names=list(filter_names),
        obs_magnitudes=obs_mag,
        obs_uncertainties=obs_err,
        R=float(R),
        d=float(fit_params.d.value if fit_params.d.is_fixed else np.nan),
        mag_system=mag_system,
        n_walkers=n_walkers,
        n_steps=n_steps,
        n_burn=n_burn,
        n_thin=n_thin,
        acceptance_fraction=acc,
        autocorr_time=tau,
        param_names=list(fit_params.free_names),
        fixed_params=fixed_params,
    )

sed_model.params

sed_model.params

sed_model.params

Shared parameter definitions used by both the forward and inverse models.

This is the single place that knows what a "stellar parameter" is, whether it is being sampled or held fixed, and what its physical bounds are. Having it here — rather than scattered across forward.py and inverse.py — is what makes the module genuinely bidirectional: both directions speak the same language about the same quantities.

Parameter modes

Every physical parameter (Teff, logg, [M/H], Av, distance) is described by a ParamSpec. Three modes are supported:

fixed(value) The parameter is not sampled. It is passed directly to the forward model at every likelihood evaluation. This is the default for Av and distance.

free(lo, hi) The parameter is sampled by the MCMC over the interval [lo, hi]. A flat prior is used. This is the default for Teff, logg, and [M/H] (bounds taken from the atmosphere grid at construction time).

bounded(lo, hi) Alias for free — included so user code can be explicit that it wants the parameter sampled with a hard upper and lower limit.

The degeneracy between Av and distance is real and expected: if both are free simultaneously the posterior will be correlated. That is honest — the data genuinely cannot break the degeneracy without additional information. The MCMC will show it.

Usage

::

from sed_model.params import ParamSpec, FitParams, fixed, free

params = FitParams(
    teff = free(4000, 8000),
    logg = free(3.0, 5.0),
    meta = free(-1.0, 0.5),
    a_v  = fixed(0.3),           # Av fixed, not sampled
    d    = fixed(500 * PC_TO_CM), # distance fixed
)

# In the forward model:
result = run_forward(..., fit_params=params, theta=[5800, 4.4, 0.0])
# theta contains only the FREE parameters in the order teff, logg, meta, [av, d]

# Inspect what is being sampled:
print(params.free_names)   # ['teff', 'logg', 'meta']
print(params.n_free)       # 3

ParamSpec dataclass 

ParamSpec(name, mode, value=None, lo=None, hi=None)

Specification for a single physical parameter.

Attributes:

Name Type Description
name str

Human-readable parameter name (e.g. 'Teff', 'Av').
mode {'fixed', 'free'}

Whether the parameter is held constant or sampled.
value float or None

Fixed value. Only meaningful when mode == 'fixed'.
lo, hi float or None

Sampling bounds. Only meaningful when mode == 'free'.

contains 

contains(v)

True if v is within the sampling bounds (free) or equals the fixed value.

Source code in sed_model/params.py 
121
122
123
124
125
def contains(self, v: float) -> bool:
    """True if *v* is within the sampling bounds (free) or equals the fixed value."""
    if self.is_fixed:
        return True  # fixed values are always "valid" from a prior standpoint
    return self.lo <= v <= self.hi

FitParams dataclass 

FitParams(teff, logg, meta, a_v=(lambda: fixed(0.0, name='a_v'))(), d=(lambda: fixed(PC_TO_CM, name='d'))())

Complete specification of all physical parameters for a fit.

Parameters:

Name Type Description Default
teff ParamSpec

Effective temperature in Kelvin.

 required
logg ParamSpec

Log surface gravity (log10 g / cm s-2).

 required
meta ParamSpec

Metallicity [M/H] in dex.

 required
a_v ParamSpec

V-band extinction in magnitudes. Default: fixed(0.0).

 (lambda: fixed(0.0, name='a_v'))()
d ParamSpec

Distance in cm. Default: fixed(1 pc).

 (lambda: fixed(PC_TO_CM, name='d'))()
Notes

The canonical theta vector used by the MCMC sampler contains only the free parameters, in the order: teff, logg, meta, a_v, d (skipping any that are fixed). Use pack / unpack to convert between the full parameter space and the reduced theta vector.

free_names property 

free_names

Names of the free parameters in canonical order.

n_free property 

n_free

Number of free parameters.

fixed_names property 

fixed_names

Names of the fixed parameters.

pack 

pack(teff, logg, meta, a_v=None, d=None)

Pack full physical values into the reduced free-parameter theta vector.

Only the free parameters are included, in canonical order. Fixed parameters are ignored (the values from the spec are used at unpack time).

Source code in sed_model/params.py 
235
236
237
238
239
240
241
242
243
244
245
246
def pack(self, teff: float, logg: float, meta: float,
         a_v: Optional[float] = None, d: Optional[float] = None) -> np.ndarray:
    """Pack full physical values into the reduced free-parameter theta vector.

    Only the free parameters are included, in canonical order.
    Fixed parameters are ignored (the values from the spec are used at
    unpack time).
    """
    full = dict(teff=teff, logg=logg, meta=meta,
                a_v=a_v if a_v is not None else self.a_v.value,
                d=d if d is not None else self.d.value)
    return np.array([full[p] for p in self.free_names], dtype=np.float64)

unpack 

unpack(theta)

Expand the reduced theta vector into a full {param: value} dict.

Fixed parameters are filled in from their specs. Free parameters are taken from theta in canonical order.

Source code in sed_model/params.py 
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
def unpack(self, theta: np.ndarray) -> dict:
    """Expand the reduced theta vector into a full {param: value} dict.

    Fixed parameters are filled in from their specs.  Free parameters
    are taken from theta in canonical order.
    """
    if len(theta) != self.n_free:
        raise ValueError(
            f"theta has {len(theta)} elements but {self.n_free} are free "
            f"({self.free_names})"
        )
    result = {}
    free_iter = iter(theta)
    for p in _PARAM_ORDER:
        spec = self._spec(p)
        result[p] = next(free_iter) if spec.is_free else spec.value
    return result

in_prior 

in_prior(theta)

Return True if all free parameters in theta are within their bounds.

Source code in sed_model/params.py 
266
267
268
269
270
271
272
273
274
275
def in_prior(self, theta: np.ndarray) -> bool:
    """Return True if all free parameters in theta are within their bounds."""
    free_iter = iter(theta)
    for p in _PARAM_ORDER:
        spec = self._spec(p)
        if spec.is_free:
            v = next(free_iter)
            if not (spec.lo <= v <= spec.hi):
                return False
    return True

initial_ball 

initial_ball(n_walkers, centre=None, scatter=0.02, rng=None)

Generate initial walker positions in a Gaussian ball.

Parameters:

Name Type Description Default
n_walkers int

Number of walkers.

 required
centre dict or None

Central values keyed by parameter name. If None, the midpoint of the free bounds is used for free parameters.

 None
scatter float

Width of the ball as a fraction of the parameter range.

 0.02
rng Generator or None
None

Returns:

Name Type Description
pos (ndarray, shape(n_walkers, n_free))
Source code in sed_model/params.py 
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
def initial_ball(self, n_walkers: int,
                 centre: Optional[dict] = None,
                 scatter: float = 0.02,
                 rng: Optional[np.random.Generator] = None) -> np.ndarray:
    """Generate initial walker positions in a Gaussian ball.

    Parameters
    ----------
    n_walkers : int
        Number of walkers.
    centre : dict or None
        Central values keyed by parameter name.  If None, the midpoint
        of the free bounds is used for free parameters.
    scatter : float
        Width of the ball as a fraction of the parameter range.
    rng : Generator or None

    Returns
    -------
    pos : ndarray, shape (n_walkers, n_free)
    """
    if rng is None:
        rng = np.random.default_rng()

    c = []
    s = []
    lo_arr = []
    hi_arr = []
    for p in self.free_names:
        spec = self._spec(p)
        mid = 0.5 * (spec.lo + spec.hi)
        c.append(centre[p] if (centre and p in centre) else mid)
        rng_width = spec.hi - spec.lo
        s.append(scatter * rng_width)
        lo_arr.append(spec.lo)
        hi_arr.append(spec.hi)

    c  = np.array(c)
    s  = np.array(s)
    lo = np.array(lo_arr)
    hi = np.array(hi_arr)

    pos = c + s * rng.standard_normal((n_walkers, self.n_free))
    pos = np.clip(pos, lo, hi)
    return pos

fixed 

fixed(value, name='')

Return a fixed ParamSpec.

Source code in sed_model/params.py 
137
138
139
def fixed(value: float, name: str = "") -> ParamSpec:
    """Return a fixed ParamSpec."""
    return ParamSpec(name=name, mode='fixed', value=float(value))

free 

free(lo, hi, name='')

Return a free (sampled) ParamSpec with flat prior over [lo, hi].

Source code in sed_model/params.py 
142
143
144
def free(lo: float, hi: float, name: str = "") -> ParamSpec:
    """Return a free (sampled) ParamSpec with flat prior over [lo, hi]."""
    return ParamSpec(name=name, mode='free', lo=float(lo), hi=float(hi))

fit_params_from_grid 

fit_params_from_grid(grid, a_v=None, d_cm=None, *, teff=None, logg=None, meta=None)

Build a FitParams from grid bounds plus optional user overrides.

By default Teff/logg/[M/H] are free over the full grid, Av is fixed to zero, and distance is fixed to 1 pc. For any parameter, pass either a float to fix it or a (lo, hi) tuple to sample it.

Examples:

Fit Teff/logg/meta, with Av fixed::

fit_params_from_grid(grid, a_v=0.3, d_cm=500*PC_TO_CM)

Fit only Teff and Av::

fit_params_from_grid(
    grid,
    teff=(5200, 6400),
    logg=4.4,
    meta=0.0,
    a_v=(0.0, 1.5),
    d_cm=500*PC_TO_CM,
)
Source code in sed_model/params.py 
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
def fit_params_from_grid(
    grid,
    a_v:  Union[float, Tuple[float, float], ParamSpec, None] = None,
    d_cm: Union[float, Tuple[float, float], ParamSpec, None] = None,
    *,
    teff: Union[float, Tuple[float, float], ParamSpec, None] = None,
    logg: Union[float, Tuple[float, float], ParamSpec, None] = None,
    meta: Union[float, Tuple[float, float], ParamSpec, None] = None,
) -> FitParams:
    """Build a FitParams from grid bounds plus optional user overrides.

    By default Teff/logg/[M/H] are free over the full grid, Av is fixed to
    zero, and distance is fixed to 1 pc.  For any parameter, pass either a
    float to fix it or a ``(lo, hi)`` tuple to sample it.

    Examples
    --------
    Fit Teff/logg/meta, with Av fixed::

        fit_params_from_grid(grid, a_v=0.3, d_cm=500*PC_TO_CM)

    Fit only Teff and Av::

        fit_params_from_grid(
            grid,
            teff=(5200, 6400),
            logg=4.4,
            meta=0.0,
            a_v=(0.0, 1.5),
            d_cm=500*PC_TO_CM,
        )
    """
    teff_default = free(grid.teff_bounds[0], grid.teff_bounds[1], name='teff')
    logg_default = free(grid.logg_bounds[0], grid.logg_bounds[1], name='logg')
    meta_default = free(grid.meta_bounds[0], grid.meta_bounds[1], name='meta')
    av_default   = fixed(0.0, name='a_v')
    d_default    = fixed(PC_TO_CM, name='d')

    return FitParams(
        teff=_spec_from_user_value('teff', teff, teff_default),
        logg=_spec_from_user_value('logg', logg, logg_default),
        meta=_spec_from_user_value('meta', meta, meta_default),
        a_v=_spec_from_user_value('a_v', a_v, av_default),
        d=_spec_from_user_value('d', d_cm, d_default),
    )

sed_model.io

sed_model.io

sed_model.io

Result containers and persistence for SED_Model outputs.

InverseResult Holds the full emcee sampler output, summary statistics, and the input observations. Can be saved to / loaded from a NumPy .npz archive so posteriors are portable without pickling the sampler.

ForwardResult serialisation Thin helpers to write a ForwardResult SED to a two-column CSV.

InverseResult dataclass 

InverseResult(samples, log_prob, filter_names, obs_magnitudes, obs_uncertainties, R, d, mag_system, n_walkers, n_steps, n_burn, n_thin, acceptance_fraction, autocorr_time=None, param_names=(lambda: ['teff', 'logg', 'meta'])(), fixed_params=dict())

Output of the inverse (photometry → stellar parameters) inference.

Attributes:

Name Type Description
samples (ndarray, shape(n_samples, n_free))

Flattened posterior samples after burn-in. Columns follow param_names (the free parameters in canonical order: teff, logg, meta, a_v, d — skipping fixed ones).
log_prob (ndarray, shape(n_samples))

Log-posterior value for each sample.
filter_names list of str

Names of filters used in the fit, in order.
obs_magnitudes (ndarray, shape(n_filters))

Observed magnitudes supplied by the user.
obs_uncertainties (ndarray, shape(n_filters))

Per-filter magnitude uncertainties (1-sigma).
R float

Stellar radius used (cm).
d float

Distance used (cm).
mag_system str

Photometric system ('Vega', 'AB', 'ST').
n_walkers int

Number of emcee walkers.
n_steps int

Total steps per walker (including burn-in).
n_burn int

Steps discarded as burn-in.
n_thin int

Thinning factor applied to the chain.
acceptance_fraction (ndarray, shape(n_walkers))

Per-walker acceptance fraction.
autocorr_time (ndarray or None, shape(n_free))

Estimated integrated autocorrelation time for each parameter. None if estimation failed (chain too short).

summary 

summary()

Return median and 1-sigma credible interval for each sampled parameter.

Source code in sed_model/io.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
def summary(self) -> dict[str, dict[str, float]]:
    """Return median and 1-sigma credible interval for each sampled parameter."""
    if len(self.param_names) != self.samples.shape[1]:
        raise ValueError(
            f"param_names has {len(self.param_names)} entries but samples has "
            f"{self.samples.shape[1]} columns"
        )

    result = {}
    for i, label in enumerate(self.param_names):
        col = self.samples[:, i]
        med = float(np.percentile(col, 50))
        lo  = float(np.percentile(col, 15.865))
        hi  = float(np.percentile(col, 84.135))
        result[label] = {
            "median":       med,
            "lo":           lo,
            "hi":           hi,
            "lower_1sigma": med - lo,
            "upper_1sigma": hi - med,
        }
    return result

map_estimate 

map_estimate()

Return the maximum a-posteriori sample as a {parameter: value} dict.

Source code in sed_model/io.py 
114
115
116
117
118
def map_estimate(self) -> dict[str, float]:
    """Return the maximum a-posteriori sample as a {parameter: value} dict."""
    idx = int(np.argmax(self.log_prob))
    s = self.samples[idx]
    return {name: float(s[i]) for i, name in enumerate(self.param_names)}

print_summary 

print_summary()

Print a formatted parameter summary for arbitrary free parameters.

Source code in sed_model/io.py 
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
def print_summary(self) -> None:
    """Print a formatted parameter summary for arbitrary free parameters."""
    s = self.summary()
    map_params = self.map_estimate()
    labels = {
        "teff": "Teff   [K]",
        "logg": "logg      ",
        "meta": "[M/H]     ",
        "a_v":  "Av    [mag]",
        "d":    "d      [cm]",
    }

    print(f"\n{'─'*60}")
    print("  SED_Model  —  Posterior Summary")
    print(f"{'─'*60}")
    print(f"  Filters : {', '.join(self.filter_names)}")
    print(f"  System  : {self.mag_system}")
    print(f"  Sampled : {', '.join(self.param_names)}")
    if self.fixed_params:
        fixed = ', '.join(f"{k}={v:.6g}" for k, v in self.fixed_params.items())
        print(f"  Fixed   : {fixed}")
    print(f"  Samples : {len(self.samples)} "
          f"(walkers={self.n_walkers}, steps={self.n_steps}, "
          f"burn={self.n_burn}, thin={self.n_thin})")
    af_mean = float(np.mean(self.acceptance_fraction))
    print(f"  Mean acceptance fraction : {af_mean:.3f}")
    if self.autocorr_time is not None:
        tau = ' / '.join(f"{x:.1f}" for x in np.ravel(self.autocorr_time))
        print(f"  Autocorr times : {tau}")

    print(f"{'─'*60}")
    for param in self.param_names:
        v = s[param]
        label = labels.get(param, param)
        print(f"  {label:<12s} :  {v['median']:>10.3f}"
              f"  +{v['upper_1sigma']:.3f}  -{v['lower_1sigma']:.3f}")

    map_text = ', '.join(f"{k}={v:.6g}" for k, v in map_params.items())
    print(f"\n  MAP : {map_text}")
    print(f"{'─'*60}\n")

save 

save(path)

Save the posterior to a compressed NumPy archive (.npz).

Parameters:

Name Type Description Default
path str or Path

Output file. The .npz extension is added if absent.

 required
Source code in sed_model/io.py 
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
def save(self, path: str | Path) -> None:
    """Save the posterior to a compressed NumPy archive (.npz).

    Parameters
    ----------
    path : str or Path
        Output file.  The ``.npz`` extension is added if absent.
    """
    path = Path(path)
    if path.suffix != ".npz":
        path = path.with_suffix(".npz")

    meta_dict = {
        "filter_names":    self.filter_names,
        "R":               self.R,
        "d":               self.d,
        "mag_system":      self.mag_system,
        "n_walkers":       self.n_walkers,
        "n_steps":         self.n_steps,
        "n_burn":          self.n_burn,
        "n_thin":          self.n_thin,
        "param_names":     self.param_names,
        "fixed_params":    self.fixed_params,
    }

    arrays = dict(
        samples=self.samples,
        log_prob=self.log_prob,
        obs_magnitudes=self.obs_magnitudes,
        obs_uncertainties=self.obs_uncertainties,
        acceptance_fraction=self.acceptance_fraction,
        meta_json=np.array([json.dumps(meta_dict)]),
    )
    if self.autocorr_time is not None:
        arrays["autocorr_time"] = self.autocorr_time

    np.savez_compressed(path, **arrays)

load classmethod 

load(path)

Load an InverseResult saved with :meth:save.

Parameters:

Name Type Description Default
path str or Path

Path to a .npz file produced by :meth:save.

 required
Source code in sed_model/io.py 
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
@classmethod
def load(cls, path: str | Path) -> "InverseResult":
    """Load an InverseResult saved with :meth:`save`.

    Parameters
    ----------
    path : str or Path
        Path to a ``.npz`` file produced by :meth:`save`.
    """
    path = Path(path)
    if path.suffix != ".npz":
        path = path.with_suffix(".npz")

    data = np.load(path, allow_pickle=False)
    meta = json.loads(str(data["meta_json"][0]))

    autocorr = (
        data["autocorr_time"]
        if "autocorr_time" in data
        else None
    )

    return cls(
        samples=data["samples"],
        log_prob=data["log_prob"],
        filter_names=meta["filter_names"],
        obs_magnitudes=data["obs_magnitudes"],
        obs_uncertainties=data["obs_uncertainties"],
        R=float(meta["R"]),
        d=float(meta["d"]),
        mag_system=meta["mag_system"],
        n_walkers=int(meta["n_walkers"]),
        n_steps=int(meta["n_steps"]),
        n_burn=int(meta["n_burn"]),
        n_thin=int(meta["n_thin"]),
        acceptance_fraction=data["acceptance_fraction"],
        autocorr_time=autocorr,
        param_names=meta.get("param_names", ["teff", "logg", "meta"]),
        fixed_params=meta.get("fixed_params", {}),
    )

save_sed 

save_sed(result, path)

Write the observed-frame SED from a ForwardResult to a CSV.

Columns: wavelength (Å), surface_flux, observed_flux (erg/s/cm²/Å).

Parameters:

Name Type Description Default
result ForwardResult
required
path str or Path
required
Source code in sed_model/io.py 
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def save_sed(result: ForwardResult, path: str | Path) -> None:
    """Write the observed-frame SED from a ForwardResult to a CSV.

    Columns: wavelength (Å), surface_flux, observed_flux (erg/s/cm²/Å).

    Parameters
    ----------
    result : ForwardResult
    path : str or Path
    """
    path = Path(path)
    header = (
        f"# SED_Model SED output\n"
        f"# Teff={result.teff:.1f} K  logg={result.logg:.3f}  "
        f"[M/H]={result.meta:.3f}\n"
        f"# R={result.R:.4e} cm  d={result.d:.4e} cm\n"
        f"# wavelength_AA, surface_flux_erg_s_cm2_AA, observed_flux_erg_s_cm2_AA"
    )
    data = np.column_stack([
        result.wavelengths,
        result.surface_flux,
        result.observed_flux,
    ])
    np.savetxt(path, data, delimiter=",", header=header, comments="")

save_magnitudes 

save_magnitudes(result, path)

Write the synthetic magnitudes from a ForwardResult to a CSV.

Columns: filter_name, magnitude, band_flux (erg/s/cm²/Å).

Parameters:

Name Type Description Default
result ForwardResult
required
path str or Path
required
Source code in sed_model/io.py 
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
def save_magnitudes(result: ForwardResult, path: str | Path) -> None:
    """Write the synthetic magnitudes from a ForwardResult to a CSV.

    Columns: filter_name, magnitude, band_flux (erg/s/cm²/Å).

    Parameters
    ----------
    result : ForwardResult
    path : str or Path
    """
    path = Path(path)
    lines = [
        "# SED_Model synthetic photometry",
        f"# Teff={result.teff:.1f} K  logg={result.logg:.3f}  "
        f"[M/H]={result.meta:.3f}  system={result.mag_system}",
        "# filter_name, magnitude, band_flux_erg_s_cm2_AA",
    ]
    for name, mag in result.magnitudes.items():
        bf = result.band_fluxes.get(name, float("nan"))
        lines.append(f"{name},{mag:.6f},{bf:.6e}")

    path.write_text("\n".join(lines) + "\n")

sed_model.sed_extinction

sed_model.sed_extinction

sed_extinction.py

Standalone interstellar dust extinction module for SED fitting.

All extinction laws are implemented in pure NumPy — no compiled Cython dependency. The numerical coefficients are taken directly from the extinction package (Barbary 2016, which itself faithfully reproduces the original papers), plus a native implementation of Gordon et al. (2023).

Usage

Apply extinction to a spectrum before filter convolution::

from sed_extinction import ExtinctionModel

ext = ExtinctionModel(law='fitzpatrick99', r_v=3.1, a_v=0.3,
                      distance_pc=100.0)
flux_obs = ext.apply(wavelength_aa, flux_intrinsic)

Or use it purely as a forward-model modifier inside a fitter::

ext = ExtinctionModel(law='ccm89', r_v=3.1, a_v=0.0,   # disabled
                      distance_pc=1.0, enabled=False)

# later, at evaluation time:
flux_obs = ext.apply(wave, flux)   # returns flux unchanged when disabled
Supported laws

'ccm89' Cardelli, Clayton & Mathis (1989) ApJ 345 245 'odonnell94' O'Donnell (1994) ApJ 422 158 (CCM89 with revised optical) 'fitzpatrick99' Fitzpatrick (1999) PASP 111 63 (R_V-dependent spline) 'fm07' Fitzpatrick & Massa (2007) ApJ 663 320 (R_V = 3.1 fixed) 'calzetti00' Calzetti et al. (2000) ApJ 533 682 (starburst galaxies) 'gordon23' Gordon et al. (2023) ApJ 950 86 (piecewise, MW/SMC/LMC)

Distance scaling

The module also applies distance dimming so that F_obs = F_intrinsic * (R_star / distance)^2 but ONLY when scale_distance is True (default False). Both Av and distance are not fitted by default — they are fixed parameters passed at construction time. To use them in a fit set enabled=True and pass the appropriate a_v / distance_pc when constructing the model for each likelihood evaluation.

Convention

All wavelengths are in Angstroms. Fluxes are in any linear unit (erg/s/cm²/Å, normalised, etc.) — only the ratio matters.

The extinction is applied as: F_obs(λ) = F_intrinsic(λ) × 10^(−0.4 × A(λ)) where A(λ) = a_v × k(λ) and k(λ) is the normalised extinction curve returned by the chosen law.

ExtinctionConfig dataclass 

ExtinctionConfig(enabled=False, law='fitzpatrick99', r_v=3.1, a_v=0.0, gordon23_env='mw', distance_pc=1.0, scale_distance=False)

Configuration container for extinction + distance settings.

All parameters are fixed (not fitted) unless you rebuild the object for each likelihood evaluation in your fitter.

Parameters:

Name Type Description Default
enabled bool

If False, ExtinctionModel.apply() is a no-op. Default False.

 False
law str

Name of the extinction law. One of AVAILABLE_LAWS.

 'fitzpatrick99'
r_v float

R_V value (ignored by 'fm07' which fixes R_V = 3.1, and by 'calzetti00' which uses its own R_V = 4.05 unless overridden).

 3.1
a_v float

V-band extinction in magnitudes. Default 0.0.

 0.0
gordon23_env str

Only used when law='gordon23'. One of 'mw', 'lmc', 'smc'.

 'mw'
distance_pc float

Distance in parsecs. Default 1.0 (no dilution applied unless scale_distance is True and a stellar radius is provided).

 1.0
scale_distance bool

If True, multiply the flux by (R_star_cm / distance_cm)^2 where R_star_cm must be supplied to ExtinctionModel.apply(). Default False.

 False

ExtinctionModel 

ExtinctionModel(config=None, **kwargs)

Apply extinction and (optionally) distance scaling to a spectrum.

Parameters:

Name Type Description Default
config ExtinctionConfig

Full configuration object. If None, a default (disabled) config is used.

 None
**kwargs

Convenience keyword arguments forwarded to ExtinctionConfig. These override fields in config if both are supplied.

 {}

Examples:

Disabled (default) — passes flux through unchanged::

ext = ExtinctionModel()
flux_out = ext.apply(wave, flux)

Enabled, Fitzpatrick99, Av=0.5 at 200 pc::

ext = ExtinctionModel(enabled=True, law='fitzpatrick99',
                      a_v=0.5, r_v=3.1, distance_pc=200.)
flux_obs = ext.apply(wave, flux)

Gordon+2023 SMC curve::

ext = ExtinctionModel(enabled=True, law='gordon23',
                      a_v=0.8, gordon23_env='smc')
flux_obs = ext.apply(wave, flux)
Source code in sed_model/sed_extinction.py 
708
709
710
711
712
713
714
715
716
717
718
719
def __init__(self,
             config: Optional[ExtinctionConfig] = None,
             **kwargs) -> None:
    if config is None:
        config = ExtinctionConfig(**kwargs)
    elif kwargs:
        # merge kwargs on top of config
        import dataclasses
        d = dataclasses.asdict(config)
        d.update(kwargs)
        config = ExtinctionConfig(**d)
    self.config = config

extinction_curve 

extinction_curve(wave)

Return A(λ) in magnitudes for the configured law and A_V.

Parameters:

Name Type Description Default
wave array_like

Wavelengths in Angstroms.

 required

Returns:

Name Type Description
A_lambda ndarray

Extinction in magnitudes. Returns zeros if enabled is False.
Source code in sed_model/sed_extinction.py 
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
def extinction_curve(self, wave: np.ndarray) -> np.ndarray:
    """Return A(λ) in magnitudes for the configured law and A_V.

    Parameters
    ----------
    wave : array_like
        Wavelengths in Angstroms.

    Returns
    -------
    A_lambda : ndarray
        Extinction in magnitudes.  Returns zeros if ``enabled`` is False.
    """
    wave = np.asarray(wave, dtype=np.float64)
    if not self.config.enabled or self.config.a_v == 0.0:
        return np.zeros(len(wave))

    law = self.config.law
    av  = self.config.a_v
    rv  = self.config.r_v

    if law == 'ccm89':
        return ccm89(wave, av, rv)
    elif law == 'odonnell94':
        return odonnell94(wave, av, rv)
    elif law == 'fitzpatrick99':
        return fitzpatrick99(wave, av, rv)
    elif law == 'fm07':
        return fm07(wave, av)
    elif law == 'calzetti00':
        return calzetti00(wave, av, rv)
    elif law == 'gordon23':
        return gordon23(wave, av, rv if rv != 3.1 else None,
                        self.config.gordon23_env)
    else:
        raise ValueError(f"Unknown law: {law}")

apply 

apply(wave, flux, r_star_cm=None)

Apply extinction (and optionally distance scaling) to a flux.

Parameters:

Name Type Description Default
wave array_like

Wavelengths in Angstroms.

 required
flux array_like

Intrinsic flux in any linear unit.

 required
r_star_cm float

Stellar radius in cm. Required only when config.scale_distance is True.

 None

Returns:

Name Type Description
flux_out ndarray

Processed flux. If enabled is False this is identical to the input flux (no copy is made if no operation is needed).
Source code in sed_model/sed_extinction.py 
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
def apply(self,
          wave: np.ndarray,
          flux: np.ndarray,
          r_star_cm: Optional[float] = None) -> np.ndarray:
    """Apply extinction (and optionally distance scaling) to a flux.

    Parameters
    ----------
    wave : array_like
        Wavelengths in Angstroms.
    flux : array_like
        Intrinsic flux in any linear unit.
    r_star_cm : float, optional
        Stellar radius in cm.  Required only when
        ``config.scale_distance`` is True.

    Returns
    -------
    flux_out : ndarray
        Processed flux.  If ``enabled`` is False this is identical to
        the input flux (no copy is made if no operation is needed).
    """
    wave = np.asarray(wave, dtype=np.float64)
    flux = np.asarray(flux, dtype=np.float64)

    if not self.config.enabled:
        return flux

    # 1. Extinction
    a_lam = self.extinction_curve(wave)
    flux_out = apply_extinction(a_lam, flux)

    # 2. Distance dilution: F_obs = F_surface × (R / d)^2
    if self.config.scale_distance:
        if r_star_cm is None:
            warnings.warn(
                "scale_distance=True but r_star_cm was not supplied; "
                "distance scaling skipped.",
                UserWarning,
                stacklevel=2,
            )
        else:
            d_cm  = self.config.distance_pc * _PC_TO_CM
            flux_out = flux_out * (r_star_cm / d_cm)**2

    return flux_out

remove 

remove(wave, flux)

Remove (de-redden) extinction from an observed flux.

Distance scaling is NOT reversed here — this only undoes extinction.

Parameters:

Name Type Description Default
wave array_like

Wavelengths in Angstroms.

 required
flux array_like

Observed flux.

 required

Returns:

Name Type Description
flux_intrinsic ndarray

De-reddened flux.
Source code in sed_model/sed_extinction.py 
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
def remove(self,
           wave: np.ndarray,
           flux: np.ndarray) -> np.ndarray:
    """Remove (de-redden) extinction from an observed flux.

    Distance scaling is NOT reversed here — this only undoes extinction.

    Parameters
    ----------
    wave : array_like
        Wavelengths in Angstroms.
    flux : array_like
        Observed flux.

    Returns
    -------
    flux_intrinsic : ndarray
        De-reddened flux.
    """
    wave = np.asarray(wave, dtype=np.float64)
    flux = np.asarray(flux, dtype=np.float64)
    if not self.config.enabled:
        return flux
    a_lam = self.extinction_curve(wave)
    return remove_extinction(a_lam, flux)

disabled classmethod 

disabled()

Return an extinction model that does nothing (convenience factory).

Source code in sed_model/sed_extinction.py 
846
847
848
849
@classmethod
def disabled(cls) -> "ExtinctionModel":
    """Return an extinction model that does nothing (convenience factory)."""
    return cls(ExtinctionConfig(enabled=False))

from_dict classmethod 

from_dict(d)

Construct from a plain dict (e.g. loaded from a config file).

Source code in sed_model/sed_extinction.py 
851
852
853
854
@classmethod
def from_dict(cls, d: dict) -> "ExtinctionModel":
    """Construct from a plain dict (e.g. loaded from a config file)."""
    return cls(ExtinctionConfig(**d))

ccm89 

ccm89(wave, a_v, r_v=3.1)

Cardelli, Clayton & Mathis (1989) extinction in magnitudes.

Parameters:

Name Type Description Default
wave array_like

Wavelengths in Angstroms. Valid range: ~910–33 000 Å.

 required
a_v float

V-band extinction in magnitudes.

 required
r_v float

Total-to-selective extinction ratio (default 3.1).

 3.1

Returns:

Name Type Description
A_lambda ndarray

Extinction in magnitudes at each wavelength.
Source code in sed_model/sed_extinction.py 
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def ccm89(wave: np.ndarray, a_v: float, r_v: float = 3.1) -> np.ndarray:
    """Cardelli, Clayton & Mathis (1989) extinction in magnitudes.

    Parameters
    ----------
    wave : array_like
        Wavelengths in Angstroms.  Valid range: ~910–33 000 Å.
    a_v : float
        V-band extinction in magnitudes.
    r_v : float
        Total-to-selective extinction ratio (default 3.1).

    Returns
    -------
    A_lambda : ndarray
        Extinction in magnitudes at each wavelength.
    """
    wave = np.asarray(wave, dtype=np.float64)
    x = _aa_to_invum(wave)
    a, b = _ccm89_ab(x)
    return a_v * (a + b / r_v)

odonnell94 

odonnell94(wave, a_v, r_v=3.1)

O'Donnell (1994) extinction — CCM89 with revised optical coefficients.

Parameters:

Name Type Description Default
wave array_like

Wavelengths in Angstroms. Valid range: ~910–33 000 Å.

 required
a_v float

V-band extinction in magnitudes.

 required
r_v float

Total-to-selective extinction ratio (default 3.1).

 3.1

Returns:

Name Type Description
A_lambda ndarray

Extinction in magnitudes at each wavelength.
Source code in sed_model/sed_extinction.py 
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def odonnell94(wave: np.ndarray, a_v: float, r_v: float = 3.1) -> np.ndarray:
    """O'Donnell (1994) extinction — CCM89 with revised optical coefficients.

    Parameters
    ----------
    wave : array_like
        Wavelengths in Angstroms.  Valid range: ~910–33 000 Å.
    a_v : float
        V-band extinction in magnitudes.
    r_v : float
        Total-to-selective extinction ratio (default 3.1).

    Returns
    -------
    A_lambda : ndarray
        Extinction in magnitudes at each wavelength.
    """
    wave = np.asarray(wave, dtype=np.float64)
    x = _aa_to_invum(wave)
    a, b = _od94_ab(x)
    return a_v * (a + b / r_v)

fitzpatrick99 

fitzpatrick99(wave, a_v, r_v=3.1)

Fitzpatrick (1999) R_V-dependent dust extinction.

Optical/IR: cubic spline through the Fitzpatrick (1999) knots updated by E. Fitzpatrick to match the IDL astrolib FM_UNRED routine. UV (< 2700 Å): analytic Fitzpatrick & Massa (1990) parametrization.

Parameters:

Name Type Description Default
wave array_like

Wavelengths in Angstroms. Valid range: 910–60 000 Å.

 required
a_v float

V-band extinction in magnitudes.

 required
r_v float

Total-to-selective extinction ratio (default 3.1).

 3.1

Returns:

Name Type Description
A_lambda ndarray

Extinction in magnitudes at each wavelength.
Source code in sed_model/sed_extinction.py 
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
def fitzpatrick99(wave: np.ndarray, a_v: float, r_v: float = 3.1) -> np.ndarray:
    """Fitzpatrick (1999) R_V-dependent dust extinction.

    Optical/IR: cubic spline through the Fitzpatrick (1999) knots updated
    by E. Fitzpatrick to match the IDL astrolib FM_UNRED routine.
    UV (< 2700 Å): analytic Fitzpatrick & Massa (1990) parametrization.

    Parameters
    ----------
    wave : array_like
        Wavelengths in Angstroms.  Valid range: 910–60 000 Å.
    a_v : float
        V-band extinction in magnitudes.
    r_v : float
        Total-to-selective extinction ratio (default 3.1).

    Returns
    -------
    A_lambda : ndarray
        Extinction in magnitudes at each wavelength.
    """
    wave  = np.asarray(wave, dtype=np.float64)
    x     = _aa_to_invum(wave)
    k_knots = _f99_kknots(r_v)

    # Optical/IR knots via spline (x <= 1e4/2700)
    opt_mask = x <= _F99_XKNOTS[7]
    uv_mask  = ~opt_mask

    result = np.empty_like(x)

    if opt_mask.any():
        k_opt = _natural_cubic_spline(_F99_XKNOTS, k_knots, x[opt_mask])
        result[opt_mask] = a_v / r_v * (k_opt + r_v)

    if uv_mask.any():
        xm  = x[uv_mask]
        c2  = -0.824 + 4.717 / r_v
        c1  =  2.030 - 3.007 * c2
        xm2 = xm**2
        d   = xm2 / ((xm2 - _F99_X0**2)**2 + xm2 * _F99_GAMMA**2)
        k   = c1 + c2*xm + _F99_C3*d
        bump = xm >= _F99_C5
        if bump.any():
            y = xm[bump] - _F99_C5
            k[bump] += _F99_C4*(0.5392*y**2 + 0.05644*y**3)
        result[uv_mask] = a_v * (1.0 + k/r_v)

    return result

fm07 

fm07(wave, a_v)

Fitzpatrick & Massa (2007) extinction (R_V = 3.1 fixed).

Defined from 910 Å to 6 µm.

Parameters:

Name Type Description Default
wave array_like

Wavelengths in Angstroms.

 required
a_v float

V-band extinction in magnitudes.

 required

Returns:

Name Type Description
A_lambda ndarray

Extinction in magnitudes at each wavelength.
Source code in sed_model/sed_extinction.py 
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
def fm07(wave: np.ndarray, a_v: float) -> np.ndarray:
    """Fitzpatrick & Massa (2007) extinction (R_V = 3.1 fixed).

    Defined from 910 Å to 6 µm.

    Parameters
    ----------
    wave : array_like
        Wavelengths in Angstroms.
    a_v : float
        V-band extinction in magnitudes.

    Returns
    -------
    A_lambda : ndarray
        Extinction in magnitudes at each wavelength.
    """
    wave  = np.asarray(wave, dtype=np.float64)
    x     = _aa_to_invum(wave)
    ebv   = a_v / _FM07_R_V

    opt_mask = x <= _FM07_XKNOTS_raw[8]   # <= 1e4/2700 Å (spline region)
    uv_mask  = ~opt_mask

    result = np.empty_like(x)

    if opt_mask.any():
        k_opt = _natural_cubic_spline(_FM07_XKNOTS_raw, _FM07_KKNOTS, x[opt_mask])
        result[opt_mask] = ebv * (k_opt + _FM07_R_V)

    if uv_mask.any():
        xm   = x[uv_mask]
        xm2  = xm**2
        d    = xm2 / ((xm2 - _FM07_X02)**2 + xm2*_FM07_G2)
        k    = _FM07_C1 + _FM07_C2*xm + _FM07_C3*d
        bump = xm > _FM07_C5
        if bump.any():
            y = xm[bump] - _FM07_C5
            k[bump] += _FM07_C4*y**2
        result[uv_mask] = a_v * (1.0 + k/_FM07_R_V)

    return result

calzetti00 

calzetti00(wave, a_v, r_v=4.05)

Calzetti et al. (2000) starburst galaxy attenuation law.

Valid range: 1200–22 000 Å.

Parameters:

Name Type Description Default
wave array_like

Wavelengths in Angstroms.

 required
a_v float

V-band attenuation in magnitudes (a_v = r_v * E(B−V)_s).

 required
r_v float

R_V for the starburst law (Calzetti default: 4.05).

 4.05

Returns:

Name Type Description
A_lambda ndarray

Attenuation in magnitudes at each wavelength.
Source code in sed_model/sed_extinction.py 
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
def calzetti00(wave: np.ndarray, a_v: float, r_v: float = 4.05) -> np.ndarray:
    """Calzetti et al. (2000) starburst galaxy attenuation law.

    Valid range: 1200–22 000 Å.

    Parameters
    ----------
    wave : array_like
        Wavelengths in Angstroms.
    a_v : float
        V-band attenuation in magnitudes (a_v = r_v * E(B−V)_s).
    r_v : float
        R_V for the starburst law (Calzetti default: 4.05).

    Returns
    -------
    A_lambda : ndarray
        Attenuation in magnitudes at each wavelength.
    """
    wave  = np.asarray(wave, dtype=np.float64)
    w_um  = wave * 1e-4   # Angstroms → microns

    k = np.zeros_like(w_um)

    # UV-optical: 0.12–0.63 µm
    m = (w_um >= 0.12) & (w_um < 0.63)
    if m.any():
        w = w_um[m]
        k[m] = (2.659*(-2.156 + 1.509/w - 0.198/w**2 + 0.011/w**3) + r_v)

    # Optical-NIR: 0.63–2.2 µm
    m = (w_um >= 0.63) & (w_um <= 2.2)
    if m.any():
        w = w_um[m]
        k[m] = 2.659*(-1.857 + 1.040/w) + r_v

    ebv_s = a_v / r_v
    return ebv_s * k

gordon23 

gordon23(wave, a_v, r_v=None, environment='mw')

Gordon et al. (2023) dust extinction law.

Native implementation of the piecewise FM-style law from Gordon, K. D. et al. 2023, ApJ, 950, 86.

Parameters:

Name Type Description Default
wave array_like

Wavelengths in Angstroms.

 required
a_v float

V-band extinction in magnitudes.

 required
r_v float or None

Total-to-selective extinction ratio. If None, the environment default from Gordon+2023 is used (MW: 3.17, LMC: 3.41, SMC: 2.74).

 None
environment ('mw', 'lmc', 'smc')

Which average extinction curve to use (default 'mw').

 'mw'

Returns:

Name Type Description
A_lambda ndarray

Extinction in magnitudes at each wavelength.
Source code in sed_model/sed_extinction.py 
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
def gordon23(wave: np.ndarray, a_v: float, r_v: Optional[float] = None,
             environment: str = 'mw') -> np.ndarray:
    """Gordon et al. (2023) dust extinction law.

    Native implementation of the piecewise FM-style law from
    Gordon, K. D. et al. 2023, ApJ, 950, 86.

    Parameters
    ----------
    wave : array_like
        Wavelengths in Angstroms.
    a_v : float
        V-band extinction in magnitudes.
    r_v : float or None
        Total-to-selective extinction ratio.  If None, the environment
        default from Gordon+2023 is used (MW: 3.17, LMC: 3.41, SMC: 2.74).
    environment : {'mw', 'lmc', 'smc'}
        Which average extinction curve to use (default 'mw').

    Returns
    -------
    A_lambda : ndarray
        Extinction in magnitudes at each wavelength.
    """
    env = environment.lower()
    if env not in _G23_PARAMS:
        raise ValueError(f"environment must be one of {list(_G23_PARAMS)}; got '{env}'")

    p   = _G23_PARAMS[env]
    c1, c2, c3, c4, c5 = p['c1'], p['c2'], p['c3'], p['c4'], p['c5']
    x0, gamma           = p['x0'], p['gamma']
    c6, alpha           = p['c6'], p['alpha']
    rv  = r_v if r_v is not None else p['r_v']

    wave = np.asarray(wave, dtype=np.float64)
    x    = _aa_to_invum(wave)    # inverse microns
    k    = np.zeros_like(x)

    # --- IR: x < 1/3 µm⁻¹  (> 3 µm) — power law ---
    m = x < (1.0/3.0)
    if m.any():
        k[m] = c6 * x[m]**alpha - rv

    # --- Optical+UV: x >= 1/3 µm⁻¹ ---
    m = ~m
    if m.any():
        xm   = x[m]
        xm2  = xm**2
        g2   = gamma**2
        x02  = x0**2
        d    = xm2 / ((xm2 - x02)**2 + xm2*g2)
        k[m] = c1 + c2*xm + c3*d

        # Far-UV non-linear term
        fuv  = xm > c5
        if fuv.any():
            y = xm[fuv] - c5
            k[m][fuv] += c4*(0.5392*y**2 + 0.05644*y**3)

    # Convert k → A(λ) = a_v/r_v * (k + r_v)
    return a_v / rv * (k + rv)

apply_extinction 

apply_extinction(a_lambda, flux)

Apply extinction magnitudes to a flux array.

F_obs = F_intrinsic × 10^(−0.4 × A_lambda)

Parameters:

Name Type Description Default
a_lambda ndarray

Extinction in magnitudes at each wavelength (same shape as flux).

 required
flux ndarray

Intrinsic flux array.

 required

Returns:

Name Type Description
flux_obs ndarray

Observed (extincted) flux.
Source code in sed_model/sed_extinction.py 
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
def apply_extinction(a_lambda: np.ndarray, flux: np.ndarray) -> np.ndarray:
    """Apply extinction magnitudes to a flux array.

    F_obs = F_intrinsic × 10^(−0.4 × A_lambda)

    Parameters
    ----------
    a_lambda : ndarray
        Extinction in magnitudes at each wavelength (same shape as flux).
    flux : ndarray
        Intrinsic flux array.

    Returns
    -------
    flux_obs : ndarray
        Observed (extincted) flux.
    """
    return flux * 10.0**(-0.4 * a_lambda)

remove_extinction 

remove_extinction(a_lambda, flux)

Remove extinction from an observed flux (de-redden).

F_intrinsic = F_obs × 10^(+0.4 × A_lambda)

Parameters:

Name Type Description Default
a_lambda ndarray

Extinction in magnitudes at each wavelength.

 required
flux ndarray

Observed (extincted) flux.

 required

Returns:

Name Type Description
flux_intrinsic ndarray

De-reddened flux.
Source code in sed_model/sed_extinction.py 
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
def remove_extinction(a_lambda: np.ndarray, flux: np.ndarray) -> np.ndarray:
    """Remove extinction from an observed flux (de-redden).

    F_intrinsic = F_obs × 10^(+0.4 × A_lambda)

    Parameters
    ----------
    a_lambda : ndarray
        Extinction in magnitudes at each wavelength.
    flux : ndarray
        Observed (extincted) flux.

    Returns
    -------
    flux_intrinsic : ndarray
        De-reddened flux.
    """
    return flux * 10.0**(0.4 * a_lambda)

make_extinction_model 

make_extinction_model(enabled=False, law='fitzpatrick99', r_v=3.1, a_v=0.0, distance_pc=1.0, scale_distance=False, gordon23_env='mw')

Convenience factory — mirrors the ExtinctionConfig keyword signature.

This is the recommended way to build an extinction model inside a fitter's forward-model function, e.g.::

ext = make_extinction_model(enabled=True, law='fitzpatrick99',
                            a_v=theta['a_v'])
flux_obs = ext.apply(wave, flux_intrinsic)

Parameters:

Name Type Description Default
enabled bool

Master switch. Default False — extinction is skipped entirely.

 False
law str

Extinction law name. See AVAILABLE_LAWS.

 'fitzpatrick99'
r_v float

R_V (shape parameter). Default 3.1.

 3.1
a_v float

V-band extinction. Default 0.0.

 0.0
distance_pc float

Source distance in parsecs. Default 1 pc (absolute flux, no dilution).

 1.0
scale_distance bool

Apply (R/d)^2 dilution. Default False.

 False
gordon23_env str

Gordon+2023 environment preset: 'mw', 'lmc', or 'smc'.

 'mw'

Returns:

Type Description
ExtinctionModel
Source code in sed_model/sed_extinction.py 
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
def make_extinction_model(
    enabled: bool = False,
    law: str = 'fitzpatrick99',
    r_v: float = 3.1,
    a_v: float = 0.0,
    distance_pc: float = 1.0,
    scale_distance: bool = False,
    gordon23_env: str = 'mw',
) -> ExtinctionModel:
    """Convenience factory — mirrors the ExtinctionConfig keyword signature.

    This is the recommended way to build an extinction model inside a fitter's
    forward-model function, e.g.::

        ext = make_extinction_model(enabled=True, law='fitzpatrick99',
                                    a_v=theta['a_v'])
        flux_obs = ext.apply(wave, flux_intrinsic)

    Parameters
    ----------
    enabled : bool
        Master switch.  Default False — extinction is skipped entirely.
    law : str
        Extinction law name.  See ``AVAILABLE_LAWS``.
    r_v : float
        R_V (shape parameter).  Default 3.1.
    a_v : float
        V-band extinction.  Default 0.0.
    distance_pc : float
        Source distance in parsecs.  Default 1 pc (absolute flux, no dilution).
    scale_distance : bool
        Apply (R/d)^2 dilution.  Default False.
    gordon23_env : str
        Gordon+2023 environment preset: 'mw', 'lmc', or 'smc'.

    Returns
    -------
    ExtinctionModel
    """
    return ExtinctionModel(ExtinctionConfig(
        enabled=enabled,
        law=law,
        r_v=r_v,
        a_v=a_v,
        distance_pc=distance_pc,
        scale_distance=scale_distance,
        gordon23_env=gordon23_env,
    ))