Composition: System and building blocks

Bases: SourceModel

Flux-weighted mixture of named source models.

A System is how you describe a scene with more than one part: a star with a disk, a binary inside a ring, a companion with its own circumstellar material. Each component is given a name, and the visibility is the flux-weighted mean \(V = \sum_i f_i V_i \,/\, \sum_i f_i\), where f_i is each component's flux. Because interferometric visibilities are normalized to 1 at zero baseline, only the ratios of the fluxes can be measured: keep one reference component (usually the star) at flux=1 and fit the others relative to it.

Components are reached by name, both as attributes (system.comp.flux) and as zodiax paths (system.get("comp.flux"), system.set("comp.flux", 0.02)). These paths are how the fitting tools in virgil.grid_fit and numpyro_model address parameters. Components keep the order in which they were given.

A System can itself be a component. Its own flux is then the total flux of the group relative to its siblings, and dra/ddec move the whole group together.

Parameters:

Examples:

A star with a faint companion:

>>> binary = System(
...     star=PointSource(),
...     comp=PointSource(dra=45.0, ddec=30.0, flux=0.01),
... )
>>> binary.comp
PointSource(flux=0.01, dra=45, ddec=30)

A star with a rim, and a companion that has its own disk:

>>> scene = System(
...     star=PointSource(),
...     rim=ModulatedGaussianRim(diam=40.0, fwhm=4.0, inc=50.0, pa=30.0, flux=0.5),
...     comp=System(
...         core=PointSource(),
...         disk=GaussianDisk(sigma=4.0, flux=0.5),
...         dra=-30.0,
...         ddec=25.0,
...         flux=0.3,
...     ),
... )
>>> list(scene.components)
['star', 'rim', 'comp']
>>> moved = scene.set(["comp.dra", "comp.ddec"], [30.0, -25.0])

Bases: SourceModel

Base class for single shapes with flux, dra and ddec.

A component on its own is normalized to unit flux. Inside a System, flux is its weight relative to the other components, and dra/ddec place its centre (milliarcseconds, positive dra to the East, positive ddec to the North).

New shapes subclass Component and implement _centred_cvis (the unit-flux visibility of the shape at the origin) and _centred_image (an un-normalized image of the shape at the origin); offsets and mixing are handled here.

Bases: Component

Unresolved point source.

Parameters:

Examples:

>>> star = PointSource()
>>> companion = PointSource(flux=0.01, dra=45.0, ddec=30.0)

Bases: Component

Circular Gaussian brightness distribution.

Parameters:

Examples:

>>> halo = GaussianDisk(sigma=8.0, flux=0.2)

Bases: Component

Elliptical Gaussian brightness distribution.

Parameters:

Examples:

>>> plume = EllipticalGaussian(fwhm=30.0, ratio=0.5, pa=20.0, flux=2.0)

Bases: Component

A Gaussian ridge bent along a circular arc, e.g. a curved shock front.

The brightness is a curve convolved with a circular Gaussian of FWHM width. The curve is a circle of radius about the centre of curvature (dra, ddec), weighted along its length by a Gaussian of FWHM length (an arc length, in milliarcseconds) peaking at position angle pa as seen from the centre. For a large radius it tends to an elliptical Gaussian with FWHMs hypot(length, width) along the arc and width across it.

The weight is a Gaussian in the arc length s from pa, taken once round the circle, -π radius <= s <= π radius: it is cut off at the antipode of pa, where the two tails meet, and never wraps round the circle more than once. The quadrature covers ±6σ (the part beyond, a fraction 2e-9 of the flux, is dropped), or the whole circle once 6σ >= π radius, i.e. length >= 0.39 π radius. The cut at the antipode matters only for longer arcs: at length = π radius it removes 2 % of the Gaussian, and as length grows the brightness tends to a uniform ring.

The visibility is the Fourier transform of the curve, a line integral evaluated by the trapezoidal rule on nodes equally spaced points spanning ±min(6σ, π radius), times the Gaussian's. It is accurate while the spacing of the points, 2 min(6σ, π radius) / (nodes - 1), is below half the shortest fringe spacing, λ / 2 B_max.

Parameters:

Examples:

>>> shock = GaussianArc(radius=30.0, width=2.0, length=40.0, pa=270.0)

Bases: Component

A thin conical shell, truncated near its apex, e.g. a dust cone.

The cone's axis points towards position angle pa and is tilted tilt degrees out of the sky plane; its half-opening angle is alpha. Its apex lies tip milliarcseconds (along the axis, in 3-D) behind (dra, ddec), opposite to pa. The emission starts a slant distance s0 from the apex along the walls and falls off as exp(-(s - s0) / length), and the shell has a Gaussian thickness of FWHM width. It is optically thin, so the sign of the tilt does not change the image of an unmodulated cone.

The cone is a stack of rings about its axis. A ring of radius ρ in the plane perpendicular to an axis tilted β out of the sky projects to an ellipse, whose visibility is J0(2π ρ q) with q² = q_perp² + (ratio · q_par sin β)², where q_par and q_perp are the spatial frequencies along and across the projected axis, times the phase of the ring's centre, which lies (s cos α - tip) cos β along the projected axis. The rings are weighted by the area element (∝ ρ) and the emissivity, and integrated over s from s0 to s0 + 5 length (the last 0.7 % of the flux is dropped) by the midpoint rule on n_rings rings.

Azimuthal modulation. az_amps and az_pas brighten one side of the walls, e.g. the leading edge of a colliding-wind shock that the orbit sweeps round. Every ring's brightness is multiplied by 1 + Σ_m A_m cos(m (φ - φ_m)) in its own azimuth φ, with A_m in az_amps and φ_m in az_pas, the same modulation along the whole length of the cone. Like those of ModulatedGaussianRim, the azimuths are measured in the plane of the ring, in the same sense as position angle. φ = pa + 90 and pa - 90 are the two walls seen across the projected axis, at those position angles on the sky. φ = pa is the side of each ring that projects towards pa on the sky when tilt > 0 (towards pa + 180 when tilt < 0). The modulation stays analytic: by the Jacobi-Anger expansion a ring adds A_m (-i)^m J_m(2π ρ q) cos(m (θ - φ_m)) to J0, with θ the direction of the spatial frequency in the ring's plane. A modulated cone no longer looks the same at tilt and -tilt: the image at -tilt is the one at tilt with az_pas mirrored to 2 pa + 180 - az_pas. That leaves the walls (pa ± 90) unchanged. Bound to an orbit's frame by Attached (bind={"az_pas": ...}), a sky position angle is deprojected into the rings' azimuth, so the bright side points at it on the sky.

Choosing n_rings. The quadrature is second order: once the rings are fine enough to resolve the fringes, the error in the visibility falls as 1 / n_rings**2, so each doubling of n_rings cuts it by about 4. It is fine enough when both the spacing of the rings' centres on the sky, 5 length cos α cos β / n_rings, and the step between their radii, 5 length sin α max(1, ratio) / n_rings, are below half the shortest fringe spacing (a cone seen down its axis, tilt 90°, has all its centres together, and only the radius step matters). That criterion only says the error is small, not that it is below your noise. For a cone with length 13.8 mas and alpha 62.5°, over baselines out to 0.3 cycles/mas, max |V(n) - V(2n)| is about 6e-4 for n = 32, 1.5e-4 for 64, 4e-5 for 128 and 1e-5 for 256, so the error of n_rings = 32 itself is about 8e-4 in visibility amplitude. That is negligible for noisy data but not for well-measured data: a high-S/N GRAVITY dataset gained about 4 in log-likelihood per epoch going from 32 to 64 rings at fixed parameters, and nearly 29 over three epochs from 24 to 64.

To check, refit or evaluate at the best fit with n_rings doubled and compare χ² (or the log-likelihood): if |Δχ²| ≳ 1 per dataset (equivalently |Δ log L| ≳ 0.5), use more rings, and double again until it is below that. Well-measured data (e.g. GRAVITY) may need 64 or more. The cost is linear in n_rings.

Parameters:

Examples:

>>> cone = TruncatedCone(tip=5.0, alpha=30.0, s0=4.0, length=10.0,
...                      width=1.0, tilt=20.0, pa=90.0)

A cone brighter on its southern wall (PA 180), half as bright on the northern one:

>>> lopsided = TruncatedCone(tip=5.0, alpha=30.0, s0=4.0, length=10.0,
...                          width=1.0, tilt=20.0, pa=90.0,
...                          az_amps=1 / 3, az_pas=180.0)

Bases: Component

Uniformly bright (tophat) circular disk, e.g. a resolved stellar photosphere.

Parameters:

Examples:

>>> photosphere = UniformDisk(diam=3.0)

Bases: _LimbDarkenedDisk

Circular disk with polynomial limb darkening of order up to 22.

The brightness is \(I(\mu) / I(1) = 1 - \sum_{n=1}^{N} u_n (1 - \mu)^n\), with \(\mu = \sqrt{1 - (r / R)^2}\) the cosine of the angle between the line of sight and the surface normal. This is the convention of jaxoplanet, starry and harmonix, so u is the same as a jaxoplanet Surface's u: u=(u1,) is the linear law, u=(u1, u2) the quadratic law, and the default u=() a uniform disk. The order is at most 22, as the Bessel functions needed are of order at most 12. To fit the quadratic law with priors that cover exactly the physical profiles, use QuadraticLimbDarkenedDisk.

The visibility is analytic (Quirrenbach et al. 1996, eq. 3; see cvis_limb_darkened_disk). harmonix (Dholakia & Pope 2025) generalizes the same result to limb-darkened spherical-harmonic maps.

Parameters:

Examples:

>>> linear = LimbDarkenedDisk(3.0, u=[0.6])
>>> quadratic = LimbDarkenedDisk(3.0, u=[0.4, 0.25])

Bases: _LimbDarkenedDisk

Circular disk with quadratic limb darkening in Kipping's (2013) \(q_1, q_2\).

The brightness is \(I(\mu) / I(1) = 1 - u_1 (1 - \mu) - u_2 (1 - \mu)^2\), with \(u_1 = 2 \sqrt{q_1}\, q_2\) and \(u_2 = \sqrt{q_1}\,(1 - 2 q_2)\) (Kipping 2013, eqs. 15-16). Every \((q_1, q_2)\) in the unit square gives a profile that is positive and decreases from the centre to the limb, and every such profile has one, so uniform priors on \([0, 1]\) for both are uninformative over exactly the physical laws. The visibility is analytic (Quirrenbach et al. 1996, eq. 3; see cvis_limb_darkened_disk).

Parameters:

Examples:

>>> import numpyro.distributions as dist
>>> star = QuadraticLimbDarkenedDisk(3.0, q1=0.4, q2=0.3)
>>> priors = {
...     "diam": dist.LogUniform(2.0, 4.0),  # a scale: log-uniform
...     "q1": dist.Uniform(0.0, 1.0),  # Kipping's uniform prior
...     "q2": dist.Uniform(0.0, 1.0),
... }

From tabulated \(u_1, u_2\) (e.g. Claret's tables):

>>> star = QuadraticLimbDarkenedDisk.from_u(3.0, u1=0.4, u2=0.25)

Bases: _LimbDarkenedDisk

Circular disk with square-root limb darkening in Kipping's (2013) \(q_1, q_2\).

The brightness is \(I(\mu) / I(1) = 1 - c (1 - \mu) - d (1 - \sqrt{\mu})\) (Díaz-Cordovés & Giménez 1992), which suits late-type stars in the near-infrared better than the quadratic law (van Hamme 1993), with \(c = \sqrt{q_1}\,(1 - 2 q_2)\) and \(d = 2 \sqrt{q_1}\, q_2\), inverting Kipping (2013), eqs. 23-24. As for QuadraticLimbDarkenedDisk, the unit square in \((q_1, q_2)\) is exactly the set of positive profiles that decrease towards the limb, so uniform priors on \([0, 1]\) are uninformative over the physical laws. The \(\sqrt{\mu}\) term needs a Bessel function of order \(5/4\) in the analytic visibility (Quirrenbach et al. 1996, eq. 3; see cvis_limb_darkened_disk).

Parameters:

Examples:

>>> star = SquareRootLimbDarkenedDisk(3.0, q1=0.5, q2=0.4)
>>> same = SquareRootLimbDarkenedDisk.from_cd(3.0, c=star.c, d=star.d)

Bases: _LimbDarkenedDisk

Limb-darkened disk with an elliptical outline, e.g. an oblate star.

A LimbDarkenedDisk compressed along its minor axis: the brightness is \(I(\mu) / I(1) = 1 - \sum_{n=1}^{N} u_n (1 - \mu)^n\), with \(\mu = \sqrt{1 - \rho^2}\) and \(\rho\) the elliptical radius (1 on the limb), so the profile follows the outline. The major axis is oriented exactly as an EllipticalGaussian's. Like every component it has unit flux on its own: compressing the disk does not change its total flux, so the visibility at zero baseline is 1 and flux is its weight in a System.

The visibility is that of the circular disk at the spatial frequencies stretched by ratio along the minor axis (the Fourier similarity theorem), so it is analytic like cvis_limb_darkened_disk. This is a geometric ellipse with a radial profile, not a model of a rapid rotator's gravity darkening; for that use GravityDarkenedStar.

Parameters:

Examples:

>>> star = EllipticalLimbDarkenedDisk(5.5, ratio=0.8, pa=30.0, u=[0.5])

Bases: Component

Rapidly rotating star: Roche shape and gravity darkening (ELR11).

The model of Espinosa Lara & Rieutord (2011, A&A 533, A43): a rigidly rotating star has the oblate Roche shape, and its local bolometric flux follows the effective gravity without a free gravity-darkening exponent \(\beta\). The brightness is proportional to that local flux (no limb darkening yet), summed over the visible triangles of a surface mesh. Ported from Shashank Dholakia's jax-interferometry (ELR_Model, commit 70689ed); the physics lives in virgil._elr.

Parameters:

star = GravityDarkenedStar(
    1.0, omega=0.9, inc=45.0, t_pole=9000.0
)
companion = PointSource(flux=BlackBody(0.01, 3000.0))
system = System(star=star, companion=companion)

Examples:

>>> star = GravityDarkenedStar(2.0, omega=0.8, inc=60.0, pa=30.0)
>>> v0 = star.model(np.zeros(1), np.zeros(1), 1.65e-6)
>>> round(float(np.abs(v0)[0]), 3)
1.0

The chromatic model's pole-to-equator contrast depends on wavelength:

>>> hot = GravityDarkenedStar(2.0, omega=0.9, inc=45.0, t_pole=9000.0)
>>> w = hot._weight(np.array([1.0e-6, 2.2e-6]))
>>> bool(w[0] > 1.0 > w[1])
True

Bases: Component

Azimuthally modulated, infinitely thin rim convolved with a 2D Gaussian that is isotropic in the plane of the rim, optionally inclined and rotated.

Parameters:

Examples:

>>> rim = ModulatedGaussianRim(
...     diam=40.0, fwhm=4.0, inc=50.0, pa=30.0, az_amps=0.7, az_pas=120.0
... )

Bases: Component

Flared, inclined scattered-light disk (Blakely et al. 2024, §III).

The geometrical disk model that Blakely et al. (2024, arXiv:2404.13032) fitted to JWST AMI data of PDS 70: a skewed Gaussian ring on a flared surface, with a forward-scattering peak on its near side. It is an abstract base: use FlaredDiskHG, FlaredDiskGaussian or FlaredDiskPowerLaw, which differ only in the azimuthal phase function. The disk contains no star or planets; compose them in a System, and any over-resolved flux with a Resolved component (the paper's \(I_o\)).

The brightness has no analytic Fourier transform, so the visibilities are the exact Fourier transform of the brightness sampled on a fixed, centred grid of npix x npix pixels of pixel_scale_mas. The grid must cover the whole disk, and its pixels must be small enough to resolve the ring and its sharp inner edge.

Parameters:

Bases: FlaredDisk

FlaredDisk with a Henyey–Greenstein phase function.

\(f(\theta) = \dfrac{1 - g^2}{4\pi\,(1 + g^2 - 2g\cos\theta)^{3/2}}\) (Blakely et al. 2024, Eq. 6).

Parameters:

Examples:

>>> disk = FlaredDiskHG(
...     g=0.3, radius=440.0, fwhm=290.0, inc=52.0, pa=160.0,
...     npix=96, pixel_scale_mas=20.0, flux=0.05,
... )

Bases: FlaredDisk

FlaredDisk with a Gaussian phase function.

\(f(\theta) = \exp\left(-\theta^2 / 2\sigma_\theta^2\right)\), with \(\theta \in (-180°, 180°]\) (Blakely et al. 2024, Eq. 7).

Parameters:

Bases: FlaredDisk

FlaredDisk with a power-law phase function.

\(f(\theta) = \cos^N(\theta / 2)\), with \(\theta \in (-180°, 180°]\) (Blakely et al. 2024, Eq. 8), the best-fitting form for PDS 70.

Parameters:

Bases: SourceModel

Fully resolved (over-resolved) flux, e.g. a large, diffuse envelope.

Its visibility is 0 on every non-zero baseline, so inside a System it only adds to the normalization, lowering every other component's visibility by the same factor.

Parameters:

Examples:

>>> scene = System(star=PointSource(), background=Resolved(flux=0.1))

Bases: SourceModel

A model rotated on the sky about the phase centre.

The rotation is by rotation_deg from North towards East, and the angle is an ordinary (fittable) parameter, unlike Image's rotation_deg, which only orients its pixel grid. Use it for a scene seen at several epochs, e.g. a spiral rotating between them (see fit with a model per dataset).

Parameters:

Examples:

A companion to the North, rotated by 90°, lands to the East:

>>> import jax.numpy as np
>>> from virgil.models import GaussianDisk, Rotated, System
>>> north = System(a=GaussianDisk(1.0), b=GaussianDisk(1.0, ddec=10.0))
>>> east = System(a=GaussianDisk(1.0), b=GaussianDisk(1.0, dra=10.0))
>>> u, v = np.array([3.0, 5.0]), np.array([1.0, -2.0])
>>> rotated = Rotated(north, 90.0).model(u, v, 1e-6)
>>> bool(np.allclose(rotated, east.model(u, v, 1e-6), atol=1e-6))
True

Bases: SourceModel

A component placed and oriented in the frame of a binary's orbit.

At each time the component is moved to its anchor on the orbit, and its angle attributes are set from angles of the binary frame (see KeplerOrbit.frame). A companion on its orbit is Attached(PointSource(flux), orbit); a disc around it in the orbital plane, brighter on the side facing the primary, is

Attached(ModulatedGaussianRim(...), orbit, bind={"pa": "node_pa", "inc": "apparent_inc", "az_pas": "towards_primary"}).

The model changes with time, so it is evaluated through at; inside a System, OIData.model evaluates every sample at its own time. Components that should share one orbit (a companion and its disc) are best built from shared parameters in a model function (see fit).

Parameters:

Bases: SourceModel

A primary and a point-source companion moving on a Keplerian orbit.

The time-dependent form of BinaryModelCartesian: at each time, at returns the BinaryModelCartesian with the companion at the orbit's sky position, so each snapshot keeps the fast binary path. It is the same scene as System(primary=PointSource(), companion=Attached( PointSource(flux), orbit)), for the common case of two unresolved stars. The primary is the scene's origin and reference (flux is the companion's flux relative to it).

Evaluate it with one snapshot per dataset or epoch through Epochs, or per sample by OIData.model.

Parameters:

Examples:

>>> from virgil.models import OrbitalBinary
>>> from virgil.orbits import KeplerOrbit
>>> orbit = KeplerOrbit(700.0, 0.0, 0.3, 50.0, 60.0, 120.0, 20.0,
...                     t_ref=60500.0)
>>> snap = OrbitalBinary(orbit, 0.1).at(60500.0)
>>> type(snap).__name__, round(float(snap.flux), 2)
('BinaryModelCartesian', 0.1)