Quickstart Guide

This guide introduces the basic workflow for using nyx to simulate the night sky brightness, then shows how to fit the model to data and save the result.

Basic Concepts

nyx models astronomical observations through a small set of composable components, all sharing a common Geometry:

  1. Instrument — telescope characteristics (effective aperture, bandpass, pixel response).

  2. Atmosphere — atmospheric scattering and absorption.

  3. Emitters — sky brightness sources (zodiacal light, airglow, moon, stars).

  4. Observation — observer location, time(s), and pointing.

These are combined into a Scene, which renders to per-pixel photon detection rates.

Working with Units

nyx is Astropy-aware: physical inputs can be supplied as Quantity objects and are converted internally.

import astropy.units as u

fov = 3.5 * u.deg          # field of view half-angle
wvls = [300, 400, 500] * u.nm

Configuring JAX

nyx does not change JAX’s global configuration on import. The render integrates over wavelength and sky position, so on GPU and TPU the reduced-precision matmuls JAX allows by default cost accuracy.

import nyx

nyx.configure()   # sets jax_default_matmul_precision to "highest"

The Geometry

The Geometry defines the resolution and grid shared by every component:

import jax.numpy as jnp
import astropy.units as u

from nyx.core.geometry import Geometry

geo = Geometry(
    wvls=jnp.linspace(300, 700, 50) * u.nm,  # wavelength grid (nm)
    nside=16,                                # HEALPix resolution of the sky hemisphere
    ngrid=2,                                 # in-scattering grid points per FOV axis
    fov=3.5 * u.deg,                         # field of view half-angle
)

Higher nside and ngrid increase accuracy at the cost of compute and memory; start small and increase as needed.

A Complete Example

The following builds a full scene and renders it to pixel rates. The emitters use convenience constructors that load their backing datasets (Leinert 1998, ESO SkyCalc, Gaia DR3, Jones 2013):

import jax.numpy as jnp
import astropy.units as u
from astropy.coordinates import EarthLocation, SkyCoord
from astropy.time import Time

import nyx
from nyx import (
    Airglow, EffectiveApertureInstrument, Geometry, Moon, Observation,
    Scene, SingleScattering, Stars, ZodiacalLight,
)

nyx.configure()

# 0. Shared render geometry
geo = Geometry(wvls=jnp.linspace(300, 700, 50) * u.nm,
               nside=16, ngrid=2, fov=3.5 * u.deg)

# 1. Instrument, from an iactrace effective-aperture scan
instrument = EffectiveApertureInstrument.from_iactrace_table(
    geo, "CT1_aperture.npz")

# 2. Atmosphere
atmosphere = SingleScattering.from_hg(geo)

# 3. Emitters, keyed by name
sources = {
    "airglow": Airglow.from_eso_skycalc(geo, sfu=100.0),
    "zodiacal": ZodiacalLight.from_leinert1998(geo),
    "GaiaDR3": Stars.from_gaia_dr3(geo, lim_mag=12),
    "moon": Moon.from_jones2013(geo),
}

# 4. Observation: location, time(s), target, geometry
location = EarthLocation.of_site("Roque de los Muchachos")
times = Time(["2026-06-05T23:00:00"])
target = SkyCoord(ra=83.6 * u.deg, dec=22.0 * u.deg)
obs = Observation(location, times, target, geo)

# 5. Build the scene
scene = Scene.build(instrument, atmosphere, sources, obs)

# 6. Render to per-pixel photon rates
rates = scene.render()

Scene.render() returns a dict keyed by instrument name; each value is a JAX array of photon detection rates per pixel [photon/s]. A single instrument is automatically named "instrument":

pixel_rates = rates["instrument"]  # shape (n_obs, n_pixels)

Fitting the Model to Data

The whole pipeline is differentiable, so you can fit model parameters to measured pixel rates with gradient-based optimisation. Parameters are Parameter objects; use freeze() / unfreeze() (with selector lambdas that address the parameter in the scene pytree) to choose what is trained.

import optimistix as optx
from nyx.core import freeze_all, unfreeze
from nyx.infer import Optimizer

# Start from everything frozen, then unfreeze what you want to fit.
scene = freeze_all(scene)
scene = unfreeze(scene, lambda s: s.instruments["instrument"].pixel_efficiency)

# Residuals against measured rates (least-squares).
def residuals(scene):
    return scene.render()["instrument"] - measured_rates

opt = Optimizer(residuals, optx.LevenbergMarquardt(rtol=1e-6, atol=1e-6))
fitted, sol = opt.run(scene)

# 1-sigma parameter uncertainties at the solution
errors = opt.errors(fitted)

For jointly fitting several targets/observations, see MultiTargetFit.

Saving a Fit

import numpy as np

from nyx.core import dump_params, set_parameters

np.savez("fit.npz", **dump_params(fitted))

# to restore them onto a scene rebuilt the same way:
restored = set_parameters(scene, dict(np.load("fit.npz")))

Next Steps