`virgil.spectra`

Wavelength-dependent fluxes. A component's flux can be a number or a spectrum; inside a System each component is weighted by its spectrum at each sample's wavelength, as in SPARCO.

Wavelength-dependent fluxes for source components.

A component's flux is either a number (achromatic) or a spectrum from this module, which gives its weight at each wavelength. Inside a System the visibility is then V(λ) = Σ f_i(λ) V_i / Σ f_i(λ), as in SPARCO (Kluska et al. 2014):

star = UniformDisk(0.5, flux=PowerLaw(1.0, index=-4.0, wavel0=1.65e-6))
disk = GaussianDisk(5.0, flux=PowerLaw(0.3, index=1.0, wavel0=1.65e-6))
ring = GaussianDisk(5.0, flux=BlackBody(0.3, temperature=1200.0))
wind = GaussianDisk(
    2.0,
    flux=Sum(
        continuum=PowerLaw(0.3, wavel0=2.15e-6),
        brg=GaussianLine(0.5, line_wavel=2.1661e-6, fwhm=1.0e-9),
    ),
)
free = PointSource(flux=Nodes(values, channel_wavelengths))

Spectrum parameters are reached by path like any other, e.g. "disk.flux.ratio", "disk.flux.index", "ring.flux.temperature", "wind.flux.brg.amplitude" or "free.flux.values" (one value per node). Flux ratios are relative: component i's fraction of the total at the reference wavelength is f_i / Σ f, as SPARCO's f_i.

Every spectrum has a reference wavelength wavel0, and its reference flux (spectrum(), used when rendering images) is its value there. Only the evaluated flux must be non-negative: a line or node excess inside a Sum may be negative (absorption) as long as the total is not. is_physical checks the total at each spectrum's characteristic wavelengths (wavel0, nodes and line centres), which is where a sum of these shapes has its minima.

`BlackBody`

Bases: Spectrum

Planck spectrum ratio * B_λ(T, λ) / B_λ(T, wavel0).

The shape of a blackbody at temperature in F_λ, normalized to ratio at wavel0, as SPARCO uses for dust and companions (e.g. Hillen et al. 2016). At long wavelengths (hc/λkT small) it tends to the Rayleigh-Jeans PowerLaw with index -4.

Parameters:

Name	Type	Description	Default
`ratio`	`float or array - like`	Flux at the reference wavelength, relative to the other components.	required
`temperature`	`float or array - like`	Temperature in kelvin.	required
`wavel0`	`float or array - like`	Reference wavelength in metres (default 1.65e-6, H band).	`1.65e-06`

Examples:

>>> round(float(BlackBody(0.2, 1500.0)(1.65e-6)), 6)
0.2

`ratio = np.asarray(ratio, dtype=float)` `instance-attribute`

`temperature = np.asarray(temperature, dtype=float)` `instance-attribute`

`wavel0 = np.asarray(wavel0, dtype=float)` `instance-attribute`

`init(ratio, temperature, wavel0=1.65e-06)`

`call(wavel=None)`

`is_physical()`

`__check_init__()`

`GaussianLine`

Bases: _Line

A Gaussian emission or absorption line.

amplitude * exp(-4 ln 2 ((λ - line_wavel) / fwhm)²): the flux at the line centre, with the full width at half maximum fwhm. Its integral over wavelength is amplitude * fwhm * sqrt(π / (4 ln 2)).

On its own a line is the whole flux of a component, so amplitude must be non-negative; inside a Sum with a continuum a negative amplitude is an absorption line.

Parameters:

Name	Type	Description	Default
`amplitude`	`float`	Flux at the line centre, relative to the other components.	required
`line_wavel`	`float`	Line centre in metres (shifted by any velocity: λ₀(1 + v/c)).	required
`fwhm`	`float`	Full width at half maximum in metres.	required
`wavel0`	`float`	Reference wavelength in metres; the line centre by default.	`None`

Examples:

>>> line = GaussianLine(0.4, line_wavel=2.1661e-6, fwhm=1.0e-9)
>>> round(float(line(2.1661e-6 + 0.5e-9)), 4)
0.2

`LorentzianLine`

Bases: _Line

A Lorentzian emission or absorption line.

amplitude / (1 + 4 ((λ - line_wavel) / fwhm)²): the flux at the line centre, with the full width at half maximum fwhm. Its integral over wavelength is amplitude * π fwhm / 2. Its wings fall off slowly, so over a wide band it adds a near-constant pedestal.

Parameters are as for GaussianLine.

Examples:

>>> line = LorentzianLine(0.4, line_wavel=2.1661e-6, fwhm=1.0e-9)
>>> round(float(line(2.1661e-6 + 0.5e-9)), 4)
0.2

`Nodes`

Bases: Spectrum

A spectrum through free values at fixed wavelengths.

Linear or natural-cubic interpolation between the nodes. Beyond the end nodes the flux is either held at the end values (outside="constant") or a fixed number (e.g. outside=0.0 for an excess that lives only in a line window, on top of a continuum in a Sum). With outside=0.0 the continuum is fixed by the channels outside the window, so the two are identifiable; give the end nodes the value 0 (or fix them there) to keep the spectrum continuous.

Fit the values with a prior of their shape; a smooth spectrum is a prior on them (e.g. a Gaussian process over wavelength), not a feature of this class.

Parameters:

Name	Type	Description	Default
`values`	`(array - like, shape(n))`	Flux at each node, relative to the other components. They may be negative (an absorption excess) as long as a component's total flux is not.	required
`wavel`	`(array - like, shape(n))`	Node wavelengths in metres: finite, positive and strictly increasing. Fixed (not fitted) in practice.	required
`kind`	`(linear, cubic)`	Interpolation; `"cubic"` is a natural cubic spline and needs at least three nodes.	`"linear"`
`outside`	`constant or float`	The flux beyond the end nodes (default `"constant"`).	`'constant'`
`wavel0`	`float`	Reference wavelength in metres; the first node by default.	`None`

Examples:

>>> spectrum = Nodes([0.2, 0.4], [2.0e-6, 2.2e-6])
>>> round(float(spectrum(2.1e-6)), 6)
0.3
>>> excess = Nodes([0.0, 0.5, 0.0], [2.16e-6, 2.166e-6, 2.172e-6], outside=0.0)
>>> float(excess(2.0e-6))
0.0

`values = np.asarray(values, dtype=float)` `instance-attribute`

`wavel = np.asarray(wavel, dtype=float)` `instance-attribute`

`kind = kind` `class-attribute` `instance-attribute`

`outside = outside if outside == 'constant' else float(outside)` `class-attribute` `instance-attribute`

`wavel0 = self.wavel[..., 0] if wavel0 is None and self.wavel.ndim == 1 and self.wavel.size else np.asarray(wavel0, dtype=float)` `instance-attribute`

`init(values, wavel, kind='linear', outside='constant', wavel0=None)`

`__check_init__()`

`PowerLaw`

Bases: Spectrum

Power-law spectrum ratio * (λ / wavel0) ** index.

Parameters:

Name	Type	Description	Default
`ratio`	`float or array - like`	Flux at the reference wavelength, relative to the other components.	required
`index`	`float or array - like`	Spectral index (default 0, i.e. achromatic). A star in the Rayleigh-Jeans regime of F_λ has index -4.	`0.0`
`wavel0`	`float or array - like`	Reference wavelength in metres (default 1.65e-6, H band).	`1.65e-06`

Examples:

>>> round(float(PowerLaw(0.2, index=-4.0, wavel0=1.6e-6)(3.2e-6)), 6)
0.0125

`ratio = np.asarray(ratio, dtype=float)` `instance-attribute`

`index = np.asarray(index, dtype=float)` `instance-attribute`

`wavel0 = np.asarray(wavel0, dtype=float)` `instance-attribute`

`init(ratio, index=0.0, wavel0=1.65e-06)`

`call(wavel=None)`

`is_physical()`

`__check_init__()`

`Spectrum`

Bases: Base

Base class for component spectra.

Subclasses implement _at(wavel), the flux at wavel (metres, any shape), and have a reference wavelength wavel0; calling the spectrum with no wavelength gives its reference flux, the value at wavel0 (as used when rendering images). _params_valid checks the parameters' domains, and _check_wavel lists the wavelengths where the flux must be non-negative.

`call(wavel=None)`

`is_physical()`

Whether the spectrum is valid, as a (traceable) boolean.

Its parameters are in their domains and its flux is non-negative at its characteristic wavelengths (wavel0, nodes, line centres).

`Sum`

Bases: Spectrum

The sum of named spectra, e.g. a continuum plus lines.

Parts are reached by name, like the components of a System: "star.flux.continuum.index", "star.flux.brg.amplitude". Only the total must be non-negative, so a part may be an absorption line or a negative node excess.

Parameters:

Name	Type	Description	Default
`parts`	`dict`	`{name: Spectrum}`, positionally; or give them by keyword.	`None`
`wavel0`	`float`	Reference wavelength in metres, where the reference flux is the sum of the parts; the first part's `wavel0` by default.	`None`

Examples:

>>> flux = Sum(
...     continuum=PowerLaw(1.0, wavel0=2.2e-6),
...     brg=GaussianLine(-0.3, line_wavel=2.1661e-6, fwhm=1.0e-9),
... )
>>> round(float(flux(2.1661e-6)), 4)  # continuum minus the absorption
0.7
>>> round(float(flux()), 4)  # the reference flux, at 2.2 µm
1.0

`names = tuple(parts)` `class-attribute` `instance-attribute`

`parts = tuple(parts.values())` `instance-attribute`

`wavel0 = np.asarray(wavel0, dtype=float)` `instance-attribute`

`components` `property`

The parts as a {name: spectrum} dictionary, in order.

`init(parts=None, /, *, wavel0=None, **named)`

`getattr(name)`

`Tabulated`

Bases: Spectrum

A free flux in every spectral channel, interpolated linearly between.

Deprecated: use Nodes, which also has cubic interpolation, a fixed value outside the nodes, and a reference flux at wavel0 like every other spectrum. Tabulated(ratio, wavel) is Nodes(ratio, wavel) except for its reference flux (the mean over the nodes) and its parameter name (ratio). It is kept so that existing scripts run unchanged, and will be removed in a later release.

For fitting a spectrum channel by channel, e.g. a companion's flux ratio across emission lines: give wavel the data's channel wavelengths and fit ratio (one value per channel) with a prior of that shape.

Parameters:

Name	Type	Description	Default
`ratio`	`(array - like, shape(n))`	Flux at each node, relative to the other components: finite and non-negative, with n >= 1.	required
`wavel`	`(array - like, shape(n))`	Node wavelengths in metres: finite, positive and strictly increasing. Beyond the end nodes the flux is constant.	required

Notes

The reference flux (wavel=None, used when rendering) is the mean over the nodes.

Examples:

>>> spectrum = Tabulated([0.2, 0.4], [2.0e-6, 2.2e-6])
>>> round(float(spectrum(2.1e-6)), 6)
0.3

`ratio = np.asarray(ratio, dtype=float)` `instance-attribute`

`wavel = np.asarray(wavel, dtype=float)` `instance-attribute`

`init(ratio, wavel)`

`call(wavel=None)`

`is_physical()`

`__check_init__()`

`flux_at(flux, wavel=None)`

Evaluate a number or a spectrum at wavel (None = reference).

`reference_flux(flux)`

The reference flux of a number or a spectrum.

For a spectrum this is flux(None), its value at wavel0 (the ratio for PowerLaw and BlackBody); only the deprecated Tabulated uses the mean over its nodes instead.

Tabulated is deprecated: it still works unchanged (and emits a DeprecationWarning), but new code should use Nodes, which adds cubic interpolation, a fixed value outside the nodes and a reference flux at wavel0. It is not exported from the top-level virgil namespace (import it from virgil.spectra) and will be removed in a later release.

virgil.spectra

BlackBody

ratio = np.asarray(ratio, dtype=float) instance-attribute

temperature = np.asarray(temperature, dtype=float) instance-attribute

wavel0 = np.asarray(wavel0, dtype=float) instance-attribute

__init__(ratio, temperature, wavel0=1.65e-06)

__call__(wavel=None)

is_physical()

__check_init__()

GaussianLine

LorentzianLine

Nodes

values = np.asarray(values, dtype=float) instance-attribute

wavel = np.asarray(wavel, dtype=float) instance-attribute

kind = kind class-attribute instance-attribute

outside = outside if outside == 'constant' else float(outside) class-attribute instance-attribute

wavel0 = self.wavel[..., 0] if wavel0 is None and self.wavel.ndim == 1 and self.wavel.size else np.asarray(wavel0, dtype=float) instance-attribute

__init__(values, wavel, kind='linear', outside='constant', wavel0=None)

__check_init__()

PowerLaw

ratio = np.asarray(ratio, dtype=float) instance-attribute

index = np.asarray(index, dtype=float) instance-attribute

wavel0 = np.asarray(wavel0, dtype=float) instance-attribute

__init__(ratio, index=0.0, wavel0=1.65e-06)

__call__(wavel=None)

is_physical()

__check_init__()

Spectrum

__call__(wavel=None)

is_physical()

Sum

names = tuple(parts) class-attribute instance-attribute

parts = tuple(parts.values()) instance-attribute

wavel0 = np.asarray(wavel0, dtype=float) instance-attribute

components property

__init__(parts=None, /, *, wavel0=None, **named)

__getattr__(name)

Tabulated

ratio = np.asarray(ratio, dtype=float) instance-attribute

wavel = np.asarray(wavel, dtype=float) instance-attribute

__init__(ratio, wavel)

__call__(wavel=None)

is_physical()

__check_init__()

flux_at(flux, wavel=None)

reference_flux(flux)

`virgil.spectra`

`BlackBody`

`ratio = np.asarray(ratio, dtype=float)` `instance-attribute`

`temperature = np.asarray(temperature, dtype=float)` `instance-attribute`

`wavel0 = np.asarray(wavel0, dtype=float)` `instance-attribute`

`init(ratio, temperature, wavel0=1.65e-06)`

`call(wavel=None)`

`is_physical()`

`__check_init__()`

`GaussianLine`

`LorentzianLine`

`Nodes`

`values = np.asarray(values, dtype=float)` `instance-attribute`

`wavel = np.asarray(wavel, dtype=float)` `instance-attribute`

`kind = kind` `class-attribute` `instance-attribute`

`outside = outside if outside == 'constant' else float(outside)` `class-attribute` `instance-attribute`

`wavel0 = self.wavel[..., 0] if wavel0 is None and self.wavel.ndim == 1 and self.wavel.size else np.asarray(wavel0, dtype=float)` `instance-attribute`

`init(values, wavel, kind='linear', outside='constant', wavel0=None)`

`__check_init__()`

`PowerLaw`

`ratio = np.asarray(ratio, dtype=float)` `instance-attribute`

`index = np.asarray(index, dtype=float)` `instance-attribute`

`wavel0 = np.asarray(wavel0, dtype=float)` `instance-attribute`

`init(ratio, index=0.0, wavel0=1.65e-06)`

`call(wavel=None)`

`is_physical()`

`__check_init__()`

`Spectrum`

`call(wavel=None)`

`is_physical()`

`Sum`

`names = tuple(parts)` `class-attribute` `instance-attribute`

`parts = tuple(parts.values())` `instance-attribute`

`wavel0 = np.asarray(wavel0, dtype=float)` `instance-attribute`

`components` `property`

`init(parts=None, /, *, wavel0=None, **named)`

`getattr(name)`

`Tabulated`

`ratio = np.asarray(ratio, dtype=float)` `instance-attribute`

`wavel = np.asarray(wavel, dtype=float)` `instance-attribute`

`init(ratio, wavel)`

`call(wavel=None)`

`is_physical()`

`__check_init__()`

`flux_at(flux, wavel=None)`

`reference_flux(flux)`