API reference

Generated from the source docstrings.

Atmosphere

pyturb.Atmosphere

A multi-layer turbulent atmosphere producing pupil OPD.

Parameters:

Name	Type	Description	Default
layers	sequence of Layer	The turbulence profile. cn2_fraction values are normalised to sum to 1 internally.	required
r0	float	Total Fried parameter [m] at wavelength and at zenith. Give either r0 or seeing.	None
seeing	float	Total seeing FWHM [arcsec] at wavelength and at zenith. Converted to r0 via the Kolmogorov relation.	None
wavelength	float	Reference wavelength [m] at which r0/seeing are defined. Default 500 nm.	5e-07
zenith_angle	float	Zenith angle [deg]. Scales r0 (cos^{3/5}) and layer ranges (sec) for the line of sight. Default 0. This is a scalar airmass/range approximation, not full slant-path geometry: the pupil grid and winds stay isotropic, so it does not apply the anisotropic pupil-to-tilted-layer coordinate transform a real slant path needs (a baseline in the zenith plane maps to a screen separation stretched by sec(z) versus the perpendicular one, making the structure function differ between axes by sec(z)^{5/3} — ~3.2x at 60 deg). Accurate as a seeing/range scaling near zenith; not a substitute for full slant-path geometry at large zenith angles.	0.0
diameter	float	Pupil diameter [m]; with n this sets pixel_scale = diameter/n. Default 8 m.	8.0
n	int	Pupil sampling in pixels across the diameter. Default 512.	512
L0	float	Override the outer scale [m] for every layer. Default: use each layer's own L0.	None
power_law	float	Power-law index of the phase PSD, PSD ~ f^{-power_law}. Default 11/3 (Kolmogorov/von Karman). Other values give non-Kolmogorov turbulence; r0 then acts as an amplitude knob rather than the strict Fried parameter. Applies to :meth:sample and to engine="spectral" :meth:frames/:meth:opd. Requires engine="spectral": the extruder's row recurrence is a closed-form conditional distribution derived specifically for power_law=11/3.	11.0 / 3.0
inner_scale	float	Inner scale [m]. 0 (default) disables it; a positive value rolls off the PSD above the corresponding spatial frequency (see :class:pyturb.PhaseScreen). Same engine restriction as power_law.	0.0
subharmonics	int	Subharmonic levels for low-frequency correction. Default 8.	8
field_of_view	float	Radius of the science/guide-star field [arcsec]. The generated screens are oversized so that off-axis footprints out to this radius sample genuinely different (non-wrapped) turbulence — set this whenever you use directions for anisoplanatism/tomography. Default 0 (screens exactly the pupil; fastest). Larger values cost a larger per-frame FFT.	0.0
tau_boil	float or sequence	Boiling (temporal decorrelation) time constant [s] of each layer's outer-scale structure, scalar or one per layer. None (default) is pure frozen flow. On the default spectral engine (boil_scaling="kolmogorov") smaller scales decorrelate faster (below), so this is the longest lifetime in the layer: with L0 = 25 m the 1-10 cm structure an AO system senses boils 40-180 times faster than tau_boil. Physically sensible values are then of order a second or more; a few tens of milliseconds makes the screen change by a large fraction of its fine structure every millisecond (see the Physics guide for the per-frame change). With boil_scaling="uniform" it is the lifetime of every scale; :func:pyturb.theory.travel_decorrelation_time gives an on-sky estimate per layer from its wind speed. When set, the layer relaxes toward fresh, statistically identical turbulence while preserving its spatial statistics, by the temporal process boil_model. Active only while stepping with :meth:frames/:meth:evolve. Works on both engines, but the two realise it differently: engine="spectral" boils per Fourier mode. With boil_scaling="kolmogorov" (default) finer spatial structure decorrelates faster than the outer scale per Kolmogorov eddy-turnover scaling, tau(f) = tau_boil * (f/f_ref)^(-2/3) for f >= f_ref = 1/L0 (the grid fundamental for Kolmogorov L0=inf), clamped to tau_boil below f_ref; with boil_scaling="uniform" every mode has tau_boil. engine="extrude" blends its ring buffer toward a fresh independently extruded screen -- a single-timescale AR(1) that decorrelates all spatial scales at tau_boil (real space has no per-mode handle). Staying non-periodic costs a modest deficit in the largest-scale power of the boiled screen, and re-extruding that fresh screen every step makes it markedly slower than spectral boiling (raises :class:ExtrudeBoilingPerformanceWarning); use "spectral" unless genuine non-periodicity is required. :func:pyturb.theory.covariance (or :meth:covariance) gives the space-time statistics of either.	None
boil_model	('ar1', 'ar2')	Temporal process of boiling. "ar1" (default): each mode is a first-order autoregression (an Ornstein-Uhlenbeck process) with retention exp(-dt/tau) per step and white innovations, so the screen's change over a short step grows as sqrt(dt/tau): it is rough in time and cannot be extrapolated from past frames. "ar2": a smooth, critically damped second-order process (a continuous-time AR(2), Matérn-3/2 correlation (1 + 2 lag/tau) exp(-2 lag/tau)), with the same integral time scale tau. Each mode carries a rate as well as its value, so the change over a short step grows as dt/tau, closer to real boiling, which is smooth at short lags. Sampled exactly for any dt. engine="spectral" only; it stores a second spectrum per layer (twice the boiling memory) and draws twice the noise per step.	"ar1"
boil_scaling	('kolmogorov', 'uniform')	Scale dependence of the boiling time constant. None (default) is the engine's own: "kolmogorov" on engine="spectral" and "uniform" on engine="extrude". "kolmogorov": the eddy-turnover tau(f) = tau_boil * max(f/f_ref, 1)^(-2/3) per mode, so fine structure boils faster than tau_boil. "uniform": every Fourier mode (main grid and subharmonics) has tau_boil, so the layer's correlation decays at one rate at every scale, as the extruder's does. Per-layer on-sky measurements (Guesalaga et al. 2014, MNRAS 440, 1925) describe the decay of a whole layer's correlation with the distance its wind carries it, not a per-scale law; pair "uniform" with :func:pyturb.theory.travel_decorrelation_time to simulate that. "kolmogorov" is rejected on engine="extrude", which has no per-mode handle. Works with both boil_model values.	"kolmogorov"
engine	('spectral', 'extrude')	Frozen-flow engine for :meth:frames/:meth:opd. "spectral" (default) is the fastest: an exact sub-pixel shift-theorem translation, all layers batched into one FFT — but the screen is periodic (it repeats after n*pixel_scale of travel). "extrude" is the Assémat–Wilson row extruder in a wind-aligned frame with rotated sub-pixel sampling: unbounded and non-periodic, the right choice for long closed-loop runs, at the cost of per-layer sampling instead of one batched FFT. Its sub-pixel interpolation is a position-dependent low-pass filter (exact at integer pixel travel, most attenuating at half-pixel travel), so its finest-scale (1-2 px) structure function deviates from theory by 5-15% and oscillates with the sub-pixel travel phase; "spectral" has no such artifact. (At larger separations the extruded field stays isotropic and on-theory to the method's usual few-percent accuracy — the row recurrence introduces no systematic large-scale along-vs-cross anisotropy.) Prefer "spectral" when fine-scale (near-Nyquist) fidelity matters more than non-periodicity. :meth:sample (Monte-Carlo) is unaffected by this choice.	"spectral"
interp	('cubic', 'linear', 'lanczos')	Sub-pixel interpolation kernel for engine="extrude". "cubic" (Catmull-Rom, default) is a good speed/quality balance; "lanczos" (6-tap Lanczos-3) has a flatter sub-Nyquist response that reduces the extruder's finest-scale structure-function deficit and its travel-phase flicker (see the engine note), at the cost of more taps per readout; "linear" (2-tap) is the lowest fidelity. Note that "cubic" and "lanczos" run through a fused CUDA/Numba readout kernel while "linear" falls back to a generic tap-broadcast gather that materialises temporaries, so despite its lower tap count "linear" is typically slower, not faster, than the default "cubic" on both CPU and GPU; pick it for its low-pass character, not for speed.	"cubic"
lgs_altitude	float	Altitude [m] of a laser guide star (e.g. 90e3 for sodium). When set, each layer's footprint is magnified by (1 - h/lgs_altitude) — the cone effect (focal anisoplanatism) that a finite-range beacon senses. None (default) is a natural-guide-star / science source at infinity. Works on both engines: the extruder samples its ring buffer on a magnified grid, and the spectral engine zoom-resamples each layer's screen about the pupil centre by the same factor (so on spectral the cone composes with tau_boil boiling). Must be greater than every layer altitude (a beacon at or below a layer is a degenerate cone).	None
dispersion	(None, 'edlen', 'ciddor')	Chromatic model for the OPD. None (default) treats the turbulence OPD as perfectly achromatic (phase at wavelength is just 2 pi * OPD / wavelength). "edlen" additionally scales the OPD by the dry-air refractivity ratio (n(wavelength)-1)/(n(lambda_ref)-1) before converting to phase — the small (~1-2 % visible-to-NIR) chromatic term that matters for high-contrast, astrometry and LGS error budgets. "ciddor" splits that scaling into a dry-air part and a water-vapour part weighted by wet_fraction (the "wet–dry" problem): because water vapour disperses differently from dry air, a dry-only correction is wrong in the mid-IR and for interferometry where the wet turbulence dominates. "ciddor" with wet_fraction=0 is identical to "edlen". Only affects output methods called with an explicit wavelength; the metre-valued OPD is unchanged. Both models use the Edlén (1966) dry-air and Ciddor (1996) water-vapour dispersion shapes, validated against the independent Ciddor standard-air formula to < 0.05 % over 0.35-1.7 µm (the dry-air ratio stays within ~0.1 % to ~2.5 µm). They are scalar models with no weather inputs: the wet/dry density split is the user-supplied wet_fraction, not a pressure/temperature/humidity computation, so treat the thermal-IR / interferometric wet-dominated regime as an order-of-magnitude estimate set by wet_fraction rather than a first-principles wet–dry model.	None
wet_fraction	float	Fraction (0–1) of the reference-wavelength turbulent refractivity carried by water vapour rather than dry air, used only by dispersion="ciddor". 0 (default) is pure dry air; typical near-surface visible/NIR values are small, but the wet term grows important in the thermal IR. Must be 0 unless dispersion="ciddor".	0.0
device	str	"cpu" (default) or "gpu".	'cpu'
dtype	str	"float32" (default) or "float64".	'float32'
seed	int	Master seed. Per-layer streams are spawned from it so results are reproducible and independent of layer count. A seed gives the same atmosphere on every CUDA GPU: the random draws are bit-identical, and the OPD agrees to float32 rounding (cuFFT/cuBLAS may pick different kernels on different GPU models). CPU (NumPy PCG64) and GPU (Philox) use different generators, so the same seed draws a different, statistically equivalent, realisation on each backend.	None
cuda_graph	bool	On the GPU, replay each spectral-engine frame (:meth:frames, :meth:evolve, single-direction :meth:opd, including the lgs_altitude cone) as a captured CUDA graph: one launch instead of a few dozen, which is what limits the frame rate on hosts with a modest CPU. Falls back to ordinary execution if the driver cannot capture. Default True; no effect on the CPU or on engine="extrude" frames.	True
oversample	float	Size of the FFT screens used by :meth:sample and the spectral engine, as a multiple of the pupil (>= 1, before the field_of_view margin). Default 1: screens exactly the pupil's size, the fastest choice. An FFT screen is periodic, so on a pupil-sized screen opposite pupil edges are neighbours: large-separation statistics come out low (structure function about -13% at D/2 and -40% edge to edge, tilt ~10-15% low, the two astigmatisms split ~0.6x/1.6x Noll; see the Validation docs) and new turbulence entering the pupil is what is leaving it. oversample=2 leaves ~-10% at r = D; 4 keeps statistics within ~2% out to D. The FFT grows as oversample**2: on an RTX 4060 at n=512, 9 layers, frames drop from ~5,300 fps to ~1,360 (x2) and ~300 (x4); at n=256, x2 is nearly free. It also multiplies :attr:time_to_wrap. No effect on engine="extrude" frames, which are non-periodic.	1.0

Feature compatibility

engine="spectral" (default) and engine="extrude" are not two interchangeable ways to get the same features faster/slower -- each supports a different subset:

• tau_boil (boiling) works on both engines, but with different spatial character: engine="spectral" boils per Fourier mode (by default finer structure faster, scale-resolved; boil_scaling= "uniform" gives every mode one rate), while engine="extrude" blends its ring buffer toward fresh extruded turbulence at a single timescale across all scales, preserving non-periodicity (see tau_boil above). boil_model="ar2" and boil_scaling="kolmogorov" need the spectral engine's per-mode state.
• non-Kolmogorov power_law/inner_scale and Kolmogorov L0=inf still require engine="spectral": they need either discrete Fourier modes or a closed-form covariance that only exists for the standard von Karman case, so none can combine with non-periodic (engine="extrude") evolution.
• lgs_altitude (the LGS cone effect) works on both engines. On the extruder the cone shrinks each layer's footprint via the ring-buffer sampling grid; on the spectral engine each layer's screen is zoom-resampled about the pupil centre by the same factor. On either engine it composes with tau_boil boiling (the cone acts on readout geometry, boiling on the stored turbulence, so they are independent). The spectral cone falls back to a per-layer transform (each layer zooms differently), so it runs below the on-axis spectral throughput.
• directions (off-axis/tomography) works with both engines, but every requested direction's radius sqrt(thx**2 + thy**2) must lie within the declared field_of_view (a ValueError is raised otherwise, in both engines) -- construct the Atmosphere with a field_of_view covering every direction you plan to request.

Examples:

>>> import pyturb
>>> atm = pyturb.Atmosphere.from_profile("two-layer", seeing=0.8,
...                                      diameter=8.0, n=256, seed=1)
>>> opd = atm.opd()                    # (256, 256) OPD [m], t = 0
>>> for t, frame in atm.frames(dt=1e-3, steps=10):
...     pass                           # closed-loop OPD [m]
>>> ensemble = atm.sample(16)          # (16, 256, 256) independent OPDs

r0 property

r0

Line-of-sight total Fried parameter [m] at the reference wavelength.

seeing property

seeing

Line-of-sight seeing FWHM [arcsec] at the reference wavelength.

The Kolmogorov value 0.98 lambda / r0. A finite outer scale makes the delivered long-exposure FWHM smaller (~17% at r0 = 0.1 m, L0 = 25 m); see :func:pyturb.theory.seeing_fwhm.

theta0 property

theta0

Isoplanatic angle [arcsec] at the reference wavelength.

tau0 property

tau0

Coherence time [s] at the reference wavelength.

greenwood_frequency property

greenwood_frequency

Greenwood frequency [Hz] at the reference wavelength.

time_to_wrap property

time_to_wrap

Seconds of on-axis wind travel before the fastest layer wraps.

Only meaningful for engine="spectral", which is periodic with period n_screen * pixel_scale metres of travel per layer; a run (via :meth:frames or :meth:opd) whose duration exceeds this re-samples the same screen realisation instead of independent turbulence for at least one layer. Returns inf for engine="extrude" (non-periodic) or if every layer has zero wind.

time property

time

Current internal clock [s] (advanced by :meth:frames).

metadata property

metadata

Provenance describing this atmosphere, for saving with output.

A flat dict of scalars/strings suitable for :func:pyturb.save headers — geometry, line-of-sight r0, the integrated seeing/theta0/tau0, and the main construction parameters — so a saved OPD records how it was made. For atmospheres built via :meth:from_profile it also records the profile name and its :func:pyturb.profile_info provenance (source, whether it is traceable to a published table, site, and the representativeness caveat). This is a descriptive summary plus, under "config", the JSON of :meth:to_config (every constructor input and the layer table), from which :meth:from_config rebuilds the atmosphere; evolved state (time, boiling) is not stored, so replay a run by stepping it the same way. L0/tau_boil are reported only when every layer shares one value (otherwise None, which :func:pyturb.save drops).

from_profile classmethod

from_profile(name, wind=None, wind_direction=None, **kwargs)

Build from a named profile (see :func:pyturb.list_profiles).

Any :class:Atmosphere keyword may be passed, e.g.::

Atmosphere.from_profile("paranal-median", seeing=0.8,
                        zenith_angle=30, diameter=8, n=512)

If neither r0 nor seeing is given and the profile's source publishes its strength (profile_info(name).conditions["r0"], the traceable site profiles), that r0 is used.

wind replaces the layers' wind speeds ("bufton", a scalar, or one per layer; see :func:pyturb.with_wind) and wind_direction their directions. Profiles whose source publishes no winds have static (0 m/s) layers and warn unless wind is given.

The profile name and its :func:pyturb.profile_info provenance (traceable-vs-representative, source, site) are recorded in :attr:metadata so a saved OPD carries where its atmosphere came from.

from_cn2 classmethod

from_cn2(heights, cn2, n_layers=10, wind='bufton', wind_direction=0.0, method='equivalent', layer_L0=25.0, source=None, **kwargs)

Build from a continuous or measured :math:C_n^2(h) profile.

The profile (heights [m] above the telescope, at zenith, and cn2 [m^-2/3] on that grid — e.g. a MASS-DIMM, SCIDAR or model profile) is compressed with :func:pyturb.discretize_cn2 (n_layers, wind, wind_direction, method; every layer gets outer scale layer_L0), then passed to the constructor with kwargs. If neither r0 nor seeing is given, r0 is the profile's own, (0.423 k^2 integral Cn2 dh)^(-3/5) at wavelength (default 500 nm). source is recorded in :attr:metadata.

import numpy as np, pyturb h = np.geomspace(10, 25e3, 2000) atm = pyturb.Atmosphere.from_cn2(h, pyturb.hufnagel_valley(h), ... n_layers=8, n=64)

to_config

to_config()

Every constructor input as a JSON-serialisable dict (see :meth:from_config).

Includes the full layer table (altitudes, normalised cn2_fraction, winds, per-layer L0), so Atmosphere.from_config(atm.to_config()) rebuilds the same atmosphere; with an integer seed it reproduces the same realisation and frames. Evolved state (time, boiling) is not included: replay by stepping the rebuilt atmosphere the same way. :attr:metadata carries this as a JSON string under "config", so a saved OPD records how to regenerate it.

from_config classmethod

from_config(config, **overrides)

Rebuild an atmosphere from :meth:to_config output (a dict or JSON string).

overrides replace individual inputs, e.g. device="gpu" to replay a CPU run on the GPU (the same seed draws a different realisation on a different backend, but the same one on any GPU).

Version-1 configs (written by pyturb 1.x, whose wind_direction was measured from axis 0 toward axis 1) are read with each layer's direction converted to (90 - wind_direction) mod 360, the same wind in the current frame (from +x = axis 1 toward +y = axis 0), so they replay the frames they recorded.

r0_at

r0_at(wavelength)

Line-of-sight total Fried parameter [m] at wavelength.

covariance

covariance(points_a, points_b, lag=0.0, wavelength=None)

Theoretical space-time covariance of this atmosphere's on-axis output.

<opd(a, t) opd(b, t + lag)> [m^2] between pupil points points_a and points_b (each one (x, y) or an (N, 2) array [m], as in :meth:opd_at), from :func:pyturb.theory.covariance with this atmosphere's layers, line-of-sight r0, per-layer L0, winds and boiling: tau_boil, boil_model and boil_scaling (scale-dependent by default on engine="spectral", single-timescale with "uniform" and always on engine="extrude"). With wavelength, phase [rad^2] at that wavelength (including dispersion). Shape (Na, Nb), or (T, Na, Nb) for a 1-D lag [s].

It is the infinite-screen statistic. The simulated screens lack power above the grid's Nyquist frequency and, on the periodic spectral engine, some of the largest scales (oversample reduces that); differences such as :meth:structure_function are much less sensitive to the large scales. Needs a finite L0 (the Kolmogorov variance is infinite).

structure_function

structure_function(points_a, points_b, lag=0.0, wavelength=None)

Theoretical structure function <[opd(a, t) - opd(b, t + lag)]^2> [m^2].

From :func:pyturb.theory.space_time_structure_function, accurate for close points and short lags. With points_a == points_b it is the mean-square change of the screen at one point over lag: for lag = dt, the per-frame change :meth:frames produces. Finite for Kolmogorov layers (L0=inf) unless they boil.

scintillation_index

scintillation_index(wavelength=None)

Theoretical scintillation index of a star seen through this atmosphere.

:func:pyturb.theory.scintillation_index (Roddier 1981, weak perturbation, point aperture) with this atmosphere's layers, zenith r0, reference wavelength and zenith angle, at wavelength [m] (default: the reference wavelength), scaled by the square of the dispersion model's OPD ratio. The reference for the variance of :meth:pupil_field's intensity, which a simulated grid reaches only when its pixel is well below the Fresnel scale (the power above its Nyquist frequency is missing).

opd

opd(t=0.0, directions=None, wavelength=None)

OPD at time t [s] under frozen flow.

Parameters:

Name	Type	Description	Default
t	float or 1-D array_like	Time since the start of the simulation [s]. An array of times returns a leading time axis ((len(t), n, n), or (len(t), len(directions), n, n)); on the GPU the spectral engine evaluates many times in one batched transform, the fastest way to generate an offline time series. Boiling is not applied to random-access times (use :meth:frames). Taylor frozen flow: each layer's phase evolves as phi(r, t) = phi_0(r + wind_vector * t) with r = (x, y) — a fixed pupil point sees the turbulence that was wind_vector * t metres further along the wind direction at t=0 (the pattern is carried past the aperture by the wind, not translated bodily along +wind_vector in pupil coordinates). engine="spectral" (default) supports arbitrary random-access t; engine="extrude" is streaming: t must not decrease relative to a previous opd/frames/ evolve call on the same object (a decreasing t raises ValueError), and :meth:reset restarts it at t=0.	0.0
directions	sequence of (thx, thy)	Off-axis directions [arcsec] from the on-axis line of sight, each (thx, thy) or (thx, thy, altitude), with thx along x (axis 1, columns) and thy along y (axis 0, rows). altitude [m, at zenith] gives that direction its own source: a laser guide star at that range (cone effect), or None for a natural star or science target at infinity, so an LGS asterism, NGS and science directions can share one call and one turbulence realisation. Two-element entries use the constructor's lgs_altitude. On engine="extrude" a source's footprint may not exceed the one the buffer was built for: build with lgs_altitude=None and give LGS directions their altitude to mix them. Each direction's radius sqrt(thx**2 + thy**2) must lie within the field_of_view declared at construction, or the screens are not guaranteed to be oversized enough and the request raises ValueError instead of silently sampling stale/wrapped data. Each layer's footprint is shifted by altitude_los * tan(theta) along the direction's own axis (thx > 0 moves it toward increasing column index) for anisoplanatism / tomography studies. If given, the result has a leading axis of length len(directions).	None
wavelength	float	If given, return phase [rad] at this wavelength; otherwise OPD [m].	None

opd_at

opd_at(x, y, t=0.0, direction=(0.0, 0.0), altitude=_DEFAULT_SOURCE, wavelength=None)

OPD at arbitrary pupil-plane coordinates (x, y) [m].

x/y (broadcastable arrays, any shape) are offsets from the pupil centre along x (axis 1, columns) and y (axis 0, rows) — e.g. DM actuator or sub-aperture positions, a sparse or segmented aperture, or several apertures. Each layer is read at its frozen-flow, off-axis-shifted and cone-magnified position with the interp kernel and the layers summed, so on the pupil grid itself this reproduces :meth:opd (exactly at its pixel centres). direction [arcsec] and altitude (None for a source at infinity; default the constructor's lgs_altitude) are as in :meth:opd. Points must lie on the simulated screen: on engine="spectral" the FFT screen (widen it with oversample or field_of_view to reach beyond the pupil, e.g. for long baselines); on engine="extrude" the extruded strip (the pupil footprint plus the field_of_view margin). Boiling state is as for :meth:opd.

pupil_field

pupil_field(t=None, direction=(0.0, 0.0), wavelength=None, weights=None)

A star's light at the pupil after Fresnel propagation through the layers.

A unit plane wave from direction enters above the highest layer, takes each layer's phase and propagates down to the next layer, and from the lowest layer to the telescope, with the paraxial angular spectrum (:func:pyturb.scintillation.propagate_layers, built on :class:aocore.AngularSpectrumPropagator), at each wavelength. The layers are this atmosphere's screens at time t (frozen flow, and the current boiling state, as for :meth:opd), at their line-of-sight ranges (altitude * sec(zenith_angle)). Diffraction between the layers turns their phase into intensity fluctuations: scintillation.

Conventions: the field takes the phase exp(+i 2 pi OPD / lambda) and propagates over z with exp(-i pi lambda z f^2), the same exp(+i k z) convention as aocore (an OPD ramp a x images at +a). An off-axis direction is handled exactly in the paraxial approximation: the field of a plane wave tilted by theta is the tilt times the on-axis propagation through the layers shifted by theta * range, the footprints :meth:opd uses for directions. The returned field is that on-axis part, relative to the unperturbed tilted wave (a constant phase is dropped).

Parameters:

Name	Type	Description	Default
t	float	Time [s]. Default: the current clock, :attr:time, so inside a :meth:frames loop or after :meth:evolve/:meth:advance the field matches the latest frame.	None
direction	(thx, thy)	Star direction [arcsec] from the line of sight, thx along x (axis 1) and thy along y (axis 0); its radius must lie within field_of_view, as for :meth:opd.	(0.0, 0.0)
wavelength	float or sequence of float	Wavelength(s) [m] to propagate. Default: the reference wavelength the atmosphere was built with. A sequence gives a leading wavelength axis on field. With dispersion the OPD at each wavelength is scaled as in :meth:opd.	None
weights	sequence of float	Weight of each wavelength in intensity (e.g. the photon spectrum over the band), normalised to sum to 1. Default equal.	None

Returns:

Type	Description
PupilField	(field, intensity, opd): the complex field(s), the weighted relative intensity and the summed OPD [m], each cropped to the (n, n) pupil and on the atmosphere's device.

Notes

Each layer is a thin screen and the propagation runs on the periodic spectral screen (n_screen pixels, see oversample), so light leaving one edge of the screen enters at the other, as the periodic turbulence implies. The screens' subharmonics (scales larger than the screen) are applied as phase at the pupil: over atmospheric ranges they do not diffract (their Fresnel number L^2 / (lambda z) is in the thousands), and the shift they give the pattern below a layer (range times their slope, typically millimetres) is left out. The grid must resolve the Fresnel scale sqrt(lambda z) of the layers and the intensity structure: a coarse grid misses the scintillation power above its Nyquist frequency. engine="spectral" only (the extruder's screens are not periodic), and only for a source at infinity (an LGS is a spherical wave). The CPU transforms run through aocore's backend, with AOCORE_FFT_WORKERS threads, not :func:set_threads.

Examples:

>>> import pyturb
>>> atm = pyturb.Atmosphere.from_profile("two-layer", r0=0.15, n=64,
...                                      diameter=2.0, seed=1)
>>> pf = atm.pupil_field(wavelength=[0.6e-6, 0.8e-6], weights=[1.0, 2.0])
>>> pf.field.shape, pf.intensity.shape
((2, 64, 64), (64, 64))

frames

frames(dt, steps, wavelength=None)

Yield (t, opd) for steps frames spaced by dt seconds.

Advances the internal clock, so successive calls continue the wind. Use :meth:reset to return to t = 0.

Yields:

Type	Description
(float, ndarray)	Simulation time [s] and the pupil OPD [m] (or phase [rad] if wavelength is given), shape (n, n) on the device.

evolve

evolve(dt, wavelength=None)

Advance the wind by dt seconds and return the new pupil OPD.

A single-step, in-seconds convenience over :meth:frames for callers driving their own loop (mirrors HCIPy's evolve_until): each layer's frozen flow is stepped by its own wind_speed — you pass a time, not a pixel count. Advances the internal clock, so successive calls continue the wind; use :meth:reset to return to t = 0.

Repeated evolve(dt) calls reproduce the same sequence as frames(dt, ...) (from its second frame on): both apply one boiling step per dt when tau_boil is set.

Parameters:

Name	Type	Description	Default
dt	float	Time step [s] (> 0).	required
wavelength	float	If given, return phase [rad] at this wavelength; otherwise OPD [m].	None

Returns:

Type	Description
ndarray	The pupil OPD [m] (or phase [rad]) at the new time, shape (n, n) on the device.

advance

advance(dt)

Advance the clock by dt [s] in one exact boiling step; returns self.

:meth:evolve without computing an OPD: read the new state with :meth:opd or :meth:pupil_field at t = atm.time. Frozen flow needs no stepping (each layer is evaluated at any time directly), and boiling is a continuous-time process that every engine samples exactly at any step: with boil_model="ar1" each mode (each layer on the extruder) keeps exp(-dt / tau) of its value and gains an independent innovation of variance 1 - exp(-2 dt / tau), and "ar2" uses its exact transition over dt. One step of s * T therefore has the statistics of s steps of T: an AR(1) sampled every s steps is an AR(1) with coefficient rho**s (rho = exp(-T / tau)) and innovation variance 1 - rho**(2 s). A history of frames tens of milliseconds apart costs one step per frame, not one per elapsed loop period. The random draws differ from those of s separate steps, so the realisation differs; :meth:frames and :meth:evolve take one step per dt.

Parameters:

Name	Type	Description	Default
dt	float	Time step [s] (> 0).	required

Examples:

>>> import pyturb
>>> atm = pyturb.Atmosphere.from_profile("two-layer", r0=0.15, n=64,
...                                      tau_boil=2.0, seed=1)
>>> history = []
>>> for step in (0.0, 0.02, 0.02, 0.03):     # frames at 0, 20, 40, 70 ms
...     if step:
...         atm.advance(step)
...     history.append(atm.opd(t=atm.time))

sample

sample(count=None, wavelength=None)

Draw statistically independent integrated OPDs.

Each call produces fresh, uncorrelated realisations of the summed atmosphere (on-axis), so the ensemble structure function matches the total r0. Layers that share an outer scale are drawn as a single aggregate screen (independent von Kármán screens with the same PSD shape add exactly, r0_agg^{-5/3} = sum_i r0_i^{-5/3}), so this costs one FFT per distinct L0 rather than one per layer -- for a profile whose layers share L0 it is ~L times faster than a per-layer sum, and distributionally identical.

Parameters:

Name	Type	Description	Default
count	int	Number of independent OPDs. If omitted a single (n, n) array is returned; otherwise the result is (count, n, n).	None
wavelength	float	If given, return phase [rad] at this wavelength; otherwise OPD [m].	None

reset

reset()

Reset to t = 0 and restore the initial turbulence. Returns self.

Rewinds the internal clock and, when boiling (tau_boil) has run, restores the pre-boil turbulence and boiling RNG so a reused atmosphere replays identically. For engine="extrude" the extruded layers are rebuilt from their seeds (wind travel is monotonic, so the run restarts identically). For engine="spectral" frozen flow never mutates the stored spectra, but boiling does, so the stacked spectra, subharmonic coefficients and boiling RNG are rebuilt from the layers' fixed realisations.

Warnings

pyturb.PeriodicWrapWarning

Bases: UserWarning

The default engine="spectral" screen has wrapped (repeated).

The spectral engine is periodic with period n_screen * pixel_scale metres of wind travel per layer. Once a layer's cumulative travel passes that period, its "new" turbulence is the same realisation seen before, biasing any temporal statistic computed over the run (PSDs, long integrations, closed-loop error budgets). Silence with warnings.filterwarnings("ignore", category=pyturb.PeriodicWrapWarning) if this is intentional (e.g. deliberately studying the periodic case), or switch to engine="extrude" for genuinely non-periodic frozen flow.

pyturb.ExtrudeBoilingPerformanceWarning

Bases: UserWarning

tau_boil on engine="extrude" is markedly slower than "spectral".

The spectral engine boils with one per-mode AR(1) update on its already- stored spectrum. The extruder has no such spectrum -- it works in real space -- so each boiling step must re-extrude an independent, fresh (rows, width) screen from scratch (seeding a few stencil rows, then running the same row recurrence forward) and blend it into the ring buffer. That extra extrusion, repeated every :meth:Atmosphere.frames / :meth:Atmosphere.evolve step, is the dominant cost: a boiling extrude frame is many times a frozen one, far more than the corresponding slowdown on the spectral engine. Prefer engine="spectral" for boiling-heavy runs unless genuine non-periodicity is required; silence with warnings.filterwarnings("ignore", category=pyturb.ExtrudeBoilingPerformanceWarning) if the cost is already accounted for.

Scintillation

pyturb.PupilField

Bases: NamedTuple

A star's light at the telescope pupil (:meth:pyturb.Atmosphere.pupil_field).

Attributes:

Name	Type	Description
field	ndarray	Complex pupil field relative to the unperturbed plane wave from the star's direction (unit amplitude without turbulence), (n, n) for one wavelength or (W, n, n) for W wavelengths.
intensity	ndarray	(n, n) relative intensity sum_w weight_w |field_w|^2 (weights normalised to sum to 1): its mean is close to 1 and its variance is the scintillation index.
opd	ndarray	(n, n) summed OPD of the layers [m], the same map as :meth:pyturb.Atmosphere.opd at that time and direction (to rounding).

pyturb.scintillation.propagate_layers

propagate_layers(opd, altitudes, pixel_scale, wavelength)

Propagate a unit plane wave through layer OPD maps down to the pupil.

The field enters above the highest layer, takes each layer's phase exp(+i 2 pi opd[l] / lambda) and propagates to the next layer below (and from the lowest layer to the pupil) with the paraxial angular-spectrum transfer function exp(-i pi lambda z f^2) (:class:aocore.AngularSpectrumPropagator, paraxial=True). The grid is periodic: each map is one period of a periodic screen, as the FFT screens of engine="spectral" are. Energy is conserved: the mean of |field|^2 over the grid is 1 to rounding.

Parameters:

Name	Type	Description	Default
opd	ndarray	(L, ny, nx) OPD [m] of each layer (NumPy or CuPy). float32 input propagates in single precision, anything else in double; the result stays on the input's device.	required
altitudes	sequence of float	Range of each layer from the pupil [m] along the line of sight (altitude times sec(zenith angle)), >= 0. Layers at one range are applied in the given order (their phases commute).	required
pixel_scale	float	Grid pitch [m].	required
wavelength	float or sequence of float	Wavelength(s) [m]. A scalar returns (ny, nx), a sequence (W, ny, nx).	required

Returns:

Type	Description
ndarray	Complex pupil field(s), relative to the unperturbed plane wave.

Notes

The CPU transforms run through aocore's backend, whose thread count is AOCORE_FFT_WORKERS (not :func:pyturb.set_threads).

Examples:

>>> import numpy as np, pyturb
>>> opd = 1e-7 * np.random.default_rng(0).standard_normal((2, 64, 64))
>>> field = pyturb.scintillation.propagate_layers(opd, [8e3, 0.0], 0.02, 1.6e-6)
>>> field.shape
(64, 64)

Phase screens

pyturb.PhaseScreen

Generator of statistically independent atmospheric phase screens.

Each call to :meth:generate returns fresh, uncorrelated realisations of von Kármán (finite L0) or Kolmogorov (L0=inf) turbulence. For temporally correlated, frozen-flow screens see :class:pyturb.InfinitePhaseScreen.

Parameters:

Name	Type	Description	Default
n	int	Screen size in pixels (screens are n x n).	required
pixel_scale	float	Size of one pixel in metres.	required
r0	float	Fried parameter in metres, defined at the wavelength at which the phase is expressed (the returned phase is in radians at that wavelength). Use :func:pyturb.r0_at_wavelength to rescale.	required
L0	float	Outer scale in metres. numpy.inf gives pure Kolmogorov statistics. Default 25 m.	25.0
subharmonics	int	Number of subharmonic levels used to restore low spatial-frequency power that the periodic FFT grid cannot represent. 0 disables. Each level extends spectral coverage by a further factor of 3 in scale; the default of 8 reproduces the Kolmogorov structure function to within a few percent at separations of a few pixels and up. The 1-2 pixel (near-Nyquist) scales show a larger deficit (order 5-10%) that comes from finite grid sampling, not from subharmonic depth, so quote the "few percent" figure with a >~2-4 pixel resolution qualifier at the shortest scale of interest.	8
power_law	float	Power-law index of the phase PSD, PSD ~ f^{-power_law}. Default 11/3 (Kolmogorov, structure function ~ r^{5/3}). Other values give non-Kolmogorov turbulence (D(r) ~ r^{power_law-2}); r0 then acts as an amplitude knob rather than the strict Fried parameter.	11.0 / 3.0
inner_scale	float	Inner scale l0 [m]. 0 (default) disables it; a positive value rolls the PSD off above fm = 5.92/(2 pi l0) (modified von Kármán / Hill), flattening the structure function at small separations.	0.0
seed	int	Seed for the random generator. A seed reproduces the same screens for a given backend and dtype. On the CPU the stream is numpy.random.default_rng(seed) (PCG64). On the GPU it is pyturb's counter-based Philox4x32-10 stream, which is identical on every CUDA device. The backends use different generators, so the same seed draws a different (statistically equivalent) realisation on CPU and GPU.	None
device	str	"cpu" (default) or "gpu" (requires CuPy).	'cpu'
dtype	str or dtype	Floating dtype of the output, default "float32" (recommended for GPU work). Use "float64" for maximum accuracy.	'float32'

Examples:

>>> import pyturb
>>> gen = pyturb.PhaseScreen(n=256, pixel_scale=0.02, r0=0.15, seed=0)
>>> phase = gen.generate()          # (256, 256) radians
>>> batch = gen.generate(32)        # (32, 256, 256), one FFT batch

generate

generate(count=None)

Generate independent phase screens.

Parameters:

Name	Type	Description	Default
count	int	Number of screens. If omitted, a single (n, n) array is returned; otherwise the result has shape (count, n, n).	None

Returns:

Type	Description
ndarray	Phase in radians, on the selected device (numpy array on CPU, cupy array on GPU).

pyturb.InfinitePhaseScreen

Endless frozen-flow phase screen, extruded row by row.

The current (n, n) screen is available as :attr:screen. Advance the wind either by whole pixels with :meth:step or by any continuous distance (in pixels) with :meth:advance; new turbulence is synthesised at the leading edge (screen[-1]) as needed and older rows are recycled, so memory stays bounded no matter how long the run. The pattern therefore moves along y (axis 0, rows) toward decreasing row index; transpose the screen for motion along x (axis 1).

Translate pixels into wind speed via your loop rate: a step of wind_speed * dt / pixel_scale pixels advances the screen by v*dt metres. :meth:advance accepts the fractional result directly.

Parameters:

Name	Type	Description	Default
n	int	Screen size in pixels.	required
pixel_scale	float	Pixel size in metres.	required
r0	float	Fried parameter in metres (at the wavelength of the returned phase).	required
L0	float	Outer scale in metres. Must be finite (the conditional-covariance construction is von Kármán); default 25 m.	25.0
stencil_rows	int	Number of edge rows the new row is conditioned on. Default 2 (per Assémat & Wilson; more rows cost setup time for marginal gain).	2
interp	('cubic', 'linear', 'lanczos')	Sub-pixel interpolation kernel used by :meth:advance. "cubic" (Catmull-Rom, default) preserves high-frequency power better than "linear"; "lanczos" (6-tap Lanczos-3) is flatter still just below Nyquist, cutting the finest-scale structure-function deficit and the travel-phase flicker further at the cost of more taps. All three are exact at integer offsets. Unused by :meth:step.	"cubic"
seed	Any	As for :class:pyturb.PhaseScreen.	None
device	Any	As for :class:pyturb.PhaseScreen.	None
dtype	Any	As for :class:pyturb.PhaseScreen.	None

Examples:

>>> import pyturb
>>> layer = pyturb.InfinitePhaseScreen(n=128, pixel_scale=0.05,
...                                    r0=0.15, L0=25, seed=0)
>>> for _ in range(100):
...     phase = layer.step()        # advance wind by one whole pixel
>>> phase = layer.advance(0.37)     # ...and by 0.37 of a pixel (sub-pixel)

screen property

screen

Current (n, n) phase screen in radians (device array).

travel property

travel

Total wind travel so far, in pixels (metres = travel * pixel_scale).

step

step(steps=1)

Advance the wind by steps whole pixels and return the screen.

Each step shifts the screen one row along axis 0 and extrudes a new statistically consistent row at screen[-1]. For fractional (wind v*dt) motion use :meth:advance.

advance

advance(pixels)

Advance the wind by pixels (any non-negative float) and return it.

The pupil is interpolated at the exact sub-pixel offset (see interp); integer offsets are reproduced exactly. Successive calls accumulate, so advance(0.5) twice lands on the same screen as step(1).

pyturb.FourierFlowScreen

One frozen (Taylor) layer evaluated at arbitrary continuous offsets.

Parameters:

Name	Type	Description	Default
template	PhaseScreen	A configured generator supplying the PSD amplitude filters, grid, backend and dtype. The layer inherits n, pixel_scale, r0, L0 and device from it.	required
seed	int	Seed for the fixed coefficient realisation. Reuse the same seed to reproduce a layer exactly on the same backend (see seed in :class:pyturb.PhaseScreen for the CPU/GPU streams).	None

Examples:

>>> import pyturb
>>> from pyturb.flow import FourierFlowScreen
>>> template = pyturb.PhaseScreen(n=256, pixel_scale=0.02, r0=0.15)
>>> layer = FourierFlowScreen(template, seed=0)
>>> a = layer.translate(0.0, 0.0)          # (256, 256) radians
>>> b = layer.translate(1.3, -0.4)         # 1.3 m along x, -0.4 m along y

reseed

reseed(seed=None)

Draw a fresh fixed realisation of the layer's Fourier coefficients.

With no argument, advances the existing random stream; pass seed to restart it. Use this to get a statistically independent layer while keeping the same PSD/grid configuration.

translate

translate(sx, sy)

Return the screen blown by (sx, sy) metres, shape (n, n).

sx is displacement along x (axis 1, columns), sy along y (axis 0, rows): the returned screen is phi_0(r + s), so a positive sx moves the pattern toward decreasing column index. Values are phase in radians at the layer's reference wavelength; the array lives on the layer's device.

pyturb.phase_covariance

phase_covariance(r, r0, L0)

Von Kármán phase covariance :math:C_\phi(r) in rad^2.

Parameters:

Name	Type	Description	Default
r	array_like	Separation(s) in metres.	required
r0	float	Fried parameter in metres.	required
L0	float	Outer scale in metres (must be finite).	required

Notes

C(0) = 0.0863 (L0/r0)^(5/3) is the total von Kármán phase variance.

Profiles

pyturb.profiles.Layer dataclass

One turbulent layer.

Parameters:

Name	Type	Description	Default
altitude	float	Height above the telescope [m], at zenith.	required
cn2_fraction	float	Fraction of the total integrated turbulence (:math:C_n^2\,dh) in this layer. Fractions across a profile are normalised to sum to 1 by :class:pyturb.Atmosphere; relative values are what matter. A layer with fraction 0 is allowed and contributes nothing (Atmosphere leaves it out of the simulated layers).	required
wind_speed	float	Wind speed [m/s].	10.0
wind_direction	float	Direction the wind blows from [deg] (the meteorological convention), measured from +x (axis 1, columns) toward +y (axis 0, rows). Arbitrary (need not be axis-aligned). The turbulence pattern therefore moves across the pupil along -wind_vector: with wind_direction=0 it travels toward decreasing column index and new turbulence enters at the high-column edge; with 90 it travels toward decreasing row index. Every engine uses this convention.	0.0
L0	float	Outer scale [m] for this layer. Default 25 m.	25.0

wind_vector property

wind_vector

(vx, vy) [m/s] along x (axis 1) and y (axis 0), pointing upwind.

The vector points where the wind comes from, and the frozen-flow pattern moves along -wind_vector: the phase at pupil position r = (x, y) and time t is phi(r, t) = phi_0(r + wind_vector * t).

pyturb.get_profile

get_profile(name)

Return a fresh list of :class:Layer for a named profile.

Names: "single-layer", "two-layer", "paranal-median", "mauna-kea", "keck", "las-campanas", "cerro-pachon", "armazones", "hv57", "paranal-p01" ... "paranal-p14", and the traceable site profiles "tmt-<site>-<good|typical|bad>", "cerro-pachon-<class>", "siding-spring-gl-<class>-fa-<class>", "sutherland-median" and "mt-graham-<class>" (sources in :func:profile_info). See :func:list_profiles.

"mauna-kea", "keck", "las-campanas" and the fourteen "paranal-pNN" profiles (Garcia-Rissmann et al. 2015: Paranal seeing classes 0.4"-1.4", each "good", "median" or "bad", with their published r0, tau0 and probability in profile_info(name).conditions) are traceable to a specific published table (see the source comments for the citations); "paranal-median", "cerro-pachon", "armazones", "single-layer" and "two-layer" are representative/illustrative rather than a specific cited site survey. Wind direction is illustrative in every profile — none of the cited sources tabulate it. Call :func:profile_info for the machine-readable provenance (source, traceable flag, site, caveat) of any named profile.

pyturb.list_profiles

list_profiles()

Names accepted by :func:get_profile and Atmosphere.from_profile.

pyturb.profile_info

profile_info(name)

Return the :class:ProfileInfo provenance record for a named profile.

profile_info("mauna-kea").traceable True profile_info("paranal-median").traceable False

pyturb.ProfileInfo dataclass

Machine-readable provenance for a named profile.

The prose in each profile-building function's comment, made queryable and serialisable so a saved OPD can carry where its atmosphere came from and how much to trust it.

Attributes:

Name	Type	Description
name	str	The profile key (as passed to :func:get_profile).
traceable	bool	True if the Cn2/altitude layers reproduce a specific published table (or an analytic model); False for a representative discretisation that only follows the general shape of a site's statistics.
source	str	Citation for a traceable profile, or a short description of what a representative one is modelled on.
site	str	Observatory/site (or "illustrative" for teaching profiles).
outer_scale	float or None	The uniform per-layer outer scale L0 [m], or None if it varies.
wind_direction_measured	bool	Always False here: none of the sources tabulate per-layer wind direction, so every profile's directions are illustrative placeholders.
caveat	str	One-line summary of the representativeness/uncertainty to carry with results.
conditions	mapping or None	Observing conditions the source publishes with the profile, when it does (e.g. r0 [m, 0.5 um], tau0 [s], seeing_class [arcsec], probability, quality); None otherwise. Pass r0=info.conditions["r0"] to reproduce the published case.
origin	str	How the numbers were obtained: "table" (transcribed from a printed table), "dataset" (computed by pyturb from the authors' published data) or "figure" (digitised from a plotted curve, so accurate only to the plot's resolution; see the caveat). Representative and illustrative profiles report "table" for their own definitions.

pyturb.discretize_cn2

discretize_cn2(heights, cn2, n_layers=10, wind='bufton', L0=25.0, method='equivalent', wind_direction=0.0)

Bin a continuous :math:C_n^2(h) profile into equivalent layers.

Each output layer carries its bin's integrated :math:C_n^2\,dh, so the total turbulence is always preserved. The methods differ in how the bin edges are chosen and in the single height (and wind) assigned to each bin:

• method="equivalent" (default) — contiguous log-spaced bin edges, each bin represented by the moment-conserving equivalent layer of Fusco (1999):

.. math:: h_{\mathrm{eq}} = \Big(\frac{\int C_n^2\,h^{5/3}\,dh} {\int C_n^2\,dh}\Big)^{3/5},

the :math:h^{5/3} moment that sets :math:\bar h and hence the isoplanatic angle :math:\theta_0. When wind is given on the input grid the per-layer speed uses the matching :math:v^{5/3} moment, so the coherence time :math:\tau_0 is conserved too. - method="centroid" — log-spaced edges, but each bin uses the simpler :math:C_n^2-weighted mean height (first moment). Preserves the turbulence centroid but not :math:\theta_0/:math:\tau_0 exactly. - method="optimal_grouping" — the bin edges are optimised (not fixed log-spaced) to minimise the total :math:C_n^2-weighted within-group variance of :math:h^{5/3}, so each output layer groups input layers that are close in the tomographic (:math:\theta_0) sense. Chosen by exact dynamic-programming partition; each group then uses the same moment-conserving equivalent height/wind as "equivalent". This is the profile-compression choice for MCAO/atmospheric-tomography reconstruction layers, where where the layers sit matters more than an even altitude spacing (Saxenhuber et al. 2017).

Parameters:

Name	Type	Description	Default
heights	array_like	Sorted altitude grid [m].	required
cn2	array_like	:math:C_n^2 on that grid [m^{-2/3}].	required
n_layers	int	Number of output layers.	10
wind	str, float, or array_like	"bufton" for the Bufton model, a scalar for uniform wind, or a per-layer array of speeds [m/s]. An array matching heights is treated as wind on the input grid and (for method="equivalent") moment-averaged per bin so that :math:\tau_0 is conserved.	'bufton'
L0	float	Outer scale assigned to every layer [m].	25.0
method	('equivalent', 'centroid', 'optimal_grouping')	Edge selection and height/wind assignment rule (see above).	"equivalent"
wind_direction	float or array_like	Wind direction [deg] (the from-direction, see :class:Layer): a scalar for every layer (default 0), an array matching heights (each bin gets the :math:C_n^2-weighted circular mean, so 350 and 10 degrees average to 0, not 180), or one value per output layer.	0.0

pyturb.with_wind

with_wind(layers, speed, direction=None)

Return copies of layers with new wind speeds (and optionally directions).

speed is "bufton" (the Bufton model, :func:bufton_wind, at each layer's altitude), a scalar for every layer, or one value per layer [m/s]. direction [deg, the from-direction, see :class:Layer] is a scalar or one value per layer; None keeps each layer's direction. For profiles whose source publishes no winds (see profile_info(name).conditions).

pyturb.hufnagel_valley

hufnagel_valley(h, wind_rms=21.0, ground=1.7e-14)

Hufnagel-Valley :math:C_n^2(h) model [m^{-2/3}].

Parameters:

Name	Type	Description	Default
h	array_like	Altitude(s) above the telescope [m].	required
wind_rms	float	High-altitude wind pseudo-parameter :math:v [m/s]. 21 with ground=1.7e-14 gives the standard "HV 5/7" profile (r0 ~ 5 cm, isoplanatic angle ~ 7 urad at 500 nm).	21.0
ground	float	Ground-layer coefficient :math:A [m^{-2/3}].	1.7e-14

pyturb.bufton_wind

bufton_wind(h)

Bufton wind-speed model [m/s] versus altitude [m] (peaks near 9.4 km).

Integrated quantities

pyturb.isoplanatic_angle

isoplanatic_angle(layers, r0)

Isoplanatic angle :math:\theta_0 = 0.314\, r_0 / \bar h [rad].

r0 and the layer altitudes must be expressed along the same line of sight (both at zenith, or both scaled for the airmass).

pyturb.coherence_time

coherence_time(layers, r0)

Atmospheric coherence time :math:\tau_0 = 0.314\, r_0 / \bar v [s].

pyturb.greenwood_frequency

greenwood_frequency(layers, r0)

Greenwood frequency :math:f_G = 0.134 / \tau_0 [Hz].

pyturb.mean_turbulence_height

mean_turbulence_height(layers)

Effective turbulence height :math:\bar h = (\sum f_i h_i^{5/3})^{3/5} [m].

pyturb.effective_wind_speed

effective_wind_speed(layers)

Effective wind speed :math:\bar v = (\sum f_i v_i^{5/3})^{3/5} [m/s].

Analysis

pyturb.analysis

Analysis and validation utilities: Zernikes, temporal PSDs, decorrelation.

These turn a pile of phase screens into the diagnostics AO people actually check turbulence against:

• :func:zernike_basis / :func:zernike_decompose — a Noll-ordered Zernike basis on a circular pupil (evaluated by aobasis) and a least-squares projection onto it.
• :func:noll_variance / :func:noll_residual_variance — the Kolmogorov Zernike-mode variances and post-correction residuals of Noll (1976), the textbook thing to validate a decomposition against.
• :func:temporal_psd / :func:fit_power_law — a one-sided temporal power spectrum and a log-log slope fit (frozen-flow gives power-law regimes).
• :func:differential_variance — angular (anisoplanatism) decorrelation.

Everything works on NumPy or CuPy input (device arrays are brought to the host for the reductions). Reference: Noll, R. J. (1976), JOSA 66, 207.

noll_to_zernike

noll_to_zernike(j)

Radial/azimuthal orders (n, m) for Noll single index j (>= 1).

Even j give the cosine terms (m > 0, cos(m theta)) and odd j the sine terms (m < 0, sin(|m| theta)), with theta measured from +x (axis 1) toward +y (axis 0): j = 2 is (1, 1), tip along x. This is the mapping :func:zernike_basis uses; the value is aobasis.noll_to_nm(j).

Parameters:

Name	Type	Description	Default
j	int	Noll index, a Python or NumPy integer >= 1.	required

Returns:

Type	Description
(n, m) : tuple of int

zernike_basis

zernike_basis(n_modes, n_pixels, diameter_pixels=None)

Noll-ordered Zernike basis over a circular pupil.

The polynomials are evaluated by aobasis.zernike_modes_on_mask (the stack's one Zernike implementation) at the pupil pixel centres, which sit at i - (n_pixels - 1) / 2 pixels from the grid centre along each axis (aocore.coordinate_grid).

Axes follow the array convention (y, x): x runs along axis 1 (columns), y along axis 0 (rows), both increasing with the index, and the azimuth theta = atan2(y, x) is measured from +x toward +y. So Noll j = 2 (tip, cos theta) varies along the columns and j = 3 (tilt, sin theta) along the rows.

Parameters:

Name	Type	Description	Default
n_modes	int	Number of modes, starting at Noll j = 1 (piston). At most the number of pixels inside the pupil.	required
n_pixels	int	Grid size; the returned array is (n_modes, n_pixels, n_pixels).	required
diameter_pixels	float	Pupil diameter in pixels (default n_pixels). The pupil is the disc of pixel centres within diameter_pixels / 2 of the grid centre; values are 0 outside it.	None

Returns:

Name	Type	Description
basis	ndarray	(n_modes, n_pixels, n_pixels) float64. Each mode carries the Noll normalisation (unit RMS over the continuous unit disc), and the pupil mask is basis[0] != 0 (piston).

zernike_decompose

zernike_decompose(phase, n_modes, basis=None)

Least-squares Zernike coefficients of phase over the pupil.

Parameters:

Name	Type	Description	Default
phase	ndarray	(n, n) screen or (count, n, n) stack (NumPy or CuPy).	required
n_modes	int	Number of Noll modes to fit.	required
basis	ndarray	A precomputed :func:zernike_basis (reused across many calls to avoid rebuilding it). Must match n_modes and the screen size.	None

Returns:

Name	Type	Description
coeffs	ndarray	(n_modes,) or (count, n_modes) coefficients in the same units as phase.

noll_variance

noll_variance(j, diameter, r0)

Kolmogorov variance of Zernike mode j [rad^2] (Noll 1976).

The per-mode variance is Delta_{j-1} - Delta_j in (D/r0)^{5/3} units. j = 1 (piston) has no finite variance, so j >= 2.

noll_residual_variance

noll_residual_variance(j, diameter, r0)

Residual wavefront variance [rad^2] after correcting the first j Zernike modes (Noll 1976).

temporal_psd

temporal_psd(series, dt, window='hann')

One-sided temporal power spectral density of a time series.

Parameters:

Name	Type	Description	Default
series	array_like	Samples along the last axis (e.g. a pupil pixel or a Zernike coefficient over frames); leading axes are treated as independent series and averaged.	required
dt	float	Sample spacing [s].	required
window	(str, array_like or None)	Taper applied to each (mean-removed) series before the transform. "hann" (default) or any :func:scipy.signal.get_window name, an explicit array of length n, or None for the untapered periodogram. Frozen-flow spectra are steep (f^{-8/3} for one pupil point, steeper for low-order modes), and an untapered periodogram of a finite, non-periodic record leaks power from low to high frequencies: it flattens the measured slope and raises the level. Use None only for records that are exactly periodic.	'hann'

Returns:

Type	Description
freq, psd : ndarray	Positive frequencies [Hz] (excluding DC) and the averaged PSD, scaled by the window's power so that sum(psd) * df approximates the series variance for any window.

fit_power_law

fit_power_law(freq, psd, fmin=None, fmax=None)

Fit psd ~ freq**slope over [fmin, fmax] (log-log least squares).

Returns:

Type	Description
slope, amplitude : float	Such that psd ≈ amplitude * freq**slope in the band.

differential_variance

differential_variance(reference, other)

Variance of other - reference [same units squared].

With reference the on-axis OPD/phase and other an off-axis one, this is the angular (anisoplanatism) error; it grows as (theta/theta0)^{5/3} and reaches ~1 rad^2 at the isoplanatic angle.

Interop

pyturb.interop.HCIPyLayer

A :class:pyturb.Atmosphere presented as an HCIPy atmospheric layer.

Parameters:

Name	Type	Description	Default
atmosphere	Atmosphere	The turbulence to apply. Its pupil sampling must match the HCIPy grid: atmosphere.n points across with spacing atmosphere.pixel_scale (e.g. hcipy.make_pupil_grid(atm.n, atm.diameter)).	required
pupil_grid	CartesianGrid	If given, checked against the atmosphere's sampling at construction.	None

Notes

Field ordering: an HCIPy pupil field reshaped to 2-D is indexed [y, x], the same (y, x) convention as pyturb's OPD arrays, which are flattened onto it row-major; pyturb's x (axis 1) is HCIPy's x and its y (axis 0) HCIPy's y, so direction tuples (thx, thy) carry over unchanged. pyturb's Layer.wind_vector points where the wind comes from (the pattern moves along -wind_vector); HCIPy 0.7's InfiniteAtmosphericLayer velocity behaves the same way, so its velocity=(vx, vy) is wind_direction = degrees(atan2(vy, vx)). HCIPy's FiniteAtmosphericLayer uses a different sign and axis order; see the Interop docs.

The OPD at the current time is evaluated once and reused for every wavelength (WFS and science), converted with the atmosphere's dispersion model. engine="extrude" atmospheres only move forward in time; call :meth:reset to start over.

Examples:

>>> import hcipy, pyturb
>>> atm = pyturb.Atmosphere.from_profile("paranal-median", seeing=0.8,
...                                      n=128, diameter=8.0, engine="extrude")
>>> grid = hcipy.make_pupil_grid(128, 8.0)
>>> layer = pyturb.interop.HCIPyLayer(atm, grid)
>>> layer.t += 1e-3
>>> wf = layer(hcipy.Wavefront(hcipy.make_circular_aperture(8.0)(grid), 1.6e-6))

t property writable

t

Simulation time [s]; setting it moves the turbulence to that time.

evolve_until

evolve_until(t)

Move to time t [s] (frozen flow; never backwards on the extruder).

reset

reset(make_independent_realization=False)

Return to t = 0 with the same turbulence.

make_independent_realization is accepted for HCIPy compatibility; for a new realisation build a new Atmosphere with another seed.

opd

opd()

OPD [m] at the current time, flattened in HCIPy field order.

phase_for

phase_for(wavelength)

Phase [rad] at wavelength [m], flattened in HCIPy field order.

forward

forward(wavefront)

Return a copy of wavefront with the turbulent phase applied.

backward

backward(wavefront)

Return a copy of wavefront with the turbulent phase removed.

Theory

pyturb.theory

Closed-form and numerically integrated reference curves for turbulence.

These are the curves pyturb validates itself against, made public so you can check your own AO simulation against the same numbers, including at finite outer scale L0, where the Kolmogorov tables (Noll 1976) overstate the low-order power. Everything is phase in rad^2 at the wavelength r0 is quoted at, unless a function says otherwise.

• :func:structure_function — von Kármán / Kolmogorov phase structure function D(r).
• :func:zernike_variance — variance of Noll mode j for a pupil of diameter D, any L0.
• :func:image_motion_variance — one-axis tip/tilt image motion [arcsec^2].
• :func:seeing_fwhm — long-exposure FWHM including the outer-scale reduction (Tokovinin 2002).
• :func:pixel_temporal_psd and :func:zernike_temporal_psd — frozen-flow temporal spectra of one pupil point and of a Zernike mode.
• :func:differential_phase_variance — point-wise anisoplanatic phase variance between two directions through a layered profile.
• :func:covariance and :func:space_time_structure_function — the space-time covariance and structure function between pupil points at a time lag, for a layered profile under frozen flow, optionally with pyturb's boiling model (what :class:pyturb.Atmosphere simulates; see :meth:pyturb.Atmosphere.covariance).
• :func:travel_decorrelation_time — an on-sky boiling time constant per layer from its wind speed (Guesalaga et al. 2014), for tau_boil.
• :func:gems_boil_time — the outer-scale tau_boil per layer for the default scale-dependent boiling, calibrated on the same GeMS measurements.
• :func:scintillation_transfer and :func:scintillation_index — the weak-perturbation intensity transfer function of a layer and the scintillation index of a profile (Roddier 1981), the references for :meth:pyturb.Atmosphere.pupil_field.

References

• Noll, R. J. (1976), JOSA 66, 207.
• Conan, J.-M., Rousset, G. & Madec, P.-Y. (1995), JOSA A 12, 1559.
• Tokovinin, A. (2002), PASP 114, 1156.
• Sasiela, R. J. (1994), Electromagnetic Wave Propagation in Turbulence.
• Guesalaga, A., Neichel, B., Cortés, A., Béchet, C. & Guzmán, D. (2014), "Using the Cn2 and wind profiler method with wide-field laser-guide-stars adaptive optics to quantify the frozen-flow decay", MNRAS 440, 1925-1933, doi:10.1093/mnras/stu366.
• Roddier, F. (1981), "The effects of atmospheric turbulence in optical astronomy", Progress in Optics 19, 281-376.

structure_function

structure_function(r, r0, L0=np.inf)

Phase structure function D(r) [rad^2] at separation r [m].

6.88 (r/r0)^{5/3} for Kolmogorov (L0=inf); 2 [C(0) - C(r)] from :func:pyturb.phase_covariance for finite L0.

zernike_variance

zernike_variance(j, diameter, r0, L0=np.inf)

Variance [rad^2] of Noll Zernike mode j (>= 2) over a circular pupil.

integral PSD(f) |Q_j(f)|^2 d^2f with the von Kármán PSD, evaluated numerically. For L0=inf it reproduces Noll's Kolmogorov values (:func:pyturb.analysis.noll_variance); a finite outer scale lowers the low orders most (tip/tilt by ~40% at L0/D = 3).

image_motion_variance

image_motion_variance(diameter, r0, wavelength, L0=np.inf)

One-axis image-motion (tilt) variance [arcsec^2] for a pupil of diameter D.

r0 is at wavelength [m]. For Kolmogorov turbulence this is the textbook 0.182 (D/r0)^{5/3} (lambda/D)^2 (Zernike tilt); a finite L0 reduces it substantially for large telescopes.

seeing_fwhm

seeing_fwhm(r0, wavelength=5e-07, L0=np.inf)

Long-exposure seeing FWHM [arcsec], with the outer-scale correction.

0.98 lambda/r0 (as :func:pyturb.seeing_from_r0) times sqrt(1 - 2.183 (r0/L0)^{0.356}) (Tokovinin 2002, accurate for L0/r0 > 20). At L0 = 25 m and r0 = 0.1 m the delivered FWHM is ~17% below the Kolmogorov value; :attr:pyturb.Atmosphere.seeing is the Kolmogorov value.

pixel_temporal_psd

pixel_temporal_psd(freq, r0, wind_speed)

One-sided temporal PSD [rad^2/Hz] of one pupil point under frozen flow.

Kolmogorov, single layer: 0.0774 r0^{-5/3} V^{5/3} f^{-8/3}.

zernike_temporal_psd

zernike_temporal_psd(freq, j, diameter, r0, wind_speed, wind_direction=0.0, L0=np.inf)

One-sided temporal PSD [rad^2/Hz] of Noll mode j for one frozen layer.

PSD_j(nu) = 2/V integral PSD(f) |Q_j(f)|^2 df_perp with f_parallel = nu / V (Conan, Rousset & Madec 1995), integrated numerically. wind_direction [deg] is in the same frame as :attr:pyturb.Layer.wind_direction and :func:pyturb.analysis.zernike_basis (from +x = axis 1 toward +y = axis 0), so wind_direction=0 is wind along x, where Noll j = 2 (tip) is the along-wind tilt. The frame matters for non-symmetric modes: tilt along the wind has a different spectrum from tilt across it. Integrates (over nu) to :func:zernike_variance. For a profile, sum the layers' spectra, each with its own r0_i and wind.

differential_phase_variance

differential_phase_variance(theta, layers, r0, L0=None)

Point-wise phase variance [rad^2] between directions theta [arcsec] apart.

sum_i D_i(h_i * theta), each layer's structure function at its own footprint separation, with r0_i^{-5/3} = f_i r0^{-5/3}. r0 and the layer altitudes must be along the same line of sight. For Kolmogorov turbulence (L0=numpy.inf) this is (theta/theta0)^{5/3}, ~1 rad^2 at the isoplanatic angle; a finite outer scale lowers it. L0=None (default) uses each layer's own L0. No aperture averaging: a pupil-averaged, piston-removed differential variance is smaller.

covariance

covariance(points_a, points_b, lag=0.0, *, layers, r0, L0=None, tau_boil=None, boil_model='ar1', boil_scaling='kolmogorov', units='phase', wavelength=5e-07)

Space-time covariance <phi(a, t) phi(b, t + lag)> of a layered atmosphere.

The statistics :class:pyturb.Atmosphere simulates, for estimator design (MMSE, predictive, Kalman): von Kármán layers under frozen flow, each following pyturb's convention phi(r, t) = phi_0(r + wind_vector * t), optionally boiling. Layer i contributes f_i B_i(|b - a + wind_vector_i * lag|) (f_i its normalised cn2_fraction) with

B_i(rho) = 2 pi int f PSD_i(f) R(|lag| / tau_i(f)) J0(2 pi f rho) df,

the von Kármán covariance (:func:pyturb.phase_covariance) when the layer does not boil or lag = 0. R is the per-mode retention of boil_model: exp(-s) for "ar1", (1 + 2s) exp(-2s) for "ar2". With boil_scaling="kolmogorov" (the spectral engine) a mode at spatial frequency f has tau(f) = tau_boil * (f L0)^(-2/3) for f >= 1/L0 and tau_boil below, so tau_boil is the lifetime of the outer-scale structure and fine structure decorrelates faster; the Hankel transform is then evaluated numerically (relative accuracy ~1e-7). With "uniform" (the extruder, and the spectral engine's boil_scaling="uniform") every scale has tau_boil, and the covariance is the frozen-flow one times R(|lag| / tau_boil).

Parameters:

Name	Type	Description	Default
points_a	array_like	Pupil-plane points [m], one (x, y) or an (N, 2) array, with x along axis 1 (columns) and y along axis 0 (rows), as in :meth:pyturb.Atmosphere.opd_at.	required
points_b	array_like	Pupil-plane points [m], one (x, y) or an (N, 2) array, with x along axis 1 (columns) and y along axis 0 (rows), as in :meth:pyturb.Atmosphere.opd_at.	required
lag	float or 1-D array_like	Time [s] from points_a to points_b. A 1-D array adds a leading lag axis to the result.	0.0
layers	sequence of Layer	The profile: cn2_fraction (normalised to sum to 1), wind vector and L0 of each layer. Altitudes are not used (one source at infinity).	required
r0	float	Total Fried parameter [m] along the line of sight, at wavelength.	required
L0	float	Outer scale [m] for every layer, overriding each layer's own. Must be finite (the Kolmogorov variance is infinite); see :func:space_time_structure_function for L0=inf.	None
tau_boil	float or sequence	Boiling time constant [s], scalar or one per layer (numpy.inf for a frozen layer), as :class:pyturb.Atmosphere's tau_boil. None (default) is pure frozen flow.	None
boil_model	('ar1', 'ar2')	Temporal boiling model, as :class:pyturb.Atmosphere's boil_model.	"ar1"
boil_scaling	('kolmogorov', 'uniform')	Scale dependence of boiling, as :class:pyturb.Atmosphere's boil_scaling: "kolmogorov" (default, the spectral engine's default) or "uniform" (the extruder, or the spectral engine with boil_scaling="uniform").	"kolmogorov"
units	('phase', 'opd')	"phase" (default): rad^2 at wavelength. "opd": m^2.	"phase"
wavelength	float	Wavelength [m] r0 is given at. Default 500 nm.	5e-07

Returns:

Type	Description
ndarray	(Na, Nb), or (T, Na, Nb) for T lags.

Notes

The result is a stationary, infinite-screen statistic. A simulated screen differs in what its finite grid cannot hold: power above the grid's Nyquist frequency, and on the periodic spectral engine the largest scales (see oversample). Differences such as the structure function are far less sensitive to the large scales than the covariance itself.

Examples:

>>> import numpy as np, pyturb
>>> layers = [pyturb.Layer(0.0, 0.7, 10.0, 0.0), pyturb.Layer(8e3, 0.3, 25.0, 90.0)]
>>> pts = np.array([[0.0, 0.0], [0.5, 0.0]])
>>> c = pyturb.theory.covariance(pts, pts, lag=2e-3, layers=layers, r0=0.15,
...                              tau_boil=1.0, units="opd")
>>> c.shape
(2, 2)

space_time_structure_function

space_time_structure_function(points_a, points_b, lag=0.0, *, layers, r0, L0=None, tau_boil=None, boil_model='ar1', boil_scaling='kolmogorov', units='phase', wavelength=5e-07)

Space-time structure function <[phi(a, t) - phi(b, t + lag)]^2>.

2 [C(0) - B] with B from :func:covariance (same arguments and model), computed without that subtraction, so it stays accurate for close points and short lags. Unlike :func:covariance it is finite for Kolmogorov layers (L0=inf), as long as they do not boil: with boiling, modes at every scale decorrelate while carrying unbounded power. At lag=0 it is :func:structure_function summed over the layers; at a = b it is the temporal structure function of one pupil point (the mean-square change over lag).

travel_decorrelation_time

travel_decorrelation_time(speed, per_metre=0.157, per_second=0.365)

Boiling time constant [s] of a layer from the distance its wind carries it.

tau = 1 / (per_metre * |speed| + per_second). Guesalaga et al. (2014) followed the turbulence of individual layers with the GeMS wide-field laser-guide-star AO system and found each layer's correlation decaying roughly linearly with lag, rho(lag) ~ 1 - r lag, with a rate that grows with wind speed, r = 0.157 v + 0.365 s^-1 (v in m/s): a layer loses correlation mainly with the distance it has travelled, plus a slower part that does not depend on the wind. The defaults are that fit; 1 / r is the time constant whose AR(1) retention exp(-lag/tau) has the same initial slope.

The fit covers 74 GeMS layer sequences with winds of 0-24 m/s, as a whole-layer, roughly linear decay; it says nothing about how the decay depends on spatial scale. Use it with boil_scaling="uniform", which gives every scale this time constant. With the default boil_scaling="kolmogorov", tau_boil is the lifetime of the outer scale only and the scales an AO system senses would boil one to two orders of magnitude faster than measured. boil_model="ar2" keeps the same integral time scale but starts flat (quadratic) instead of linear at short lags. Treat values outside the fitted speed range as extrapolation.

Parameters:

Name	Type	Description	Default
speed	float or array_like	Layer wind speed [m/s]; the sign is ignored. An array gives one time constant per layer, e.g. [layer.wind_speed for layer in layers].	required
per_metre	float	Decorrelation per metre of travel [m^-1]. Default 0.157.	0.157
per_second	float	Wind-independent decorrelation rate [s^-1]. Default 0.365.	0.365

Returns:

Type	Description
float or ndarray	tau [s], a float for a scalar speed, else an array of its shape. inf (frozen flow) where the rate is zero.

References

Guesalaga, A., Neichel, B., Cortés, A., Béchet, C. & Guzmán, D. (2014), "Using the Cn2 and wind profiler method with wide-field laser-guide-stars adaptive optics to quantify the frozen-flow decay", MNRAS 440, 1925-1933, doi:10.1093/mnras/stu366: per-layer frozen-flow decay measured on GeMS.

Examples:

>>> import pyturb
>>> layers = pyturb.get_profile("two-layer")
>>> tau = pyturb.theory.travel_decorrelation_time(
...     [layer.wind_speed for layer in layers])
>>> atm = pyturb.Atmosphere(layers, r0=0.15, n=64, tau_boil=tau,
...                         boil_scaling="uniform")

gems_boil_time

gems_boil_time(speed)

Outer-scale boiling lifetime [s] for boil_scaling="kolmogorov", from GeMS.

The tau_boil per layer for the spectral engine's default, scale-dependent boiling (each mode at spatial frequency f has tau(f) = tau_boil * max(f L0, 1)^(-2/3)), derived from the decorrelation Guesalaga et al. (2014) measured on GeMS. They followed each layer in the space-time cross-correlation of tip/tilt-removed laser-guide-star slopes (8 m pupil, 16x16 subapertures of 0.5 m), deconvolved by the zero-lag autocorrelation, and took the slope of the energy in a 3x3-subaperture window around the correlation peak beyond 0.1 s, finding a rate 0.157 v + 0.365 s^-1 for wind speed v (see :func:travel_decorrelation_time). The lifetimes here, 2.5, 1.4 and 0.74 s at 5, 10 and 20 m/s (2.50, 1.44, 0.74 s before rounding), are those for which the analytic windowed energy of the von Kármán spectrum (L0 = 25 m, AR(1) modes, not deconvolved) has that rate as the slope of its logarithm over 0.1-0.5 s. In between, 1 / tau is linear in speed; outside 5-20 m/s the end values are held.

The measured rate is a linear slope of the deconvolved statistic, not a logarithmic one. Under that definition, AR(1) modes with these lifetimes decay at about 0.9, 0.6 and 0.5 times the published rate at 5, 10 and 20 m/s, and no AR(1) law reaches the published rate above about 6 m/s: over 0.1-0.5 s the linear slope of any mixture of exponential decays is at most about 1.3 s^-1. The GeMS decay is linear in time and includes the dispersion of the wind speed within a layer, which AR(1) boiling of a single layer does not represent. Treat these lifetimes as an approximate, AR(1) description of the GeMS measurements, faithful in order of magnitude.

Use it with the default boil_scaling="kolmogorov" and boil_model="ar1"; with "uniform" scaling use :func:travel_decorrelation_time instead, which describes the same measurements as one rate for every scale. The calibration assumes L0 = 25 m: the outer-scale lifetime that gives the same decay of the sensed scales grows with L0 (recalibrated at 10 m/s, 0.63, 1.49 and 2.3 times these values for L0 = 10, 50 and 100 m), so treat other outer scales as an extrapolation. The calibration covers lags of 0.1-0.5 s; its extrapolation to millisecond lags is not constrained by the GeMS data.

Parameters:

Name	Type	Description	Default
speed	float or array_like	Layer wind speed [m/s]; the sign is ignored. An array gives one value per layer, e.g. [layer.wind_speed for layer in layers].	required

Returns:

Type	Description
float or ndarray	tau_boil [s], a float for a scalar speed, else an array of its shape.

References

Guesalaga, A., Neichel, B., Cortés, A., Béchet, C. & Guzmán, D. (2014), MNRAS 440, 1925-1933, doi:10.1093/mnras/stu366.

Examples:

>>> import pyturb
>>> layers = pyturb.get_profile("two-layer")
>>> tau = pyturb.theory.gems_boil_time([layer.wind_speed for layer in layers])
>>> atm = pyturb.Atmosphere(layers, r0=0.15, n=64, tau_boil=tau)

scintillation_transfer

scintillation_transfer(f, altitude, wavelength)

Weak-perturbation transfer [1/m] from a layer's OPD to relative intensity.

For a thin layer at range h above the pupil with a small OPD w, a plane wave arrives with relative intensity 1 + delta whose spectrum is delta(f) = T(f) w(f), with

T(f) = 2 k sin(pi lambda h f^2), k = 2 pi / lambda,

the Fresnel filter of first-order (Rytov/Born) perturbation theory, in the sign convention of :meth:pyturb.Atmosphere.pupil_field (phase exp(+i k w)): a positive OPD bump (a converging lens) brightens the pupil. Below the Fresnel frequency 1 / sqrt(lambda h) it tends to the achromatic transport-of-intensity limit 4 pi^2 h f^2, i.e. delta = -h laplacian(w). The layers of a profile add their contributions; for a band, average T over the wavelengths with the weights of the intensity. Multiplying |T|^2 by the OPD power spectrum and integrating gives the scintillation index (:func:scintillation_index); multiplying T by other transfer functions gives cross-spectra, e.g. for the covariance between intensity and slopes.

Parameters:

Name	Type	Description	Default
f	array_like	Spatial frequency [cycles/m] (its magnitude).	required
altitude	float or array_like	Range of the layer above the pupil [m] along the line of sight.	required
wavelength	float or array_like	Wavelength [m].	required

Returns:

Type	Description
ndarray	T [relative intensity per metre of OPD], broadcast over the inputs.

scintillation_index

scintillation_index(layers, r0, wavelength, *, r0_wavelength=5e-07, zenith_angle=0.0)

Scintillation index sigma_I^2 of a star through a layered profile (Roddier).

The variance of the relative intensity at a point in the pupil, in the weak-perturbation (Rytov) approximation for Kolmogorov turbulence:

sigma_I^2 = 19.2 lambda^(-7/6) sec(z)^(11/6) integral Cn^2(h) h^(5/6) dh

(four times the log-amplitude variance), with each layer's Cn^2 dh = f_i r0^(-5/3) / (0.423 k_0^2) from the profile's normalised cn2_fraction f_i and the total r0 at r0_wavelength (k_0 = 2 pi / r0_wavelength). It is the integral of :func:scintillation_transfer squared times the Kolmogorov OPD spectrum, so it scales as lambda^(-7/6) and is dominated by the high layers.

It describes a point source seen through a point aperture: a telescope larger than the Fresnel scale sqrt(lambda h) averages the fluctuations down. The weak-perturbation result holds while sigma_I^2 is well below 1; stronger scintillation saturates. The outer scale hardly matters (the Fresnel filter suppresses scales larger than sqrt(lambda h)). A simulation reproduces it only up to its grid: the power above the Nyquist frequency of the pupil sampling is missing (several percent unless the pixel is well below the Fresnel scale).

Parameters:

Name	Type	Description	Default
layers	sequence of Layer	The profile: altitudes [m] at zenith and cn2_fraction weights.	required
r0	float	Total Fried parameter [m] at zenith, at r0_wavelength.	required
wavelength	float or array_like	Observing wavelength(s) [m].	required
r0_wavelength	float	Wavelength [m] r0 is quoted at. Default 500 nm.	5e-07
zenith_angle	float	Zenith angle [deg]: ranges and the path through each layer grow as sec(z). Default 0.	0.0

Returns:

Type	Description
float or ndarray	sigma_I^2 (dimensionless), a float for a scalar wavelength.

References

Roddier, F. (1981), Progress in Optics 19, 281-376.

Examples:

>>> import pyturb
>>> layers = pyturb.get_profile("paranal-median")
>>> s2 = pyturb.theory.scintillation_index(layers, r0=0.15, wavelength=0.5e-6)

I/O

pyturb.save

save(path, data, **metadata)

Write data (a phase/OPD array) to path with metadata.

Parameters:

Name	Type	Description	Default
path	str	Destination. .fits/.fit writes FITS (needs astropy); any other extension writes a NumPy .npz.	required
data	ndarray	Screen or stack of screens (NumPy or CuPy; device arrays are copied to the host). Saved as float32/float64 as-is.	required
**metadata	Any	Arbitrary scalars/strings to record alongside the data, e.g. pixel_scale=0.0156, r0=0.15, L0=25.0, wavelength=500e-9, units="metres", seed=1. None values are dropped. The pyturb version is always added as pyturb_version.	{}

Examples:

>>> import pyturb
>>> atm = pyturb.Atmosphere.from_profile("two-layer", seeing=0.8)
>>> pyturb.save("opd.fits", atm.opd(), **atm.metadata)

pyturb.load

load(path)

Load an array and its metadata written by :func:save.

Returns:

Type	Description
(ndarray, dict)	The data array (on the host) and its metadata dictionary.

Utilities

pyturb.r0_from_seeing

r0_from_seeing(seeing, wavelength=5e-07)

Fried parameter (m) from seeing FWHM (arcsec) at wavelength (m).

Uses the Kolmogorov relation FWHM = 0.98 lambda / r0.

pyturb.seeing_from_r0

seeing_from_r0(r0, wavelength=5e-07)

Seeing FWHM (arcsec) from the Fried parameter (m) at wavelength (m).

pyturb.r0_at_wavelength

r0_at_wavelength(r0, wavelength_in, wavelength_out)

Rescale the Fried parameter between wavelengths (r0 ~ lambda^(6/5)).

Example: convert an r0 quoted at 500 nm to the K band, r0_at_wavelength(0.15, 500e-9, 2.2e-6).

pyturb.opd_to_phase

opd_to_phase(opd, wavelength)

Convert optical path difference [m] to phase [rad] at wavelength [m].

phase = 2 pi * opd / wavelength. OPD is achromatic, so the same OPD gives different phase at different wavelengths. Delegates to aocore.opd_to_phase; NumPy and CuPy arrays keep their device and dtype.

pyturb.phase_to_opd

phase_to_opd(phase, wavelength)

Convert phase [rad] at wavelength [m] to optical path difference [m].

opd = phase * wavelength / (2 pi). Delegates to aocore.phase_to_opd; NumPy and CuPy arrays keep their device and dtype.

pyturb.air_refractivity

air_refractivity(wavelength)

Refractivity n - 1 of standard dry air at wavelength [m].

Edlén (1966) dispersion formula for dry air at 15 °C, 101.325 kPa, 0.03 % CO2::

(n - 1) x 1e8 = 8342.13 + 2406030 / (130 - s^2) + 15997 / (38.9 - s^2)

with s = 1 / lambda in micron^-1. Turbulence OPD scales with n - 1, so the ratio of this quantity between two wavelengths is the (weak, ~1-2 % across the visible-to-NIR) chromatic correction to an otherwise achromatic path length. Temperature/pressure scale n - 1 overall and so cancel in that ratio; only the dispersion shape matters here.

Validity: the Edlén (1966) formula is established over ~0.23-1.7 µm. Its dispersion ratio (the only thing used here) agrees with the independent Ciddor (1996) standard-air formula to < 0.05 % over 0.35-1.7 µm and stays within ~0.1 % out to ~2.5 µm, so the ratio is reliable somewhat past the nominal range; the limiting factor for thermal-IR / interferometric work is not the dry-air dispersion shape but the wet/dry split and the weather-independent wet_fraction (see :func:water_vapour_refractivity and the dispersion argument of :class:pyturb.Atmosphere).

Parameters:

Name	Type	Description	Default
wavelength	float or array_like	Wavelength in metres.	required

Returns:

Type	Description
float or ndarray	n - 1 (dimensionless), same shape as wavelength.

pyturb.water_vapour_refractivity

water_vapour_refractivity(wavelength)

Dispersion of the water-vapour refractivity at wavelength [m].

The wavelength-dependent part of the pure-water-vapour refractivity from Ciddor (1996), Refractive index of air: new equations for the visible and near infrared, Applied Optics 35, 1566::

(n_wv - 1) x 1e8 = 1.022 (295.235 + 2.6422 s^2
                          - 0.032380 s^4 + 0.004028 s^6)

with s = 1 / lambda in micron^-1. The overall scale depends on the water-vapour partial pressure/temperature (a density factor) and cancels in any ratio between wavelengths, so this returns the dispersion shape only — exactly what the wet/dry chromatic OPD split needs.

Water vapour disperses differently from dry air (:func:air_refractivity): across the visible-to-NIR the two shapes track loosely, but into the mid-IR the dry term flattens while the wet term keeps more structure. That divergence is the "wet–dry" problem in interferometry, where the water-vapour contribution to the turbulent OPD dominates.

Validity: the Ciddor (1996) water-vapour dispersion is established over ~0.35-1.7 µm (extrapolated beyond). It gives only the dispersion shape; the density weighting (partial pressure, temperature) and the split between wet and dry turbulent refractivity are not modelled from weather inputs — the caller supplies a single scalar wet_fraction (see the dispersion/wet_fraction arguments of :class:pyturb.Atmosphere), which is the real accuracy limit for thermal-IR / interferometric OPD, not the dispersion formula itself.

Parameters:

Name	Type	Description	Default
wavelength	float or array_like	Wavelength in metres.	required

Returns:

Type	Description
float or ndarray	A quantity proportional to the water-vapour n - 1 (dimensionless), same shape as wavelength.

pyturb.structure_function

structure_function(phase, pixel_scale=1.0, max_separation=None)

Phase structure function averaged over the two array axes.

D(r) = <[phase(x) - phase(x + r)]^2>, estimated from pixel pairs separated along axis 0 and along axis 1 and averaged (not a full azimuthal average: diagonal separations are not used). For Kolmogorov turbulence the expectation is D(r) = 6.88 (r / r0)^(5/3).

Parameters:

Name	Type	Description	Default
phase	ndarray	A single screen (n, n) or a stack (count, n, n); NumPy or CuPy (device arrays are copied to the host).	required
pixel_scale	float	Pixel size in metres; separations are returned in the same unit.	1.0
max_separation	int	Largest pixel separation to evaluate (default n // 4).	None

Returns:

Type	Description
r, D : ndarray	Separations and the structure function estimate at each separation.

pyturb.set_fft_workers

set_fft_workers(workers)

Set the thread count for CPU (SciPy) FFTs; affects all pyturb objects.

None (default) is single-threaded; -1 uses every core; any positive integer pins that many. Other negative values follow scipy.fft wraparound semantics (os.cpu_count() + 1 + workers, e.g. -2 = all cores but one). 0 is rejected. No effect on the GPU path (CuPy). Returns the previous value.

import pyturb # doctest: +SKIP pyturb.set_fft_workers(-1) # use all cores for CPU FFTs

pyturb.get_fft_workers

get_fft_workers()

Return the current CPU FFT thread setting (see :func:set_fft_workers).

pyturb.set_threads

set_threads(threads)

Cap the CPU threads pyturb uses; returns the previous cap.

threads (an integer >= 1) limits, for this process:

• pyturb's parallel Numba kernels (numba.set_num_threads, at most the NUMBA_NUM_THREADS pool Numba started with), in every Python thread;
• the SciPy FFT workers, if threaded FFTs were enabled with :func:set_fft_workers (it never turns them on).

None removes the cap (Numba back to its full pool). Results do not depend on the setting: every parallel kernel computes each output element in the same order whatever the thread count.

It leaves the BLAS thread pool alone, because OpenBLAS's rounding depends on its thread count (a seeded screen changes in the last bit). Limit BLAS with OPENBLAS_NUM_THREADS/OMP_NUM_THREADS before starting Python, or threadpoolctl.threadpool_limits, and expect last-bit differences from runs with a different BLAS thread count.

Use it in multiprocessing workers, where every process would otherwise start a thread per core, e.g. Pool(initializer=pyturb.set_threads, initargs=(1,)). Without calling it, pyturb honours OMP_NUM_THREADS for its Numba kernels when NUMBA_NUM_THREADS is not set.

import pyturb # doctest: +SKIP pyturb.set_threads(1) # one thread per worker process

pyturb.get_threads

get_threads()

Return the current thread cap (see :func:set_threads); None if unset.

Backend

pyturb.to_numpy

to_numpy(array)

Copy an array to host memory as a numpy.ndarray.

A no-op (beyond asarray) for arrays that are already on the CPU.

pyturb.get_array_module

get_array_module(device)

Return the array module (numpy or cupy) for a device name.

Parameters:

Name	Type	Description	Default
device	str	"cpu" (alias "numpy") or "gpu" (aliases "cuda", "cupy"), optionally with a GPU number: "gpu:1". Without a number the current CuPy device is used (normally GPU 0).	required

pyturb.benchmark

pyturb.benchmark() — print your machine's throughput in a few seconds.

benchmark

benchmark(n=512, profile='paranal-median', device='cpu', seconds=1.0, seeing=0.8, engine='spectral')

Measure and print closed-loop and Monte-Carlo throughput.

Parameters:

Name	Type	Description	Default
n	int	Pupil sampling.	512
profile	str	Named turbulence profile.	'paranal-median'
device	str	"cpu" or "gpu".	'cpu'
seconds	float	Wall-clock budget per measurement.	1.0
seeing	float	Passed to :class:Atmosphere.	0.8
engine	float	Passed to :class:Atmosphere.	0.8

Returns:

Type	Description
dict	{"frames_per_s": ..., "screens_per_s": ...}.