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. When set, the layer relaxes via an AR(1) process toward fresh, statistically identical turbulence with retention exp(-dt/tau) while preserving its spatial statistics. 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, so 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. 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.	None
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.	None

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 (finer structure faster, scale-resolved), 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).
• 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.

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, not a full replayable checkpoint: the per-layer arrays (altitudes, Cn2 fractions, winds) and any evolved/boiled stochastic state are not serialised here, so it cannot by itself reconstruct a specific evolved frame. 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, **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)

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.

r0_at

r0_at(wavelength)

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

opd

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

OPD at time t [s] under frozen flow.

Parameters:

Name	Type	Description	Default
t	float	Time since the start of the simulation [s]. Taylor frozen flow: each layer's phase evolves as phi(x, t) = phi_0(x + wind_vector * t) — 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 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) 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

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.

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.

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. Screens are reproducible for a fixed seed, backend and dtype.	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.

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.	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)         # blown 1.3 m / -0.4 m

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 axis 0 (rows), sy along axis 1 (columns). Values are phase in radians at the layer's reference wavelength; the array lives on the layer's device.

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.	required
wind_speed	float	Wind speed [m/s].	10.0
wind_direction	float	Wind direction [deg], measured from axis 0 toward axis 1. Arbitrary (need not be axis-aligned).	0.0
L0	float	Outer scale [m] for this layer. Default 25 m.	25.0

wind_vector property

wind_vector

(vx, vy) wind components [m/s] along axes 0 and 1.

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". See :func:list_profiles.

"mauna-kea", "keck" and "las-campanas" are traceable to a specific published table (see each profile-building function's source comment for the citation); "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.discretize_cn2

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

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"

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 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).

zernike_basis

zernike_basis(n_modes, n_pixels, diameter_pixels=None)

Noll-ordered Zernike basis over a circular pupil.

Parameters:

Name	Type	Description	Default
n_modes	int	Number of modes, starting at Noll j = 1 (piston).	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 inscribed disc; values are 0 outside it.	None

Returns:

Name	Type	Description
basis	ndarray	(n_modes, n_pixels, n_pixels). Each mode is orthonormalised so that its variance over the pupil is 1 (Noll normalisation), 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)

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

Returns:

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

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.

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.

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).

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)

Azimuthally averaged (along both axes) phase structure function.

D(r) = <[phase(x) - phase(x + r)]^2>, estimated from pixel pairs separated along each image axis and averaged. 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.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": ...}.