# Simulating zodi with a telescope's orbit

The `zodi_rate` and `zodi_readout` entry points in `coronagraphoto` compute the local zodiacal-light
count rate on the detector for one epoch and one line of sight. To
generate a time series of the kind a paper figure or a mission yield
calculation needs, the per-frame geometry has to be threaded in from an
observatory model. This page shows the recommended composition using
`orbix.observatory.ObservatoryL2Halo` for the observatory and
`skyscapes.background.LeinertZodi` for the brightness model, and
ends with a short validation checklist.

## The composition

A single-epoch zodi simulation pulls together four pieces, each
responsible for a distinct part of the calculation. The
`ObservatoryL2Halo` instance provides the heliocentric position of the
telescope and the per-target sky-geometry angles needed for the Leinert
lookup, the `LeinertZodi` instance evaluates the surface
brightness from those angles, the optical path threads the brightness
through the coronagraph's `sky_trans` map and the detector resampling,
and `zodi_rate` returns the per-pixel count rate. Adding
Poisson shot noise on top is the job of `zodi_readout`, which composes a
detector readout step around `zodi_rate`.

The two Leinert inputs are read off the geometry by the observatory
helpers as described in the skyscapes
[Local zodi + telescope geometry][skyscapes-zodi] doc, and then handed
to `zodi_rate`:


```python
from orbix.observatory import ObservatoryL2Halo
from skyscapes.background import LeinertZodi
from coronagraphoto import zodi_rate

obs = ObservatoryL2Halo.from_default()
zodi = LeinertZodi(reference_mag_arcsec2=22.0)

mjd = 60575.25
start_time_jd = mjd + 2_400_000.5
ecl_lat = float(obs.ecliptic_latitude_deg(mjd, ra_rad, dec_rad))
helio_lon = float(obs.helio_ecliptic_longitude_deg(mjd, ra_rad, dec_rad))

rate = zodi_rate(
    zodi,
    optical_path,
    start_time_jd=start_time_jd,
    wavelength_nm=550.0,
    bin_width_nm=50.0,
    ecliptic_lat_deg=ecl_lat,
    solar_lon_deg=helio_lon,
)
```

The returned `rate` is an `(ny, nx)` array of electrons per second on
the detector. Multiplying by an exposure time and adding Poisson shot
noise yields a realistic electron image, and `zodi_readout` packages
both steps:

## Year-long simulations

A year-long animation is a simple matter of looping the single-epoch
recipe over a sequence of MJD samples and stepping the geometry along
the L2 halo trajectory. Each frame draws an independent PRNG key so the
Poisson realisations are uncorrelated, and any frame whose Leinert
lookup falls in the out-of-range region returns `NaN`, which downstream
code uses as the "target unobservable this epoch" gate:

```python
import jax
import jax.numpy as jnp
import numpy as np
from coronagraphoto import zodi_readout

obs = ObservatoryL2Halo.from_default(equinox_mjd=60575.25)
zodi = LeinertZodi(reference_mag_arcsec2=22.0)
prng_keys = jax.random.split(jax.random.PRNGKey(0), n_frames)
mjds = 60575.25 + np.linspace(0.0, 365.25, n_frames)

ra_rad = jnp.deg2rad(target_ra_deg)
dec_rad = jnp.deg2rad(target_dec_deg)

for i, mjd in enumerate(mjds):
    el = float(obs.ecliptic_latitude_deg(float(mjd), ra_rad, dec_rad))
    sl = float(obs.helio_ecliptic_longitude_deg(float(mjd), ra_rad, dec_rad))
    image = zodi_readout(
        zodi,
        optical_path,
        prng_keys[i],
        start_time_jd=float(mjd) + 2_400_000.5,
        exposure_time_s=exposure_s,
        wavelength_nm=wavelength_nm,
        bin_width_nm=bin_width_nm,
        ecliptic_lat_deg=el,
        solar_lon_deg=sl,
    )
```

A useful pattern is to run a single target through a full year, with
the L2-halo trajectory and the telescope-to-target pointing vector
overlaid for visual context, and a four-target grid version for visual
validation against the expected phase relationships across the
ecliptic.

## What changes through the year

The spatial pattern of the zodi image on the detector is constant up to
a scalar, because `zodi_rate` multiplies a uniform sky
brightness by the coronagraph's `sky_trans` map and only the scalar
Leinert factor varies with epoch. The dark-hole region of the detector,
where `sky_trans` is approximately zero, stays dark all year, and the
off-axis regions where `sky_trans` is approximately one all brighten
and dim together as `Δλ_⊙` sweeps. Per-pixel shot noise scales as
`sqrt(rate * t_exp)` of whatever the scalar is at that epoch, so the
relative noise structure is also constant up to the same scalar.

This means the natural "annual quantity" to plot is the scalar
integrated count rate, not any per-pixel metric, because the per-pixel
contrast structure on the detector is fixed.

## Validation checklist

When wiring up a new survey, mission-yield, or paper figure, sanity
checks against the geometry catch the most common mistakes. The
recommended check picks three on-ecliptic targets at ecliptic
longitudes 0°, 90°, and 180°, which in equatorial coordinates are
`(RA=0°, Dec=0°)`, `(RA=90°, Dec=+23.44°)`, and `(RA=180°, Dec=0°)`.
Compensating for obliquity is what dictates the middle target's non-zero
declination, as explained in the
[Local zodi + telescope geometry][skyscapes-zodi] doc. A year-long `zodi_readout` loop over each of those targets and the
`argmax` day of the integrated count rate will show three peaks
separated by roughly 92 days, all of comparable amplitude, because each
target undergoes solar conjunction on a different date but with similar
zodi brightness at conjunction.

A fourth target placed at high ecliptic latitude, for example
`(RA=0°, Dec=+60°)`, serves as the flat-baseline control. The peak
integrated count rate is roughly 50 to 100 times smaller than the
on-ecliptic peaks, and no frame falls into the `NaN` unobservable
window because the line of sight never approaches the Sun.

## See also

The skyscapes [Local zodi + telescope geometry][skyscapes-zodi] doc
covers the geometry side of this pipeline.
{class}`orbix.observatory.ObservatoryL2Halo` documents the L2 halo
orbit interpolator and the sky-geometry helpers used above.
{class}`skyscapes.background.LeinertZodi` documents the
Leinert+1998 surface-brightness model.

[skyscapes-zodi]: ../../../skyscapes/docs/explanation/local_zodi_geometry.md