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, flux_fortran, model_dir, lookup_nodes=None)

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.
flux_fortran (ndarray, shape(n_T, n_L, n_M, n_W))

Fortran-contiguous copy of flux for zero-copy Fortran calls. Pre-computed at load time so the hot path avoids repeated copies.
model_dir Path

Directory from which the grid was loaded.
lookup_nodes (ndarray, shape(N, 3), optional)

(Teff, logg, [M/H]) of the tabulated models, i.e. the cube nodes that hold a spectrum. MESA's Interp_rad is measured to the entries of lookup_table.csv, which lists exactly these models (a rectangular cube also has empty nodes where no model exists). Computed from flux when not supplied.

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 
88
89
90
91
92
93
94
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 
 96
 97
 98
 99
100
101
102
103
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

cell_has_empty_nodes 

cell_has_empty_nodes(teff, logg, meta, stencil=False)

True if the cell containing the query touches a node with no model.

The cell is the one MESA's interpolators use (out-of-range queries take the boundary cell, dummy and single-point axes collapse). A corner counts only if it carries interpolation weight, so a query sitting exactly on a populated node is never reported. Corners that hold no spectrum are all-zero in the cube, and interpolating through them gives unreliable fluxes.

With stencil=True the check also covers the nodes the Hermite derivative estimates read: each corner's axis-aligned neighbours, i.e. one extra node on either side of the cell along each axis. Those feed the cubic and an empty one makes it overshoot even when all eight corners are populated. Linear interpolation does not use them.

Source code in sed_model/grid.py 
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
def cell_has_empty_nodes(
    self, teff: float, logg: float, meta: float, stencil: bool = False
) -> bool:
    """True if the cell containing the query touches a node with no model.

    The cell is the one MESA's interpolators use (out-of-range queries
    take the boundary cell, dummy and single-point axes collapse).  A
    corner counts only if it carries interpolation weight, so a query
    sitting exactly on a populated node is never reported.  Corners that
    hold no spectrum are all-zero in the cube, and interpolating through
    them gives unreliable fluxes.

    With ``stencil=True`` the check also covers the nodes the Hermite
    derivative estimates read: each corner's axis-aligned neighbours, i.e.
    one extra node on either side of the cell along each axis.  Those
    feed the cubic and an empty one makes it overshoot even when all
    eight corners are populated.  Linear interpolation does not use them.
    """
    sel, live = [], []
    for axis, val in ((self.teff_grid, teff), (self.logg_grid, logg),
                      (self.meta_grid, meta)):
        n = len(axis)
        if n < 2 or np.all(axis == 0.0) or np.all(axis == 999.0) \
                or np.all(axis == -999.0):
            sel.append([0])
            live.append(0)
            continue
        live.append(n)
        if val <= axis[0]:
            sel.append([0])
        elif val >= axis[-1]:
            sel.append([n - 1])
        else:
            i = int(np.searchsorted(axis, val, side="right")) - 1
            sel.append([i] if axis[i] == val else [i, i + 1])
    has = self._node_has_model
    if not has[np.ix_(*sel)].all():
        return True
    if not stencil:
        return False
    for ax, idx in enumerate(sel):
        for nb in (idx[0] - 1, idx[-1] + 1):
            if live[ax] and 0 <= nb < live[ax]:
                ext = list(sel)
                ext[ax] = [nb]
                if not has[np.ix_(*ext)].all():
                    return True
    return False

nearest_grid_distance 

nearest_grid_distance(teff, logg, meta)

Distance in normalised parameter space from the nearest tabulated model; identical to the Interp_rad diagnostic of the MESA colors module (compute_interp_radius in colors/private/bolometric.f90).

Each axis is normalised by the range spanned by the tabulated models (floored at 1e-12), and an axis on which every model has the value 0, 999 or -999 (i.e. the axis does not exist for this grid) is ignored. The query is not clamped to the grid.

Source code in sed_model/grid.py 
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
def nearest_grid_distance(
    self, teff: float, logg: float, meta: float
) -> float:
    """Distance in normalised parameter space from the nearest tabulated
    model; identical to the ``Interp_rad`` diagnostic of the MESA colors
    module (``compute_interp_radius`` in colors/private/bolometric.f90).

    Each axis is normalised by the range spanned by the tabulated models
    (floored at 1e-12), and an axis on which every model has the value 0,
    999 or -999 (i.e. the axis does not exist for this grid) is ignored.
    The query is not clamped to the grid.
    """
    nodes = self._tabulated_nodes()
    q = (teff, logg, meta)
    total = np.zeros(len(nodes))
    for k in range(3):
        col = nodes[:, k]
        if np.all(col == 0.0) or np.all(col == 999.0) or np.all(col == -999.0):
            continue
        lo = col.min()
        rng = max(col.max() - lo, 1.0e-12)
        total += ((q[k] - lo) / rng - (col - lo) / rng) ** 2
    return float(np.sqrt(total.min()))

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 
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
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,
        flux_fortran=np.asfortranarray(flux, dtype=np.float64),
        model_dir=model_dir,
        lookup_nodes=_tabulated_nodes(teff_grid, logg_grid, meta_grid, flux),
    )

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.

Returns:

Type Description
DataFrame
Source code in sed_model/grid.py 
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
def validate_lookup_table(grid: AtmosphereGrid):
    """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.

    Returns
    -------
    pandas.DataFrame
    """
    import pandas as pd

    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 integrals as the MESA colors module (colors/private/synthetic.f90): photon counting (weight by wavelength), as MESA does.

Detector type

SED_Tools tags each filter file with the response convention of its curve, DetectorType=photon or energy, as a ;-separated field after the second comma of the column-name row::

Wavelength,Transmission,FilterID=GAIA/GAIA3.G;DetectorType=photon;...

Only the column row is read. An energy-counting curve T_E is converted to its photon-counting equivalent T_E/λ (Bessell & Murphy 2012, eq. A9), rescaled to keep its peak, so that the λ-weighted band integrals below equal ∫F T_E / ∫T_E. This is what the MESA colors module does. A missing or unrecognised tag means photon counting. A WavelengthUnit tag other than Angstrom converts the file's wavelength column to Angstrom.

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)

Zero-points are computed through the compiled cc_api Fortran extension so the Python loader uses the same numerical kernels as run_forward() and the current MESA Colors implementation. A pure-Python trapezoid fallback is kept for environments where the extension has not yet been built.

Filter dataclass 

Filter(name, path, wavelengths, transmission, vega_zero_point=-1.0, ab_zero_point=-1.0, st_zero_point=-1.0, photon_counting=True, detector_type='photon')

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.
photon_counting bool

True when the band integrals are weighted by wavelength. Always True for loaded filters: energy-counting curves are converted to their photon-counting equivalent when loaded (see detector_type).
detector_type str

'photon' or 'energy': the response convention the file's curve was tagged with. For 'energy' the transmission array holds the converted curve T_E/λ (rescaled to the original peak), not the file's values.

zero_point 

zero_point(system)

Return the zero-point for system.

Parameters:

Name Type Description Default
system str

One of 'Vega', 'AB', 'ST' (case-insensitive).

 required
Source code in sed_model/filters.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def zero_point(self, system: str) -> float:
    """Return the zero-point for *system*.

    Parameters
    ----------
    system : str
        One of ``'Vega'``, ``'AB'``, ``'ST'`` (case-insensitive).
    """
    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. A DetectorType tag (see the module docstring) selects photon or energy convention.

 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 
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
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.  A ``DetectorType`` tag
        (see the module docstring) selects photon or energy convention.
    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)
        header = _read_filter_header(p)
        wave = wave * _wavelength_factor(header, p)
        detector = _detector_type(header, p)
        if detector == "energy":
            trans = _energy_to_photon_response(wave, trans)
        pc = True
        filt = Filter(
            name=p.stem,
            path=p,
            wavelengths=wave,
            transmission=trans,
            photon_counting=pc,
            detector_type=detector,
        )
        filt.ab_zero_point = _compute_ab_zero_point(wave, trans, pc)
        filt.st_zero_point = _compute_st_zero_point(wave, trans, pc)
        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, pc
            )
        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.

Filters are returned in index-file order, or alphabetically when the index file is absent. Capture [f.name for f in filters] immediately after loading to record the stable order for use with obs_magnitudes arrays.

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 
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
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.

    Filters are returned in index-file order, or alphabetically when
    the index file is absent.  Capture ``[f.name for f in filters]``
    immediately after loading to record the stable order for use with
    ``obs_magnitudes`` arrays.

    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():
        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, bounded Hermite, 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.

Calling conventions

Classic (keyword-explicit)::

run_forward(teff=5778, logg=4.44, meta=0.0,
            R=6.957e10, d=3.086e19,
            grid=grid, filters=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.

ForwardResult dataclass 

ForwardResult(wavelengths, surface_flux, observed_flux, magnitudes, band_fluxes, bol_flux, bol_mag, nearest_grid_distance, 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.

Attributes:

Name Type Description
nearest_grid_distance float

Euclidean distance in normalised parameter space from the nearest grid point (mirrors the MESA Interp_rad diagnostic). Zero for an exact grid point, larger as the query moves away from the grid.
magnitudes dict[str, float]

Synthetic magnitudes keyed by filter name, in the order the filters list was provided to run_forward.

run_forward 

run_forward(teff=None, logg=None, meta=None, R=None, d=None, grid=None, filters=None, mag_system='AB', interp_method='hermite_bounded', 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')

Case-insensitive.

 'AB'
interp_method ('hermite_bounded', 'hermite', 'linear')

'hermite_bounded' (default, as in MESA colors) uses Hermite interpolation, clips negligible negative undershoot to zero, and falls back to trilinear interpolation where Hermite is unusable. 'hermite' is the raw Hermite interpolant; 'linear' is positivity-preserving trilinear interpolation.

 'hermite_bounded'
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 
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
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
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_bounded",
    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'}
        Case-insensitive.
    interp_method : {'hermite_bounded', 'hermite', 'linear'}
        ``'hermite_bounded'`` (default, as in MESA colors) uses Hermite
        interpolation, clips negligible negative undershoot to zero, and
        falls back to trilinear interpolation where Hermite is unusable.
        ``'hermite'`` is the raw Hermite interpolant; ``'linear'`` is
        positivity-preserving trilinear interpolation.
    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.
    """
    if R is None:
        raise ValueError("R (stellar radius in cm) is always required.")

    cc = get_cc_api(required=True)

    # ------------------------------------------------------------------
    # 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, d, grid, filters)):
            raise ValueError(
                "teff, logg, meta, d, grid, and filters must all be supplied "
                "when fit_params is not used."
            )
        d_use = float(d)
        av    = 0.0

    # ------------------------------------------------------------------
    # 1. SED interpolation (Fortran)
    # ------------------------------------------------------------------
    clamped  = not grid.in_bounds(teff, logg, meta)
    grid_dist = grid.nearest_grid_distance(teff, logg, meta)
    teff_q, logg_q, meta_q = grid.clamp(teff, logg, meta)

    method = interp_method.lower()
    hermite = method in ("hermite", "hermite_bounded")
    if grid.cell_has_empty_nodes(teff_q, logg_q, meta_q, stencil=hermite):
        warnings.warn(
            f"Teff={teff_q:.1f}, logg={logg_q:.3f}, [M/H]={meta_q:.3f} lies in a "
            "grid cell with empty nodes (no model at a corner"
            + (" or in the Hermite stencil around it" if hermite else "")
            + "); the interpolated SED and magnitudes are unreliable there.",
            UserWarning, stacklevel=2,
        )

    # grid.flux_fortran is pre-computed at load time — zero-copy to Fortran.
    # The axis arrays are already float64 C-contiguous from load_grid.
    if method == "hermite_bounded":
        surface_flux, ierr = cc.interp_sed_hermite_bounded(
            teff_q, logg_q, meta_q,
            grid.teff_grid, grid.logg_grid, grid.meta_grid,
            grid.flux_fortran, grid.wavelengths,
        )
    elif method == "hermite":
        surface_flux, ierr = cc.interp_sed_hermite(
            teff_q, logg_q, meta_q,
            grid.teff_grid, grid.logg_grid, grid.meta_grid,
            grid.flux_fortran,
        )
    elif method == "linear":
        surface_flux, ierr = cc.interp_sed_linear(
            teff_q, logg_q, meta_q,
            grid.teff_grid, grid.logg_grid, grid.meta_grid,
            grid.flux_fortran,
        )
    else:
        raise ValueError(
            f"Unknown interp_method '{interp_method}'. "
            "Choose 'hermite_bounded', 'hermite' or 'linear'."
        )

    if method == "hermite_bounded" and ierr == 2:
        raise ValueError(
            "Hermite and linear interpolation both returned an invalid SED "
            f"at Teff={teff:.1f}, logg={logg:.3f}, [M/H]={meta:.3f}."
        )
    # Hermite ierr=1 is a boundary clamp.  Linear ierr=1 only records that some
    # wavelength bins fell back to the nearest node, which says nothing about
    # the query; only ierr=2 (no usable flux at all) is flagged.
    if (method == "hermite" and ierr != 0) or (method == "linear" and ierr == 2):
        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.
    # 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:
            extinction = extinction.with_av(av)

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

    # ------------------------------------------------------------------
    # 4. Bolometric quantities (Fortran)
    # ------------------------------------------------------------------
    bol_flux, bol_mag, bol_ierr = cc.bolometric(grid.wavelengths, observed_flux)
    if bol_ierr != 0:
        clamped = True

    # ------------------------------------------------------------------
    # 5. Synthetic photometry per filter (Fortran)
    #
    # Filter arrays are float64 C-contiguous from load_filters — no copy needed.
    # ------------------------------------------------------------------
    magnitudes:  dict = {}
    band_fluxes: dict = {}

    for filt in filters:
        zp = filt.zero_point(mag_system)
        mag, band_flux, mag_ierr = cc.synthetic_magnitude(
            grid.wavelengths, observed_flux,
            filt.wavelengths, filt.transmission,
            zp,
            photon_counting=int(filt.photon_counting),
        )
        magnitudes[filt.name]  = float(mag) if mag_ierr == 0 else float('nan')
        band_fluxes[filt.name] = float(band_flux)

    return ForwardResult(
        wavelengths=grid.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),
        nearest_grid_distance=float(grid_dist),
        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_bounded', 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 
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
def run_forward_batch(
    params: np.ndarray,
    R: float,
    d: float,
    grid: AtmosphereGrid,
    filters: list,
    mag_system: str = "AB",
    interp_method: str = "hermite_bounded",
    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. Two 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_bounded', n_walkers=32, n_steps=2000, n_burn=500, n_thin=1, p0_centre=None, p0_scatter=0.02, seed=None, progress=True, verbose=True, init='auto', sys_floor=0.0, fit_jitter=False, jitter_max=0.2, priors=None, upper_limits=None, auto_stop=False, moves=None, pool=None)

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. Non-finite entries (NaN/inf) are treated as missing and dropped, with a warning.

 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_law str or None

Convenience shorthand used to build a default extinction model when both fit_params and extinction are None and Av is active. Ignored (with a warning) when fit_params is supplied.

 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_bounded', 'hermite', 'linear')
'hermite_bounded'
n_walkers int

emcee settings. n_steps is the (maximum) chain length; n_burn=None chooses the burn-in automatically (2 integrated autocorrelation times, at most half the chain).

 32
n_steps int

emcee settings. n_steps is the (maximum) chain length; n_burn=None chooses the burn-in automatically (2 integrated autocorrelation times, at most half the chain).

 32
n_burn int

emcee settings. n_steps is the (maximum) chain length; n_burn=None chooses the burn-in automatically (2 integrated autocorrelation times, at most half the chain).

 32
n_thin int

emcee settings. n_steps is the (maximum) chain length; n_burn=None chooses the burn-in automatically (2 integrated autocorrelation times, at most half the chain).

 32
p0_centre dict or None

Optional starting point keyed by parameter name, e.g. {'teff': 5800, 'logg': 4.4, 'meta': 0.0}. When given, walkers start in a small Gaussian ball around it (the legacy behaviour, equivalent to init='ball').

 None
p0_scatter float

Width of initial ball as fraction of each parameter's range (only used by init='ball' / p0_centre).

 0.02
seed int or None
None
progress bool

Show the emcee progress bar.

 True
verbose bool

Print the parameter summary and extinction configuration before sampling. Set to False for batch runs.

 True
init ('auto', 'ball')

Walker initialisation. 'auto' (default) does a stratified global search over the prior box (log-uniform in Teff and distance), refines the best distinct candidates with Levenberg-Marquardt, and starts the walkers in Laplace-approximation ellipsoids around the optima (several modes are populated in proportion to their approximate mass). With p0_centre supplied, or init='ball', the legacy Gaussian ball (centre = p0_centre or the middle of the prior box) is used.

 'auto'
sys_floor float

Systematic error floor in mag, added in quadrature to every uncertainty (default 0 = off).

 0.0
fit_jitter bool

Fit one extra common noise term sigma_sys (mag, flat prior on [0, jitter_max]) added in quadrature to all uncertainties. It appears as the last column of samples / entry of param_names.

 False
priors dict or None

Extra priors on free parameters, e.g. {'d': ('parallax', 2.0, 0.1)} (mas) or {'a_v': ('halfnormal', 0.5)}. See sed_model._inv_core.parse_priors for the supported forms; a callable f(value) -> log prior is also accepted.

 None
upper_limits array-like of bool or None

Mark bands whose obs_magnitudes entry is a magnitude lower limit (flux upper limit). Their likelihood is Phi((m_model - m_obs)/sigma).

 None
auto_stop bool

Sample in blocks and stop early once the chain is >50 autocorrelation times long and the estimate has stabilised; warns if n_steps is reached first.

 False
moves emcee move, list, str or None

Passed to emcee.EnsembleSampler. None (default) means the differential-evolution mixture 'de' (mode-hopping jumps; handles multimodal/degenerate posteriors) when init='auto', and emcee's stretch move for the legacy ball. 'stretch' forces the latter.

 None
pool multiprocessing-style pool or None

Evaluate the posterior in parallel.

 None

Returns:

Type Description
InverseResult

Includes converged, n_eff, chi2_best, n_data and init_info diagnostics in addition to the chain.
Source code in sed_model/inverse.py 
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
375
376
377
378
379
380
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
426
427
428
429
430
431
432
433
434
435
436
437
438
439
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
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
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
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_bounded",
    n_walkers:         int   = 32,
    n_steps:           int   = 2000,
    n_burn:            Optional[int] = 500,
    n_thin:            int   = 1,
    p0_centre:         Optional[dict] = None,
    p0_scatter:        float = 0.02,
    seed:              Optional[int]  = None,
    progress:          bool  = True,
    verbose:           bool  = True,
    # --- options added after 0.3.1 (all default to a sensible/legacy-safe value) ---
    init:              str   = "auto",
    sys_floor:         float = 0.0,
    fit_jitter:        bool  = False,
    jitter_max:        float = 0.2,
    priors:            Optional[dict] = None,
    upper_limits:      Optional[list] = None,
    auto_stop:         bool  = False,
    moves=None,
    pool=None,
) -> 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``.  Non-finite
        entries (NaN/inf) are treated as missing and dropped, with a warning.
    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_law : str or None
        Convenience shorthand used to build a default extinction model when
        both ``fit_params`` and ``extinction`` are None and Av is active.
        Ignored (with a warning) when ``fit_params`` is supplied.
    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_bounded', 'hermite', 'linear'}
    n_walkers, n_steps, n_burn, n_thin : int
        emcee settings.  ``n_steps`` is the (maximum) chain length;
        ``n_burn=None`` chooses the burn-in automatically (2 integrated
        autocorrelation times, at most half the chain).
    p0_centre : dict or None
        Optional starting point keyed by parameter name, e.g.
        ``{'teff': 5800, 'logg': 4.4, 'meta': 0.0}``.  When given, walkers
        start in a small Gaussian ball around it (the legacy behaviour,
        equivalent to ``init='ball'``).
    p0_scatter : float
        Width of initial ball as fraction of each parameter's range
        (only used by ``init='ball'`` / ``p0_centre``).
    seed : int or None
    progress : bool
        Show the emcee progress bar.
    verbose : bool
        Print the parameter summary and extinction configuration before
        sampling.  Set to False for batch runs.
    init : {'auto', 'ball'}
        Walker initialisation.  ``'auto'`` (default) does a stratified global
        search over the prior box (log-uniform in Teff and distance), refines
        the best distinct candidates with Levenberg-Marquardt, and starts the
        walkers in Laplace-approximation ellipsoids around the optima
        (several modes are populated in proportion to their approximate mass).
        With ``p0_centre`` supplied, or ``init='ball'``, the legacy Gaussian
        ball (centre = p0_centre or the middle of the prior box) is used.
    sys_floor : float
        Systematic error floor in mag, added in quadrature to every
        uncertainty (default 0 = off).
    fit_jitter : bool
        Fit one extra common noise term ``sigma_sys`` (mag, flat prior on
        [0, ``jitter_max``]) added in quadrature to all uncertainties.  It
        appears as the last column of ``samples`` / entry of ``param_names``.
    priors : dict or None
        Extra priors on free parameters, e.g.
        ``{'d': ('parallax', 2.0, 0.1)}`` (mas) or ``{'a_v': ('halfnormal', 0.5)}``.
        See ``sed_model._inv_core.parse_priors`` for the supported forms; a
        callable ``f(value) -> log prior`` is also accepted.
    upper_limits : array-like of bool or None
        Mark bands whose ``obs_magnitudes`` entry is a magnitude *lower limit*
        (flux upper limit).  Their likelihood is Phi((m_model - m_obs)/sigma).
    auto_stop : bool
        Sample in blocks and stop early once the chain is >50 autocorrelation
        times long and the estimate has stabilised; warns if ``n_steps`` is
        reached first.
    moves : emcee move, list, str or None
        Passed to ``emcee.EnsembleSampler``.  ``None`` (default) means the
        differential-evolution mixture ``'de'`` (mode-hopping jumps; handles
        multimodal/degenerate posteriors) when ``init='auto'``, and emcee's
        stretch move for the legacy ball.  ``'stretch'`` forces the latter.
    pool : multiprocessing-style pool or None
        Evaluate the posterior in parallel.

    Returns
    -------
    InverseResult
        Includes ``converged``, ``n_eff``, ``chi2_best``, ``n_data`` and
        ``init_info`` diagnostics in addition to the chain.
    """
    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 is not None and n_burn >= n_steps:
        raise ValueError("n_burn must be < n_steps.")
    if init not in ("auto", "ball"):
        raise ValueError("init must be 'auto' or 'ball'.")
    if sys_floor < 0:
        raise ValueError("sys_floor must be >= 0.")

    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}"
        )

    upper = (np.zeros(obs_mag.size, bool) if upper_limits is None
             else np.asarray(upper_limits, dtype=bool))
    if upper.shape != obs_mag.shape:
        raise ValueError("upper_limits must have the same length as obs_magnitudes.")

    finite = np.isfinite(obs_mag) & np.isfinite(obs_err)
    if not finite.all():
        warnings.warn(
            f"Dropping {int((~finite).sum())} band(s) with non-finite magnitude or "
            f"uncertainty: {[n for n, ok in zip(filter_names, finite) if not ok]}",
            UserWarning, stacklevel=2,
        )
        obs_mag, obs_err, upper = obs_mag[finite], obs_err[finite], upper[finite]
        filter_names = [n for n, ok in zip(filter_names, finite) if ok]
    if obs_mag.size == 0 or upper.all():
        raise ValueError("No detected (finite, non-upper-limit) bands to fit.")

    # ------------------------------------------------------------------
    # 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)
    elif extinction_law is not None:
        warnings.warn(
            "extinction_law is ignored when fit_params is provided. "
            "Configure the extinction law via the extinction parameter instead.",
            UserWarning,
            stacklevel=2,
        )

    n_dim_phys = fit_params.n_free
    if n_dim_phys == 0:
        raise ValueError(
            "All parameters are fixed — nothing to sample. "
            "Make at least one parameter free in fit_params."
        )
    n_dim = n_dim_phys + (1 if fit_jitter else 0)
    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 none provided.
    # ------------------------------------------------------------------
    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,
        )

    # ------------------------------------------------------------------
    # Optionally log what is being sampled vs fixed
    # ------------------------------------------------------------------
    if verbose:
        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")

    # ------------------------------------------------------------------
    # Posterior object (lean photometry evaluator; same kernels as run_forward)
    # ------------------------------------------------------------------
    model = _FastModel(grid, filters, list(filter_names), R, mag_system,
                       interp_method, extinction)
    problem = _Problem(model, fit_params, obs_mag, obs_err, upper=upper,
                       priors=priors, sys_floor=sys_floor,
                       fit_jitter=fit_jitter, jitter_max=jitter_max)

    # ------------------------------------------------------------------
    # Initial walker positions
    # ------------------------------------------------------------------
    rng = np.random.default_rng(seed)
    use_ball = (init == "ball") or (p0_centre is not None)
    if use_ball:
        p0 = fit_params.initial_ball(
            n_walkers, centre=p0_centre or None, scatter=p0_scatter, rng=rng
        )
        if fit_jitter:
            extra = np.clip(np.abs(rng.normal(0.0, 0.1 * jitter_max, n_walkers)) + 1e-4,
                            1e-4, jitter_max)
            p0 = np.column_stack([p0, extra])
        init_info = dict(method="ball")
    else:
        p0, init_info = initial_walkers(problem, n_walkers, rng)

    # ------------------------------------------------------------------
    # Run MCMC
    # ------------------------------------------------------------------
    if moves is None and not use_ball:
        # Differential-evolution moves mix between modes and decorrelate ~2x
        # faster than the plain stretch move (see CHANGELOG); the legacy ball
        # path keeps emcee's default stretch move.
        moves = "de"
    if isinstance(moves, str):
        moves = _named_moves(emcee, moves)
    if pool is None:
        sampler = emcee.EnsembleSampler(n_walkers, n_dim, problem.batch,
                                        moves=moves, vectorize=True)
    else:
        sampler = emcee.EnsembleSampler(n_walkers, n_dim, problem,
                                        moves=moves, pool=pool)

    done = 0
    state = p0
    converged_auto = False
    if auto_stop:
        block = max(100, n_steps // 20)
        last_tau = None
        while done < n_steps:
            k = min(block, n_steps - done)
            state = sampler.run_mcmc(state, k, progress=progress, skip_initial_state_check=done > 0)
            done += k
            try:
                tau_now = sampler.get_autocorr_time(tol=0, quiet=True)
            except Exception:
                continue
            tmax = float(np.nanmax(tau_now))
            if np.isfinite(tmax) and done > 50.0 * tmax:
                if last_tau is not None and abs(tmax - last_tau) / tmax < 0.05:
                    converged_auto = True
                    break
            if np.isfinite(tmax):
                last_tau = tmax
        if not converged_auto:
            warnings.warn(
                f"auto_stop: chain did not reach 50 autocorrelation times within "
                f"n_steps={n_steps}; increase n_steps or check for multimodality.",
                UserWarning, stacklevel=2,
            )
    else:
        sampler.run_mcmc(state, n_steps, progress=progress)
        done = n_steps

    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(tol=0, quiet=True)
        if not np.all(np.isfinite(tau)):
            tau = None
    except Exception:
        tau = None

    if n_burn is None:
        n_burn_eff = 0 if tau is None else int(min(done // 2, np.ceil(2.0 * np.nanmax(tau))))
    else:
        n_burn_eff = int(min(n_burn, done - 1))

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

    converged = None
    n_eff = None
    if tau is not None:
        n_post = done - n_burn_eff
        converged = bool(n_post > 50.0 * float(np.nanmax(tau)))
        n_eff = n_walkers * n_post / np.maximum(tau, 1.0)

    ib = int(np.argmax(flat_log_prob))
    chi2_best = float(problem.chi2(flat_samples[ib]))

    if converged is False and not (auto_stop and not converged_auto):
        warnings.warn(
            "The chain is not converged (shorter than 50 autocorrelation times "
            "after burn-in), so the posterior may be unreliable. Increase "
            "n_steps, or set auto_stop=True to sample until it is.",
            UserWarning, stacklevel=2,
        )
    if (not fit_jitter and problem.n_data > 0 and np.isfinite(chi2_best)
            and chi2_best / problem.n_data > _POOR_FIT_CHI2_PER_DATUM):
        warnings.warn(
            f"Poor fit: chi2_best/n_data = {chi2_best / problem.n_data:.1f} "
            f"(> {_POOR_FIT_CHI2_PER_DATUM:g}) with fit_jitter=False. The quoted "
            "uncertainties probably under-estimate the model/photometry "
            "mismatch; consider fit_jitter=True (or sys_floor).",
            UserWarning, stacklevel=2,
        )
    init_info = dict(init_info)
    init_info["n_evals"] = int(problem.n_evals)

    fixed_params = {name: getattr(fit_params, 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=done,
        n_burn=n_burn_eff,
        n_thin=n_thin,
        acceptance_fraction=acc,
        autocorr_time=tau,
        param_names=list(problem.names),
        fixed_params=fixed_params,
        converged=converged,
        n_eff=n_eff,
        chi2_best=chi2_best,
        n_data=problem.n_data,
        init_info=init_info,
    )

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. Two 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) Explicit alias for free — use when you want to be unambiguous that the parameter is sampled with hard upper and lower limits.

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', 'a_v').
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 always True (fixed).

Source code in sed_model/params.py 
122
123
124
125
126
def contains(self, v: float) -> bool:
    """True if *v* is within the sampling bounds (free) or always True (fixed)."""
    if self.is_fixed:
        return True
    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 
233
234
235
236
237
238
239
240
241
242
243
244
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 
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
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 = getattr(self, 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 
264
265
266
267
268
269
def in_prior(self, theta: np.ndarray) -> bool:
    """Return True if all free parameters in theta are within their bounds."""
    if len(theta) != self.n_free:
        return False
    unpacked = self.unpack(theta)
    return all(getattr(self, p).contains(unpacked[p]) for p in self.free_names)

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 
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
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 = getattr(self, p)
        mid = 0.5 * (spec.lo + spec.hi)
        c.append(centre[p] if (centre and p in centre) else mid)
        s.append(scatter * (spec.hi - spec.lo))
        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 
138
139
140
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 
143
144
145
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))

bounded 

bounded(lo, hi, name='')

Alias for free — use to be explicit about hard sampling bounds.

Source code in sed_model/params.py 
148
149
150
def bounded(lo: float, hi: float, name: str = "") -> ParamSpec:
    """Alias for ``free`` — use to be explicit about hard sampling bounds."""
    return free(lo, hi, name=name)

fit_params_from_grid 

fit_params_from_grid(grid, *, teff=None, logg=None, meta=None, a_v=None, d_cm=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. All parameters are keyword-only.

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 
374
375
376
377
378
379
380
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
def fit_params_from_grid(
    grid,
    *,
    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,
    a_v:  Union[float, Tuple[float, float], ParamSpec, None] = None,
    d_cm: 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.  All parameters
    are keyword-only.

    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(), converged=None, n_eff=None, chi2_best=None, n_data=None, init_info=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).
converged bool or None

True if the post-burn-in chain is longer than 50 integrated autocorrelation times for every parameter (emcee's rule of thumb). None if the autocorrelation time could not be estimated.
n_eff (ndarray or None, shape(n_free))

Effective number of independent samples per parameter, n_walkers * (n_steps - n_burn) / autocorr_time.
chi2_best float or None

Data chi-square (sum of ((obs-model)/sigma_eff)^2 over detections) of the maximum-posterior sample.
n_data int or None

Number of detected (non-upper-limit) bands used in the fit.
init_info dict

How the walkers were initialised (method, number of modes found, ...) and sampler bookkeeping (number of likelihood evaluations).

summary 

summary()

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

Source code in sed_model/io.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
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

credible_interval 

credible_interval(name, level=0.6827, method='hpd')

Credible interval for one sampled parameter.

method='hpd' (default) returns the shortest interval containing level of the posterior mass. Unlike the equal-tailed interval used by :meth:summary, it still contains the parameter when the posterior piles up against a prior bound (e.g. [M/H] at the grid edge). method='equal' gives the equal-tailed interval.

Source code in sed_model/io.py 
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
def credible_interval(self, name: str, level: float = 0.6827,
                      method: str = "hpd") -> tuple[float, float]:
    """Credible interval for one sampled parameter.

    ``method='hpd'`` (default) returns the shortest interval containing
    ``level`` of the posterior mass.  Unlike the equal-tailed interval
    used by :meth:`summary`, it still contains the parameter when the
    posterior piles up against a prior bound (e.g. [M/H] at the grid
    edge).  ``method='equal'`` gives the equal-tailed interval.
    """
    if name not in self.param_names:
        raise KeyError(f"'{name}' is not a sampled parameter ({self.param_names})")
    if not 0.0 < level < 1.0:
        raise ValueError("level must be in (0, 1)")
    x = np.sort(self.samples[:, self.param_names.index(name)])
    n = len(x)
    if method == "equal":
        a = 0.5 * (1.0 - level)
        return float(np.quantile(x, a)), float(np.quantile(x, 1.0 - a))
    if method != "hpd":
        raise ValueError("method must be 'hpd' or 'equal'")
    k = max(int(np.ceil(level * n)), 1)
    if k >= n:
        return float(x[0]), float(x[-1])
    width = x[k:] - x[:n - k]
    i = int(np.argmin(width))
    return float(x[i]), float(x[i + k])

map_estimate 

map_estimate()

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

Source code in sed_model/io.py 
172
173
174
175
176
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 
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
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]",
        "sigma_sys": "sigma_sys [mag]",
    }

    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}")
    if self.converged is not None:
        neff = ('' if self.n_eff is None else
                '  (n_eff = ' + ' / '.join(f"{x:.0f}" for x in np.ravel(self.n_eff)) + ')')
        print(f"  Converged (N > 50 tau) : {self.converged}{neff}")
    if self.chi2_best is not None and self.n_data is not None:
        print(f"  Best-fit chi2 : {self.chi2_best:.2f} for {self.n_data} bands")

    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 
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
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 = _normalise_npz_path(path)

    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,
        "converged":       self.converged,
        "chi2_best":       self.chi2_best,
        "n_data":          self.n_data,
        "init_info":       self.init_info,
    }

    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
    if self.n_eff is not None:
        arrays["n_eff"] = self.n_eff

    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 
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
@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 = _normalise_npz_path(path)
    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", {}),
        converged=meta.get("converged"),
        n_eff=data["n_eff"] if "n_eff" in data else None,
        chi2_best=meta.get("chi2_best"),
        n_data=meta.get("n_data"),
        init_info=meta.get("init_info", {}),
    )

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 
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
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 
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
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)
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, 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)

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')

Configuration container for extinction settings.

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'

ExtinctionModel 

ExtinctionModel(config=None, **kwargs)

Apply extinction 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::

ext = ExtinctionModel(enabled=True, law='fitzpatrick99', a_v=0.5, r_v=3.1)
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 
675
676
677
678
679
680
681
682
683
684
685
def __init__(self,
             config: Optional[ExtinctionConfig] = None,
             **kwargs) -> None:
    if config is None:
        config = ExtinctionConfig(**kwargs)
    elif kwargs:
        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 
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
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))
    fn = _LAW_DISPATCH[self.config.law]
    return fn(wave, self.config.a_v, self.config.r_v, self.config.gordon23_env)

apply 

apply(wave, flux)

Apply extinction to a flux array.

Parameters:

Name Type Description Default
wave array_like

Wavelengths in Angstroms.

 required
flux array_like

Intrinsic flux in any linear unit.

 required

Returns:

Name Type Description
flux_out ndarray

Extincted flux. If enabled is False this is identical to the input flux (no copy is made).
Source code in sed_model/sed_extinction.py 
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
def apply(self, wave: np.ndarray, flux: np.ndarray) -> np.ndarray:
    """Apply extinction to a flux array.

    Parameters
    ----------
    wave : array_like
        Wavelengths in Angstroms.
    flux : array_like
        Intrinsic flux in any linear unit.

    Returns
    -------
    flux_out : ndarray
        Extincted flux.  If ``enabled`` is False this is identical to
        the input flux (no copy is made).
    """
    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 apply_extinction(a_lam, flux)

remove 

remove(wave, flux)

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

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 
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
def remove(self, wave: np.ndarray, flux: np.ndarray) -> np.ndarray:
    """Remove (de-redden) extinction from an observed flux.

    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)

with_av 

with_av(a_v)

Return a new ExtinctionModel identical to this one but with a different A_V.

Source code in sed_model/sed_extinction.py 
755
756
757
758
def with_av(self, a_v: float) -> "ExtinctionModel":
    """Return a new ExtinctionModel identical to this one but with a different A_V."""
    import dataclasses
    return ExtinctionModel(dataclasses.replace(self.config, a_v=a_v))

disabled classmethod 

disabled()

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

Source code in sed_model/sed_extinction.py 
770
771
772
773
@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 
775
776
777
778
@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 
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
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 
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
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 
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
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)
    M       = _solve_spline_M(_F99_XKNOTS, k_knots)

    opt_mask = x <= _F99_XKNOTS[7]
    uv_mask  = ~opt_mask

    result = np.empty_like(x)

    if opt_mask.any():
        k_opt = _eval_spline(_F99_XKNOTS, k_knots, M, 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 
374
375
376
377
378
379
380
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
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[8]
    uv_mask  = ~opt_mask

    result = np.empty_like(x)

    if opt_mask.any():
        k_opt = _eval_spline(_FM07_XKNOTS, _FM07_KKNOTS, _FM07_SPLINE_M, 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 
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
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

    k = np.zeros_like(w_um)

    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)

    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 
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
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
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)
    k    = np.zeros_like(x)

    ir_mask = x < (1.0 / 3.0)
    if ir_mask.any():
        k[ir_mask] = c6 * x[ir_mask]**alpha - rv

    opt_uv_mask = ~ir_mask
    if opt_uv_mask.any():
        xm  = x[opt_uv_mask]
        xm2 = xm**2
        g2  = gamma**2
        x02 = x0**2
        d   = xm2 / ((xm2 - x02)**2 + xm2*g2)
        k_uv = c1 + c2*xm + c3*d

        fuv = xm > c5
        if fuv.any():
            y = xm[fuv] - c5
            k_uv[fuv] += c4*(0.5392*y**2 + 0.05644*y**3)

        k[opt_uv_mask] = k_uv

    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 
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
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 
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
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, 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
gordon23_env str

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

 'mw'

Returns:

Type Description
ExtinctionModel
Source code in sed_model/sed_extinction.py 
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
def make_extinction_model(
    enabled: bool = False,
    law: str = 'fitzpatrick99',
    r_v: float = 3.1,
    a_v: float = 0.0,
    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.
    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,
        gordon23_env=gordon23_env,
    ))