Simulating Atmospheric Turbulence

We’ll model how Earth’s atmosphere distorts starlight using HCIPy’s phase screen framework, from single-layer models to realistic multi-layer atmospheric profiles.

Atmospheric turbulence causes spatial and temporal variations in the refractive index along the line of sight, which translates to optical path length differences — phase aberrations — that degrade image quality. In this tutorial we’ll build up from a single turbulent layer to a full multi-layer profile, exploring how parameters like the Fried parameter \(r_0\) and wind speed affect the Point Spread Function (PSF).

import hcipy
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm

We start by defining the telescope and observing parameters. The aperture diameter sets the diffraction limit, while the wavelength determines the spatial scale of the focal plane. We’ll use a \(256 \times 256\) pupil grid with generous padding to avoid aliasing, and a focal grid that oversamples the PSF by a factor of 4 (\(q = 4\)) to resolve fine structure.

# Telescope parameters
pupil_diameter = 8.2  # meters
wavelength = 1.5e-6  # meters

# Computational grids
pupil_grid = hcipy.make_pupil_grid(256, pupil_diameter * 1.2)
spatial_resolution = wavelength / pupil_diameter
focal_grid = hcipy.make_focal_grid(q=4, num_airy=32, spatial_resolution=wavelength / pupil_diameter)

Now that we have our grids defined, we can set up the Fraunhofer propagator and the telescope aperture. The propagator will transform the electric field from the pupil plane to the focal plane.

# Propagator and aperture
prop = hcipy.FraunhoferPropagator(pupil_grid, focal_grid)
aperture = hcipy.make_vlt_aperture()(pupil_grid)

print(f"Telescope diameter: {pupil_diameter} m")
print(f"Observing wavelength: {wavelength*1e9:.0f} nm")
print(f"Pupil grid: {pupil_grid.dims} points")
print(f"Focal grid: {focal_grid.dims} points")

Telescope diameter: 8.2 m
Observing wavelength: 1500 nm
Pupil grid: (256, 256) points
Focal grid: (256, 256) points

The core of an atmospheric simulation is the phase screen — a 2D map of the optical path difference introduced by refractive index fluctuations. We’ll use HCIPy’s InfiniteAtmosphericLayer, which generates scrollable turbulence patterns using the method of Assemat et al. (2006). This type of layer synthesizes an infinitely large phase screen on-the-fly by summing Fourier modes with Kolmogorov-von Kármán statistics, so wind can blow the screen across the aperture without repeating.

The key parameters are the Fried parameter \(r_0\) (which sets the turbulence strength), the outer scale \(L_0\) (the largest eddy size), and the wind velocity (which controls the frozen-flow evolution rate).

# Atmospheric parameters
fried_parameter = 0.2  # r0 = 20 cm
outer_scale = 20  # L0 (meters)
velocity = 10  # wind speed (m/s)

# Convert Fried parameter to Cn squared
Cn_squared = hcipy.Cn_squared_from_fried_parameter(fried_parameter, wavelength=500e-9)

print("Atmospheric Parameters:")
print(f"  Fried parameter (r0): {fried_parameter} m at 500 nm")
print(f"  Outer scale (L0): {outer_scale} m")
print(f"  Wind velocity: {velocity} m/s")
print(f"  Cn squared: {Cn_squared:.2e} m^(-2/3)")

Atmospheric Parameters:
  Fried parameter (r0): 0.2 m at 500 nm
  Outer scale (L0): 20 m
  Wind velocity: 10 m/s
  Cn squared: 2.19e-13 m^(-2/3)

With these parameters defined, we can create the InfiniteAtmosphericLayer. The layer will use these atmospheric parameters to generate phase screens consistent with Kolmogorov-von Kármán turbulence statistics.

# Atmospheric layer
layer = hcipy.InfiniteAtmosphericLayer(
    pupil_grid,
    Cn_squared=Cn_squared,
    L0=outer_scale,
    velocity=velocity,
    seed=0
)

print("\nAtmospheric phase screen created")

Atmospheric phase screen created

Let’s take an instantaneous snapshot of the phase screen and compute wavefront error statistics across the aperture.

# Phase screen at observing wavelength
phase = layer.phase_for(wavelength)

# Visualize
plt.figure(figsize=(10, 8))
im = hcipy.imshow_field(phase, cmap='RdBu_r', mask=aperture, vmin=-10, vmax=10)
plt.colorbar(im, label='Phase (radians)')
plt.title(f'Atmospheric Phase Screen\nr0 = {fried_parameter}m at {wavelength*1e9:.0f}nm')
plt.xlabel('x (m)')
plt.ylabel('y (m)')
plt.axis('equal')
plt.tight_layout()
plt.show()

# Wavefront error statistics
phase_inside_aperture = phase[aperture > 0]
rms_wfe = np.sqrt(np.mean(phase_inside_aperture**2))
peak_to_valley = np.max(phase_inside_aperture) - np.min(phase_inside_aperture)

print(f"Wavefront Error Statistics (inside aperture):")
print(f"  RMS: {rms_wfe:.2f} rad ({rms_wfe * wavelength * 1e9:.1f} nm)")
print(f"  Peak-to-Valley: {peak_to_valley:.2f} rad ({peak_to_valley * wavelength * 1e9:.1f} nm)")

Wavefront Error Statistics (inside aperture):
  RMS: 2.83 rad (4247.1 nm)
  Peak-to-Valley: 12.95 rad (19418.7 nm)

The RMS wavefront error is on the order of several radians — significantly larger than the \(\lambda/14\) (≈0.45 rad) Marechal criterion for diffraction-limited imaging. These phase errors scatter light from the core of the PSF into a broad halo. Let’s quantify this by comparing the diffraction-limited PSF to the turbulence-degraded version.

# Plane wave entering the telescope
wf = hcipy.Wavefront(aperture, wavelength)

# Apply atmosphere and propagate
wf_atmos = layer(wf)
img_turbulent = prop(wf_atmos)

# Diffraction-limited comparison
img_diffraction_limited = prop(wf)

# Visualize
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Diffraction-limited PSF
im0 = hcipy.imshow_field(
    np.log10(img_diffraction_limited.intensity / img_diffraction_limited.intensity.max()),
    vmin=-5, vmax=0, cmap='inferno', grid_units=spatial_resolution, ax=axes[0]
)
axes[0].set_title('Diffraction-Limited (No Atmosphere)', fontsize=12)
axes[0].set_xlabel('x (lambda/D)')
axes[0].set_ylabel('y (lambda/D)')
plt.colorbar(im0, ax=axes[0], label='log10(Intensity)')

# Turbulent PSF
im1 = hcipy.imshow_field(
    np.log10(img_turbulent.intensity / img_turbulent.intensity.max()),
    vmin=-5, vmax=0, cmap='inferno', grid_units=spatial_resolution, ax=axes[1]
)
axes[1].set_title(f'With Atmospheric Turbulence (r0 = {fried_parameter}m)', fontsize=12)
axes[1].set_xlabel('x (lambda/D)')
axes[1].set_ylabel('y (lambda/D)')
plt.colorbar(im1, ax=axes[1], label='log10(Intensity)')

plt.tight_layout()
plt.show()

# Strehl ratio
strehl_dl = np.max(img_diffraction_limited.intensity) / np.sum(img_diffraction_limited.intensity)
strehl_turb = np.max(img_turbulent.intensity) / np.sum(img_turbulent.intensity)
print(f"\nImage Quality Comparison:")
print(f"  Diffraction-limited peak intensity: {strehl_dl:.4f}")
print(f"  Turbulent peak intensity: {strehl_turb:.4f}")
print(f"  Strehl ratio: {strehl_turb/strehl_dl:.2%}")

Image Quality Comparison:
  Diffraction-limited peak intensity: 0.0460
  Turbulent peak intensity: 0.0034
  Strehl ratio: 7.41%

Even with \(r_0 = 20\) cm — which is moderate seeing — the Strehl ratio drops to just a few percent. This makes it clear why every major observatory needs adaptive optics to recover diffraction-limited performance.

The atmosphere isn’t static. Wind moves turbulent layers across the aperture, causing the PSF to evolve on millisecond timescales. We’ll simulate this frozen-flow evolution by advancing the InfiniteAtmosphericLayer in time and recording frames.

# Animation parameters
duration = 1.0  # seconds
fps = 40
num_frames = int(duration * fps)
time_step = duration / num_frames

writer = hcipy.FFMpegWriter('atmosphere_evolution.mp4', framerate=10)

print(f"Creating animation with {num_frames} frames...")

Creating animation with 40 frames...

Now we’ll loop through the time steps, evolving the atmosphere, computing the PSF, and recording frames to a video file. At each step we display the phase screen alongside the PSF with the current Strehl ratio overlaid.

for t in tqdm(np.arange(num_frames) * time_step):
    layer.evolve_until(t)

    wf_atmos = layer(wf)
    img = prop(wf_atmos)

    phase = layer.phase_for(wavelength)

    fig, axes = plt.subplots(1, 2, figsize=(14, 6))

    # Phase pattern
    im0 = hcipy.imshow_field(phase, cmap='RdBu_r', mask=aperture, vmin=-10, vmax=10, ax=axes[0])
    axes[0].set_title(f'Phase Screen (t = {t*1000:.0f} ms)', fontsize=12)
    axes[0].set_xlabel('x (m)')
    axes[0].set_ylabel('y (m)')
    plt.colorbar(im0, ax=axes[0], label='Phase (rad)')

    # Point spread function
    im1 = hcipy.imshow_field(np.log10(img.intensity / img.intensity.max()), vmin=-4, cmap='inferno', grid_units=spatial_resolution, ax=axes[1])
    axes[1].set_title('Point Spread Function', fontsize=12)
    axes[1].set_xlabel('x (lambda/D)')
    axes[1].set_ylabel('y (lambda/D)')
    plt.colorbar(im1, ax=axes[1], label='log10(Intensity)')

    # Strehl annotation
    strehl = np.max(img.intensity) / np.sum(img.intensity)
    fig.suptitle(f'Strehl Ratio: {strehl/strehl_dl:.1%}', fontsize=14, y=0.98)

    plt.tight_layout(rect=[0, 0, 1, 0.96])

    writer.add_frame(fig)
    plt.close(fig)

writer.close()
writer

100%|██████████| 40/40 [00:32<00:00,  1.25it/s]

The animation shows the phase screen drifting across the pupil while the PSF core jumps and speckles shift. Notice how the Strehl ratio fluctuates frame-to-frame — it can vary by a factor of 2 or more as different turbulent eddies cross the aperture.

In reality, atmospheric turbulence is concentrated in discrete layers at different altitudes — the boundary layer, the jet stream, and the tropopause — each with its own wind speed and direction. HCIPy provides make_mauna_kea_atmospheric_layers which returns a realistic profile based on measurements at Mauna Kea. We’ll randomize the wind directions of each layer to produce a more realistic, non-degenerate speckle pattern, then combine them into a MultiLayerAtmosphere.

# Mauna Kea atmospheric profile
layers = hcipy.make_mauna_kea_atmospheric_layers(pupil_grid, outer_scale=outer_scale)

# Randomize layer orientations
for layer in layers:
    angle = np.random.uniform(0, 2 * np.pi)
    velocity = np.linalg.norm(layer.velocity)
    layer.velocity = [np.cos(angle) * velocity, np.sin(angle) * velocity]

# Multi-layer atmosphere
atmos = hcipy.MultiLayerAtmosphere(layers, scintillation=True)

print(f"Mauna Kea Atmospheric Profile:")
print(f"  Number of layers: {len(layers)}")
print(f"\nLayer details:")
for i, layer in enumerate(layers):
    print(f"  Layer {i+1}: {layer.height:5.0f}m altitude, "
          f"Cn squared={layer.Cn_squared:.2e}, velocity={layer.velocity} m/s")

Mauna Kea Atmospheric Profile:
  Number of layers: 6

Layer details:
  Layer 1:   500m altitude, Cn squared=5.00e-14, velocity=[ 5.65557002 -8.24709208] m/s
  Layer 2:  1000m altitude, Cn squared=1.93e-14, velocity=[-9.35552189 -3.53188478] m/s
  Layer 3:  2000m altitude, Cn squared=1.46e-14, velocity=[ 9.84856669 -1.73370535] m/s
  Layer 4:  4000m altitude, Cn squared=3.19e-14, velocity=[-6.14066207 -7.89254517] m/s
  Layer 5:  8000m altitude, Cn squared=7.34e-14, velocity=[-5.5589905 -8.3124981] m/s
  Layer 6: 16000m altitude, Cn squared=2.96e-14, velocity=[ 9.51284607 -3.08314119] m/s

The Mauna Kea profile includes layers spanning altitudes from a few hundred meters to nearly 20 km, with the strongest turbulence in the boundary layer and the jet stream. Each layer contributes different amounts to the total wavefront error.

With a single thin layer, the atmosphere only imprints phase errors. But when light propagates through multiple layers at different altitudes, diffraction converts phase variations into amplitude variations — this is scintillation (the same phenomenon that makes stars twinkle). HCIPy’s MultiLayerAtmosphere handles this Fresnel propagation between layers when scintillation=True. The scintillation index — the normalized variance of the intensity — quantifies how strongly the atmosphere modulates the amplitude of the wavefront.

# Set turbulence strength
atmos.Cn_squared = hcipy.Cn_squared_from_fried_parameter(0.1, wavelength=550e-9)
atmos.reset()

# Input wavefront
wf_in = hcipy.Wavefront(pupil_grid.ones(), wavelength=wavelength)

We’ve set up a strong turbulence case (\(r_0 = 10\) cm) and a unit-amplitude input wavefront. Now let’s animate the evolution of both the amplitude (scintillation) and phase patterns over one second.

# Animation
duration = 1  # seconds
fps = 40
num_frames = int(duration * fps)
time_step = duration / num_frames

writer = hcipy.FFMpegWriter('scintillation_evolution.mp4', framerate=10)

print(f"Creating scintillation animation with {num_frames} frames...")

Creating scintillation animation with 40 frames...

The animation loop evolves the multi-layer atmosphere, computes the resulting wavefront, and records the intensity and phase patterns along with the scintillation index.

for t in tqdm(np.arange(num_frames) * time_step):
    atmos.evolve_until(t)
    wf_out = atmos(wf_in)

    fig, axes = plt.subplots(1, 2, figsize=(14, 6))

    # Intensity (scintillation)
    im0 = hcipy.imshow_field(wf_out.intensity * aperture, cmap='gray', ax=axes[0], vmin=0, vmax=3)
    axes[0].set_title(f'Intensity (t = {t*1000:.0f} ms)', fontsize=12)
    axes[0].set_xlabel('x (m)')
    axes[0].set_ylabel('y (m)')
    plt.colorbar(im0, ax=axes[0], label='Intensity')

    # Phase
    im1 = hcipy.imshow_field(wf_out.phase * aperture, cmap='RdBu_r', ax=axes[1], vmin=-3, vmax=3)
    axes[1].set_title('Phase', fontsize=12)
    axes[1].set_xlabel('x (m)')
    axes[1].set_ylabel('y (m)')
    plt.colorbar(im1, ax=axes[1], label='Phase (rad)')

    # Scintillation statistics
    intensity_in_aperture = wf_out.intensity[aperture > 0]
    scintillation_index = np.std(intensity_in_aperture) / np.mean(intensity_in_aperture)

    fig.suptitle(f'Scintillation Index: {scintillation_index:.3f}', fontsize=14, y=0.98)
    plt.tight_layout(rect=[0, 0, 1, 0.96])

    writer.add_frame(fig)
    plt.close(fig)

writer.close()
writer

100%|██████████| 40/40 [00:25<00:00,  1.56it/s]

The intensity across the pupil fluctuates significantly as different turbulent cells drift through the beam. A scintillation index above 0.1 indicates strong amplitude modulation, which conventional adaptive optics cannot correct with a deformable mirror alone (it requires a second DM for amplitude control).

We can also sweep across different seeing conditions without recreating the layers. The atmos.Cn_squared property scales the turbulence strength of all layers simultaneously. Let’s compare three representative seeing values: 1.0 arcsec (poor), 0.5 arcsec (median), and 0.25 arcsec (excellent).

# Reset atmosphere
atmos.reset()

# Seeing conditions
seeing_arcsec = [1.0, 0.5, 0.25]

We’ll create a figure with three subplots and populate each one with the PSF for a different seeing value, annotating each with the Strehl ratio.

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

for i, seeing in enumerate(seeing_arcsec):
    r0 = hcipy.seeing_to_fried_parameter(seeing)
    atmos.Cn_squared = hcipy.Cn_squared_from_fried_parameter(r0)
    atmos.reset()

    wf_out = atmos(wf_in)
    img = prop(wf_out)

    im = hcipy.imshow_field(
        np.log10(img.intensity / img.intensity.max()),
        vmin=-3, cmap='inferno', grid_units=spatial_resolution, ax=axes[i]
    )
    axes[i].set_title(f'r0 = {r0*100:.0f} cm at 500nm\n({seeing:.2f} arcsec seeing)', fontsize=11)
    axes[i].set_xlabel('x (lambda/D)')
    axes[i].set_ylabel('y (lambda/D)')

    # Strehl ratio
    strehl = np.max(img.intensity) / np.sum(img.intensity)
    axes[i].text(0.02, 0.98, f'Strehl: {strehl:.2%}',
                 transform=axes[i].transAxes, verticalalignment='top',
                 bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

plt.suptitle('Effect of Seeing on Image Quality', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

Better seeing (larger \(r_0\)) produces tighter PSFs with higher Strehl ratios — at 0.25 arcsec the Strehl is nearly 5× better than at 1.0 arcsec. Poor seeing smears the image into a broad halo with little distinguishable structure. The multi-layer model also produces more realistic speckle patterns than a single layer, with the characteristic broken spoke-like structure seen in real coronagraphic images.

# Cleanup
import os

if os.path.exists('atmosphere_evolution.mp4'):
    os.remove('atmosphere_evolution.mp4')
if os.path.exists('scintillation_evolution.mp4'):
    os.remove('scintillation_evolution.mp4')