Simulating observations

OITOOLS can generate synthetic OIFITS datasets from a parametric model or image, either by copying the UV coverage of an existing file or by building observations from scratch using array geometry and observation times.

From an existing OIFITS file

The simplest approach reuses the UV coverage and noise properties of a real dataset:

using OITOOLS
simulate_from_oifits("data/2004-data1.oifits", "data/sim.oifits";
                     image="data/2004true.fits", pixsize=0.101)

A flat-dict parametric model can be used instead of an image:

params = Dict("star,ud" => 3.0, "star,f" => 1.0)
model = dict_to_model(params, String[])

simulate_from_oifits("data/2004-data1.oifits", "data/sim.oifits";
                     flat_model=model, flat_params=Float64[])

See example_simulate_observations_from_OIFITS.jl.

From observation times and array geometry

To simulate a full night of observations at a specific interferometer, you need four configuration objects: facility, target, combiner, and wavelength setup. These are read from TOML files shipped with OITOOLS in src/configs/.

Configuration files

The .toml extension is optional — OITOOLS resolves built-in config names automatically.

Facility — array layout, telescope positions, atmospheric conditions:

facility = read_facility_file("CHARA")

Target — celestial coordinates and proper motion:

target = read_obs_file("default_obs")

You can also query SIMBAD directly:

ra, dec = ra_dec_from_simbad("Vega")   # decimal degrees, matching TargetConfig.raep0
target = TargetConfig(target="Vega", raep0=ra, decep0=dec)

simbad_target fills in the rest of TargetConfig from a single request, and carries the magnitudes simulate needs for its noise model:

t = simbad_target("Vega")
target = TargetConfig(target = t.main_id, raep0 = t.ra, decep0 = t.dec,
                      pmra = t.pmra, pmdec = t.pmdec,
                      parallax = t.plx, spectyp = t.sptype)

simulate(facility, target, combiner, wavelength, dates, "vega.oifits";
         mag = t.mags["H"], mag_ao = t.mags["V"])

Bands SIMBAD has no measurement for come back as NaN, never 0.0: zero is Vega-bright, and an unmeasured band written as zero would silently become the brightest in the row.

Combiner — beam combiner properties (throughput, read noise, calibration errors):

combiner = read_comb_file("MIRCX")

Wavelength — spectral setup for a given combiner mode:

wave = read_wave_file("MIRCX_LOWH")

Simulating from an image

using Dates

# Observation times: every 15 minutes over a 5.5-hour window
dates = collect(DateTime(2024,8,13,3,0,0):Minute(15):DateTime(2024,8,13,8,30,0))

facility = read_facility_file("CHARA")
target   = read_obs_file("default_obs")
combiner = read_comb_file("MIRCX")
wave     = read_wave_file("MIRCX_LOWH")

simulate(facility, target, combiner, wave, dates, "sim_image.oifits";
         image="data/2004true.fits", pixsize=0.101)

Simulating from a parametric model

params = Dict(
    "star,ud"    => 3.0,
    "star,f"     => 0.7,
    "disk,f"     => "1 - \$star,f",
    "disk,diamout" => 10.0,
    "disk,profile" => "exp(-(\$R/3.0)^2)",
)
model = dict_to_model(params, String[])

simulate(facility, target, combiner, wave, dates, "sim_model.oifits";
         flat_model=model, flat_params=Float64[])

See example_simulate_observations_from_model.jl and example_simulate_observations_from_image.jl.

Polychromatic simulation

example_simulate_polychromatic_disk.jl demonstrates simulating a chromatic, time-variable disk with an off-axis companion. The companion introduces wavelength-dependent differential phases, producing non-zero OIVIS signals. The example writes OIVIS2, OIVIS (with differential phases), OIT3, and OI_FLUX tables.

For an image cube, the spectrum matters: each plane is normalised to unit total flux (as it must be, for the visibilities to be correct), but the plane-to-plane totals are captured first and used to weight the photon count per channel and to fill OI_FLUX.

Keyword arguments to simulate

Coordinates follow the OIFITS standard: target.raep0 and target.decep0 are in degrees.

Noise model

Photons per telescope, per spectral channel, per frame:

N = F0(λ) · 10^(-m(λ)/2.5) · A_tel · δλ · DIT
    · T_atm(λ)             atmospheric transmission (`atm_transmission`, 1.0 by default)
    · facility.throughput  telescopes and beam train only
    · combiner.transmission   end-to-end: array + instrument, excluding QE/Strehl/atmosphere
    · flux_frac            split between interferometric and photometric channels
    · QE
    · S(λ, elevation, m_ao)   Strehl ratio

Matching ASPRO

The CHARA combiner configs (MIRCX.toml, MYSTIC.toml, SPICA.toml) are transcribed field for field from ASPRO 2's aspro-conf/…/CHARA.xml, and CombinerConfig mirrors ASPRO's FocalInstrumentSetup, so numbers can be copied across directly:

Because transmission is end-to-end, FacilityConfig.throughput is 1.0 for CHARA.

Two places where ASPRO cannot be reproduced exactly:

• SPICA's Strehl. ASPRO's NoiseService ignores the AO model for SPICA and hardcodes 0.25 / 0.15 / 0.10 by seeing, even though its own published Strehl plots use the AO model — ASPRO is internally inconsistent here. SPICA.toml sets strehl_model = "fixed_spica" to match ASPRO's noise; set strehl_model = "ao" for the physical model (~0.28 at 1″).
• SPICA's fringe tracker. ASPRO marks FT_SPICA as required, which permits much longer integrations than dit. Not modelled, so faint-end SPICA is pessimistic here.

S comes from strehl_ratio, a port of JMMC's Band.strehl, and reproduces ASPRO 2's published CHARA Strehl curves to a median 0.7%. It needs an [ao] block in the facility config; without one the code falls back to the seeing-limited coupling min(1,(r₀/D)²), which underestimates an AO-equipped array by roughly 5× in H and 20× in R.

Noise is drawn once, on the complex visibility, and every observable is derived from that one perturbation. VISAMP² == VIS2 exactly and the closure phase is exactly the sum of the three baseline phases — which is not true if each observable is noised independently.

T3AMP and T3PHI error bars are estimated by sampling (n_samples, default 100). The analytic closure errors are a small-error expansion and are only adequate while every baseline is well detected; measured reduced chi² against the true model, on a resolved disc:

Sampling costs no measurable time, so it is the default; set n_samples=0 for the analytic form.

Run demos/validate_noise_model.jl to print the Strehl comparison against ASPRO and the predicted σ(V²)/σ(CP) against magnitude for MIRC-X, MYSTIC and SPICA.

Observability filtering (opt-in)

simulate is a pure uv-coverage simulator: by default every epoch you pass is used, whether or not the target was above the horizon or within the delay lines. That is usually what you want when generating data to test image reconstruction.

To apply real constraints, either filter first:

dates_ok, mask, report = observable_epochs(facility, target, dates;
                                           min_elevation = 30.0,
                                           pops = [1,3,5,2,4,1])   # from best_pop
simulate(facility, target, combiner, wave, dates_ok, "sim.oifits"; flat_model=model)

or pass the same options through:

simulate(facility, target, combiner, wave, dates, "sim.oifits";
         flat_model=model, observability=(min_elevation=30.0,))

POP configurations are never chosen for you — omit pops and no delay-line check is done at all; run best_pop yourself if you want a recommendation.

Observation planning

OITOOLS provides tools for checking delay-line feasibility and producing Gantt charts for a given target and night:

The high-level entry point is obs_plan, which computes the night and renders the Gantt chart in one call:

facility = read_facility_file("CHARA")
ra, dec  = ra_dec_from_simbad("Vega")        # decimal degrees
config   = [1, 1, 1, 1, 1, 2]                # 0=unused, 1=use, 2=reference cart
pop      = [1, 1, 1, 1, 1, 1]                # one POP per telescope, 1:5

obs_plan("Vega", facility, ra, dec, DateTime(2026, 6, 3), pop, config;
         alt_limit = 30.0, savefile = "vega.png")

The pieces underneath, if you want them separately:

# Dark window, in decimal UT hours. `zenith` is in degrees: 102 is nautical twilight.
dusk_rise, dusk_set = sunrise_sunset(DateTime(2026, 6, 3), facility.lat, facility.lon)

# Everything for one night: LST, hour angle, altitude, azimuth, Moon separation
obs = night_observability(facility, ra, dec, DateTime(2026, 6, 3); alt_limit = 30.0)

# Altitude/azimuth directly — note the argument order (dec, lat, ha) and that `ha` is in hours
altitude, azimuth = alt_az(dec, facility.lat, obs.ha)

# Delay-line feasibility. An arbitrary POP choice often yields no usable time at all --
# that is what best_pop is for.
d = in_delay(facility, dec, obs.ha, config, pop)
results = best_pop(facility, dec, obs.ha, config; n_best = 5)
print_pop_results(facility, config, results)

See example_chara_plan.jl.